Understanding ray against ellipsoid collision

So I have looked around for a better explanation for Ray against an Ellipsoid collision detection as well as its normal (on the collision point) calculation. I found this useful website and it has an implementation:

// ellipsoid centered at the origin with radii ra
vec2 eliIntersect( in vec3 ro, in vec3 rd, in vec3 ra )
{
    vec3 ocn = ro/ra;
    vec3 rdn = rd/ra;
    float a = dot( rdn, rdn );
    float b = dot( ocn, rdn );
    float c = dot( ocn, ocn );
    float h = b*b - a*(c-1.0);
    if( h<0.0 ) return vec2(-1.0); //no intersection
    h = sqrt(h);
    return vec2(-b-h,-b+h)/a;
}

vec3 eliNormal( in vec3 pos, in vec3 ra )
{
    return normalize( pos/(ra*ra) );
}

I just want to understand how this works and also how it would work if it's not centered at the origin?

Solution

An ellipsoid not only has an arbitrary center, but also arbitrary axis directions, pairwise orthogonal. I am discussing the general case.

You can deal with the ellipsoid center by subtracting the center coordinates from the vector ro.

When this is done, if necessary you can deal with the axis direction by forming a 3x3 matrix with the normalized axis vectors and appling the inverse matrix to the new ro and to rd. If the ellipsoid is axis-aligned, the matrix is unit and you don't have to care.

Now, we write the vector equation of the ray, ro + t.rd (with t≥0 for a half-ray) and plug it in the implicit equation of the ellipsoid,

x²/ra.x² + y²/ra.y² + z²/ra.z² = 1.

This gives

(ro.x + t rd.x)²/ra.x² + (ro.y + t rd.y)²/ra.y² + (ro.y + t rd.y)²/ra.z² = 1.

After expansion, which I will not detail, we get a quadratic equation in t. If the discriminant is negative, there are no real roots, so the ray misses the ellipsoid. Otherwise we get two intersection points, which should be on the positive side of the ray.

From t and the ray equation, you obtain the coordinates of an intersection. Don't forget to apply the matrix and add the ellipsoid center to return to the original coordinates.

Now for the normal vector, we compute the gradient of the implicit equation, giving the components x/ra.x², y/ra.y², z/ra.z². You can compute these from the x, y, z of the intersection. Finally, apply the matrix if necessary and normalize the vector.