I think you could try something like the following.
Let an offset ellipsoid be:
E(u,v) = [a * cos(u) * sin(v); b * sin(u) * sin(v); c * cos(u)] + [x_off; y_off; z_off]
where
u and v are like spherical coordinate angles.
a, b, and c are ellipsoid deformation parameters.
x_off, y_off, and z_off are the offsets defining the center of the ellipsoid.
Let the unit sphere at the origin be:
S(theta,phi) = [cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi)]
where
theta and phi are spherical coordinate angles.
Then you want to solve something like:
E(u,v) = g(u,v) * S'
=> S' = E/g
where
g(u,v) scales the vector E(u,v) (which gives each point on the ellipsoid E) back to a vector S' (a point on the surface of the sphere S).
S' is the subset of points on the surface of the unit sphere S that are "occluded" by E. This means |E/g| = 1 for all (u,v).
So basically you have the surface of the ellipsoid E parameterized by (u,v), which gives a vector in 3-space as E(u,v), and the surface of the unit sphere S parameterized by (theta,phi), which gives a vector in 3-space as S(theta,phi). You then want to find the subset of points S' in S that can be scaled by a real number in order to produce a vector on the surface E. We have to use a different real number, in general, for each surface point E to it's map on S because the distance changes across the curved surfaces. Therefore, we use the real function g(u,v) over the parameterization of the ellipsoid E.
This might give some strange results if the ellipsoid is intersecting with the unit sphere, so the above simply assumes that |E| > |S| for all u and v. There is still the tricky part of solving for g(u,v), but this is the basic setup.