For an ellipsoid of the form
$$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1 $$
with orientation vector $\vec r$ and centre at point $\vec p$, how to find whether a point $\vec q$ is inside the ellipsoid or not?
Note: The geometry is a oblate spheroid with a=b and therefore one axis is sufficient to define orientation.
So you have an ellipsoid whose equation in "nice" $xyz$-coordinates is $$ \frac{x^2+y^2}{a^2} + \frac{z^2}{c^2} = 1 $$ byt your test point $\vec q$ is given in "non-nice" $x'y'z'$-coordinates.
You also also know the $x'y'z'$-coordinates of the point $\vec p$ whose $xyz$-coordinates are $(0,0,0)$, and the $x'y'z'$-coordinates of a displacement vector $\vec r$ that is parallel to the "nice" $z$-axis.
Now start by computing $\vec q-\vec p$ and then find the nice $z$ coordinate by $$ z^2 = \left(\pm\frac{(\vec q-\vec p)\cdot \vec r}{|\vec r|}\right)^2 = \frac{((\vec q-\vec p)\cdot \vec r)^2}{\vec r\cdot \vec r} $$ This gives you half of what you need to insert into the equation of the ellipsoid. The other half we get from Pythagoras: $$ x^2+y^2 = |\vec q-\vec p|^2 - z^2 = (\vec q-\vec p)\cdot(\vec q-\vec p) - z^2 $$ and then it's just a matter of inserting into the equation and see whether the LHS is less than $1$ or not.
(If you need to do this many times, you can save some effort by pre-scaling $\vec r$ to have unit length, of course).
A part of this answer to a related question by the same poster applies to this question here as well, so I'll post it here in a suitably edited fashion.
First off, your question says knowing $\vec r$ is enough, so I assume that even your rotated spheroid is still centered at the origin. If this is not the case, the first step would be a translation to make it centered around the origin.
If you take the spheroid $\frac{x^2+y^2}{a^2}+\frac{z^2}{c^2}=1$ and rotate it so that the original $z$ axis aligns with a vector $\vec r=(s,t,u)$, what do you get? Assume for the moment that $r$ has unit length. Then instead of $z$ you write $(sx+ty+uz)$ since that is the dot product between $\vec r$ and $\vec v=(x,y,z)$ and the dot product gives the length of the orthogonal projection. For the projection onto the orthogonal complement of $\vec r$ you can use the cross product. That's because as the dot product is the cosine of the angle times the product of the lengths, the length of the cross product is the sine of the angle times the product of the lenths.
$$\vec r\times \vec v= \begin{pmatrix}s\\t\\u\end{pmatrix}\times \begin{pmatrix}x\\y\\z\end{pmatrix}= \begin{pmatrix}tz-uy\\ux-sz\\sy-tx\end{pmatrix}$$
This vector is perpendicular to $\vec r$ and $\vec v$, and its length equals that of the orthogonal projection of $\vec v$ into the orthogonal complement of $\vec r$. You can use its squared length instead of $x^2+y^2$ in your formula.
$$\frac{(tz-uy)^2+(ux-sz)^2+(sy-tx)^2}{a^2}+\frac{(sx+ty+uz)^2}{c^2}=1$$
Now drop the assumption that $\vec r$ has unit length. If it does not, then your numerators will be too big by a factor of $\lVert r\rVert^2$, so you need to divide by that factor.
$$\frac{(tz-uy)^2+(ux-sz)^2+(sy-tx)^2}{a^2(s^2+t^2+u^2)}+\frac{(sx+ty+uz)^2}{c^2(s^2+t^2+u^2)}=1$$
If you want a polynomial equation, i.e. want to avoid the divisions, then you can multiply by the common denominator.
\begin{multline*} c^2\bigl((tz-uy)^2+(ux-sz)^2+(sy-tx)^2\bigr) + a^2(sx+ty+uz)^2 \\ = a^2c^2(s^2+t^2+u^2) \end{multline*}
Comparing this to the currently accepted answer by Henning Makholm, the two are equivalent. My $\vec v$ is his $(\vec p-\vec q)$ since he considers that translation in all his formulas. Apart from that one can verify (e.g. by brute-force computation, preferrably using some CAS) that
$$\vec v\cdot\vec v-\frac{\left(\vec v\cdot\vec r\right)^2}{\vec r\cdot\vec r} = \frac{\left(\vec r\times\vec v\right)\cdot\left(\vec r\times\vec v\right)} {\vec r\cdot\vec r}$$
where the left hand side is what Henning used, and the right hand side is what I used. If you have $\lVert\vec r\rVert=\lVert\vec v\rVert=1$ this equation simply corresponds to $1-\cos^2\alpha=\sin^2\alpha$ where $\alpha$ is the angle between $\vec r$ and $\vec v$. I know that at least for me it was easier to actually write the equation down at the coordinate level using the right hand side, but that might be a matter of practice.
By applying the reverse rotation and translation you will find the coordinates of $q$ in the ellipsoid-centered coordinates. Then check that the LHS $\le1$.
