WLOG I assume that $z_2=0$.
Let a ray be cast to the viewing plane in an arbitrary direction $\vec d=(x,y,z_1)$, so that a point along it is $t\vec d$. If we plug the coordinates to the sphere equation,
$$(t\vec d-\vec c)^2=r^2,$$
we get a second degree equation in $t$
$$\vec d^2t^2-2\vec d\cdot\vec c t+\vec c^2-r^2=0.$$
The ray passes through the apparent outline when the equation has a double root, i.e. when
$$(\vec d\cdot\vec c)^2-\vec d^2(\vec c^2-r^2)=0,$$
or
$$(xx_0+yy_0+z_1z_0)^2-(x^2+y^2+z_1^2)(x_0^2+y_0^2+z_0^2-r^2)=0.$$
This is the implicit equation of the searched ellipse. You can reduce it to the canonical form by well-kown techniques of translation and rotation (also).
The numerical result that you got is (almost) correct, but it looks like you might not be interpreting it correctly. Szeliski glosses over this stuff rather quickly.
First, observe that the projection that you’ve constructed in your test leaves a point’s $z$-coordinate unchanged, so the inverse map must do the same. You should be getting something like $(-5.178,-4.144,10.,1.)$ for the inverse image of your test point. Looks like the translation by $1$ ended up in the wrong place in your world-to-camera matrix. Even so, the correct $y$-coordinate in the world frame is also negative, so that still needs explanation.
The image coordinate system used with these matrices is also homogeneous, so that after applying the projection the result still needs to be normalized by dividing through by the $z$-coordinate. Thus, the test point $(200,200,10,1)$ actually corresponds to the point $(20,20)$ in the image. According to the calibration matrix, the camera’s axis is at about $(166,108)$ in the image, so after translating back this test point is firmly in the third quadrant. The rest of the inverse map is scaling and a small translation in $y$, so it’s not really surprising that both the $x$ and $y$ world coordinates end up being negative.
Lastly, the inverse map gives you a point on the ray, but to get the direction of the ray you have to subtract the camera’s world coordinate position—the starting point of the ray—which gives $(-5.178,-3.144,10.)$ for the direction vector of the ray in the world frame. If you had applied the inverse map to the homogeneous direction vector $(200,200,10,0)$ instead of the point $(200,200,10,1)$ you would’ve gotten this value directly. As well, if you back-map the equivalent image point $(20,20,1,1)$ and compute the resulting direction vector (or map the vector $(20,20,1,0)$ directly), you’ll end up with a scalar multiple of the above direction vector.
You don’t need to do anything special to incorporate a mapped point’s discrepancy, which is just the reciprocal of its (signed) distance from the camera. The discrepancy shows up as the $w$-coordinate of the projected point after normalization. Using the same test point of $(200,200,10,1)$, its normalized equivalent is $(20,20,1,1/10)$, so $d=0.10$ for this point. Back-mapping this gives $(-0.508,-0.414,1.0,0.1)$ and normalizing by dividing by $w$ produces the same answer as at top.
Best Answer
You can construct a line through the center of the sphere if you know the center of the projection and if you can identify the axes of the ellipse.
Let $O$ be the center of the projection and let $P$ and $Q$ be the endpoints of the major axis of the ellipse. Construct the lines $OP$ and $OQ$. These lines lie along the surface of the cone. The plane containing these lines also contains the axis of the cone. Bisect the angle $\angle POQ$. Then the angle bisector is the axis of the cone and passes through both the center of the projection and the center of the sphere.
The angle bisector of $\angle POQ$ intersects the major axis of the ellipse at a point $M$ that divides the major axis into two segments whose lengths have the ratio $\lvert MP\rvert:\lvert MQ\rvert = \lvert OP\rvert:\lvert OQ\rvert$. This implies that $M$ is the center of the ellipse only if $\lvert OP\rvert = \lvert OQ\rvert$, in which case the ellipse is a circle.
If you know the radius and location of the sphere then it is possible to find the point $M$ using that information.
If you know only the ellipse then it is not possible to identify the point $M$ where the axis of the cone intersects the plane of the ellipse. This is because there are infinitely many points that could be the center of a projection mapping a sphere to the given ellipse. See the answers to From ellipse equation to circular cone axis.