Your intuition of setting the two equations equal is correct and that is how you solve for the intersection. I'll provide a full explanation, with code examples.
A common way of representing a plane $\mathbf{P}$ is in point-normal form, $\mathbf{\vec{n}} \cdot (X-Y)=0$, where $\mathbf{\vec{n}}$ is the plane normal and both $X$ and $Y$ are points that lie in the plane. This can be rewritten into constant-normal form by distributing the dot product and rearranging terms to obtain: $\mathbf{\vec{n}} \cdot X = d$, where $d = \mathbf{\vec{n}} \cdot Y$ which is equal to the distance from the origin when $\mathbf{\vec{n}}$ is unit-length. Below is a simple data structure that you might use to represent a plane, and the signature of a constructor that will compute the plane from three points in $\mathbb{R^3}$. Implementation is left as an exercise to the reader ;).
struct Plane {
Vector3 n; // normal
float d; // distance from origin
Plane(); // default constructor
Plane(Vector3 a, Vector3 b, Vector3 c); // plane from 3 points
Vector3 intersectLine(Vector3 a, Vector3 b); // we'll get to this later
};
Given two points, $A$ and $B$, a line can be represented parametrically by adding to one point the vector formed by the two points, scaled by a parameter $t$. In symbols, $L(t) = A + t(B-A)$. Using your intuition, we insert this equation (whose output is a point), into $X$ in the constant-normal plane representation: $\mathbf{\vec{n}} \cdot [A + t(B-A)] = d$. We want to know how many copies of $(B-A)$ we need to add to $A$ to get to a point that lies within the plane, in other words we want to solve for $t$. Doing some fancy algebra, we obtain: $t = \frac{d-\mathbf{\vec{n}} \cdot A}{\mathbf{\vec{n}} \cdot (B-A)}$. We can (finally) stick this expression for $t$ back into the equation for our line to obtain: $I = A+\frac{d - (\mathbf{\vec{n}} \cdot A)}{\mathbf{\vec{n}} \cdot (B-A)}(B-A).$
Armed with this equation, we can now implement a nice function that will tell what we want to know:
Vector3 Plane::intersectLine(Vector3 a, Vector3 b) {
Vector3 ba = b-a;
float nDotA = Vector3::dotProduct(n, a);
float nDotBA = Vector3::dotProduct(n, ba);
return a + (((d - nDotA)/nDotBA) * ba);
}
Hopefully this works for you, and hopefully I didn't fudge any of the details! If you plan to be doing a lot of this sort of geometric computing it's worthwhile to pick up Christer Ericson's Real-time Collision Detection, which is an excellent reference source for this sort of thing. Alternatively, you could snag some already-constructed classes from something like OGRE3D, if you're not particularly interested in creating your own.
[Collecting the discussion from the comments into an answer.]
You can certainly expand the expression that you’ve got and collect terms to get a quadratic equation in $\lambda$: Setting $\mathbf p = \mathbf o-\mathbf c$ and $\mathbf Q=\mathbf R^T\mathbf A\mathbf R$, you have $$(\lambda\mathbf l+\mathbf p)^T\mathbf Q(\lambda\mathbf l+\mathbf p) = (\mathbf l^T\mathbf Q\mathbf l)\lambda^2+2(\mathbf l^T\mathbf Q\mathbf p)\lambda +\mathbf p^T\mathbf Q\mathbf p=1,$$ which has real solutions when $(\mathbf l^T\mathbf Q\mathbf p)^2\ge(\mathbf l^T\mathbf Q\mathbf l)(\mathbf p^T\mathbf Q\mathbf p-1)$.
On the other hand, you’re essentially constructing the ellipsoid by scaling, rotating and translating the unit sphere. Since you’re not interested in the actual points of intersection, you could apply the inverse transformation to the line and then use the point-line distance formula to determine the line’s distance from the origin. If this is no greater than one, then there is at least one intersection.
If $\mathbf x_1$ and $\mathbf x_2$ are two points on a line, then the distance of the point $\mathbf x_0$ from this line is given by the formula (from Wolfram MathWorld) $${\lvert (\mathbf x_1-\mathbf x_0)\times(\mathbf x_2-\mathbf x_0)\rvert \over \lvert\mathbf x_2-\mathbf x_1\rvert}.$$ If $\mathbf x_0$ is the origin, this reduces to $${\lvert \mathbf x_2\times\mathbf x_1\rvert \over \lvert \mathbf x_2-\mathbf x_1\rvert}.\tag1$$ Now, since in your question you have $\mathbf R^T\mathbf A\mathbf R$ for the matrix of the ellipse, the matrix $\mathbf R$ already represents the inverse rotation. That is, it aligns the ellipsoid with the coordinate axes. (If instead you meant for it to represent a rotation that maps the standard basis onto the principal axes, then the ellipsoid’s matrix should be $\mathbf R\mathbf A\mathbf R^T$ instead and the inverse map below should use $\mathbf R^{-1}=\mathbf R^T$ instead of $\mathbf R$.) Setting $\mathbf S=\mathbf A^{1/2} = \operatorname{diag}(1/a,1/b,1/c)$, the required inverse map is therefore $\mathbf x\mapsto \mathbf S\mathbf R(\mathbf x-\mathbf c)$. A convenient choice for the two points on the line are the images of $\mathbf o$ and $\mathbf o+\mathbf l$. Plugging them into (1), we get $${\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times \mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert \over \lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)-\mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert} = {\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert \over \lvert \mathbf S\mathbf R\,\mathbf l\rvert}.$$ Since we’re comparing a ratio of nonnegative numbers to $1$, we can square and rearrange to get $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert^2 \le \lvert \mathbf S\mathbf R\,\mathbf l\rvert^2.$$ This looks like it could be more efficient to compute that the quadratic equation’s discriminant.
We can even go a bit farther. Using $\lvert\mathbf u\times\mathbf v\rvert^2=\lvert\mathbf u\rvert^2\lvert\mathbf v\rvert^2-(\mathbf u\cdot\mathbf v)^2$, $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l\rvert^2=\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2\lvert\mathbf S\mathbf R\,\mathbf l\rvert^2-\left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2$$ so $$\left(\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2-1\right) \lvert\mathbf S\mathbf R\,\mathbf l\rvert^2 \le \left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2,$$ but $\mathbf Q=(\mathbf S\mathbf R)^T(\mathbf S\mathbf R)$, so this is just the original inequality derived from the quadratic equation, refactored.
Note that if you already have unit direction vectors for the principal axes, there’s no reason for the purposes of this calculation to decompose the rotation into three primitive rotations as in your question: $\mathbf R$ is simply the matrix with these vectors as its rows. Note also that since $\mathbf S$ is diagonal, you can implement multiplication of a vector by $\mathbf S$ with three multiplications.
Best Answer
Translate the ellipsoid to origin, by subtracting $(x_3, y_3, z_3)$ from $P1$ and $P2$, so the line is parametrised as $$\vec{p} = \vec{v}_0 + t \vec{v}_1$$ where $$\begin{aligned} \vec{v}_0 = (\chi_0, \gamma_0, \zeta_0) &= (x_1 - x_3, y_1 - y_3, z_1 - z_3) \\ \vec{v}_1 = (\chi_1, \gamma_1, \zeta_1) &= (x_2 - x_1, y_2 - y_1, z_2 - z_1) \\ \end{aligned}$$ and the ellipsoid is $$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{b^2} = 1$$ Substituting $\vec{p}$ into the ellipsoid yields $$\left(\frac{\chi_0 + t \chi_1}{a}\right)^2 + \left(\frac{\gamma_0 + t \gamma_1}{b}\right)^2 + \left(\frac{\zeta_0 + t \zeta_1}{c}\right)^2 = 1$$ Expand the terms and collect the coefficients for powers of $t$ and you get $$\begin{aligned} \left( \displaystyle \frac{\chi_1^2}{a^2} + \frac{\gamma_1^2}{b^2} + \frac{\zeta_1^2}{c^2} \right) & t^2 ~ + \\ \left( \displaystyle \frac{2 \chi_0 \chi_1}{a^2} + \frac{2 \gamma_0 \gamma_1}{b^2} + \frac{2 \zeta_0 \zeta_1}{c^2} \right) & t ~ + \\ \left( \frac{\chi_0^2}{a^2} + \frac{\gamma_0^2}{b^2} + \frac{\zeta_0^2}{c^2} - 1 \right) & ~ = 0 \\ \end{aligned}$$ This is a simple quadratic equation in $t$, which can have 0, 1, or 2 real roots $t$.
Remember that we only translated the coordinate system so that they were relative to the center of the ellipsoid, but we didn't scale or rotate the coordinate system. So, after you have found $t$, you can find the point it refers to, in the original coordinate system, using your original parametrised line, $$P = P1 + t ( P2 - P1 )$$ i.e. $$\left\lbrace ~ \begin{aligned} x &= x_1 + t \chi_1 \\ y &= y_1 + t \gamma_1 \\ z &= z_1 + t \zeta_1 \\ \end{aligned} \right.$$ for each intersection point $t$.