You can reduce it to an algebraic problem as follows:
The definition of a surface of revolution is that there is some axis such that rotations about this axis leave the surface invariant. Let's denote the action of a rotation by angle $\theta$ about this axis by the map $\vec x\mapsto \vec R_\theta(\vec x)$. With your surface given as a level set $f(\vec{x})=0$, the invariance condition is that $f(\vec R_\theta(\vec x))=0$ whenever $f(\vec{x})=0$, for all $\theta$.
In particular, we can differentiate this with respect to $\theta$ to get $\vec\nabla f(\vec R_\theta(\vec x))\cdot \frac{\partial}{\partial\theta}\vec R_\theta(\vec x)=0$, which at $\theta=0$ gives $\vec\nabla f(\vec x)\cdot \vec k(\vec x)=0$, where $\vec k(\vec x)=\left.\frac{\partial}{\partial\theta}\vec R_\theta(\vec x)\right|_{\theta=0}$ is the vector field representing the action of an infinitesimal rotation. (If the language of differential geometry is familiar, this is a Killing vector field).
So with this language established, what we need to check is whether there is any Killing field $\vec k(\vec x)$ associated to a rotation, which is orthogonal to the gradient of $f$ everywhere on the surface (i.e., whenever $f(\vec x)=0$). In fact, this will be not just necessary, but also sufficient, since (intuitively) any rotation can be built from many small ones.
Luckily, it's quite straightforward to write down the most general Killing field: if $\vec x_0$ is a point on the axis of rotation, and $\vec a$ a unit vector pointing along the axis, we have $\vec k(\vec x)=\vec a \times (\vec x-\vec x_0)$. (Note that this is a degenerate parametrization, since we can pick any point on the axis, corresponding to shifting $\vec x_0$ by a multiple of $\vec a$, and also send $\vec a$ to $-\vec a$, to give the same rotation).
To summarize: the question has been recast as "are there any solutions $\vec x_0, \vec a$ to $\vec a \times (\vec x-\vec x_0)\cdot\vec\nabla f(\vec x)=0$, which hold for all $\vec x$ such that $f(\vec x)=0$?".
(You could also write the equation as $\det[\vec a, \, \vec x-\vec x_0, \,\vec\nabla f(\vec x)]=0$).
For your example, I got Mathematica to use this method. I let $\vec x_0=(x_0,y_0,0)$, taking $z_0=0$ to remove some degeneracy, and $\vec a=(a_x,a_y,a_z)$, giving 5 unknowns for the axis. I then found four points on the surface, by solving $f(x,y,z)=0$ with $y=z=0$ and $x=z=0$. Then I got it to solve $\vec a \times (\vec x-\vec x_0)\cdot\vec\nabla f(\vec x)=0$ for the four points, and $|\vec a|^2=1$ (5 equations for the 5 unknowns), getting a unique solution up to the sign of $\vec a$. I then substituted the solution into $\vec a \times (\vec x-\vec x_0)\cdot\vec\nabla f(\vec x)$ for general $\vec x$, getting zero, so the solution indeed gave an axis of revolution. It is simplified in this case because all the level sets of $f$ give a surface of revolution about the same axis, so this last step did not require assuming $f(\vec x)=0$.
TL;DR: Just add zeros.
The general quadric equation can be broken down into the sum of a quadratic form, a linear term and a constant: $$\mathbf x^T\mathtt A\mathbf x+2\mathbf b^T\mathbf x+c=0.$$ The coefficients can be packaged up into a single homogeneous matrix $\mathtt Q$ by setting $\hat{\mathbf x}=(\mathbf x^T; 1)^T$ so that the equation becomes $$\hat{\mathbf x}^T\mathtt Q\hat{\mathbf x} = \begin{bmatrix}\mathbf x&1\end{bmatrix} \begin{bmatrix}\mathtt A & \mathbf b \\ \mathbf b^T & c\end{bmatrix} \begin{bmatrix}\mathbf x\\1\end{bmatrix} = 0.$$ (This equation is homogeneous in the components of $\hat{\mathbf x}$, so we can relax the restriction that the last component be $1$.)
For your covariance matrix, the ellipsoid equation is $\mathbf x^T\mathtt M\mathbf x=1$, therefore $\mathbf b=0$ and $c=-1$, i.e., $$\mathtt Q = \begin{bmatrix}\mathtt M & \mathbf 0 \\ \mathbf 0^T & -1\end{bmatrix}.$$ The projection matrix $\mathtt P$ for orthogonal projection onto the $x$-$y$ plane simply deletes the third component of a vector, so the matrix $\mathtt C$ that represents the projected outline of the ellipsoid will also have this form, but with a $2\times2$ matrix $\mathtt M'$. Computing the matrix of the outline thus comes down to inverting the upper-left $2\times2$ submatrix of $\mathtt M^{-1}$ and the Cartesian equation of the resulting ellipse is $(x,y)\,\mathtt M'(x,y)^T=1$. If you’re so inclined, you could even develop a closed-form formula for the coefficients of this equation from $\mathtt M^{-1}=\operatorname{cof}(M)^T/\det(M)$. The resulting expressions turn out to be quite simple.
Best Answer
It took a while, but I have an answer.
This answer describes the transformation of a hyperbolic paraboloid into canonical form. The procedure is similar for other quadric shapes.
An affine transform $\mathbf{T}$ is of the form: \begin{equation} \mathbf{T} = \begin{bmatrix} T_{11} & T_{12} & T_{13} & T_{14} \\ T_{21} & T_{22} & T_{23} & T_{24} \\ T_{31} & T_{32} & T_{33} & T_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation}
Every quadric may be transformed via affine transformations into either a diagonal matrix or a near diagonal matrix, i.e. \begin{equation} \mathbf{X}^T \mathbf{T}^T \mathbf{A} \mathbf{T} \mathbf{X} = \mathbf{X}^T \mathbf{D} \mathbf{X} \end{equation} where $\mathbf{D}$ is a (near-)diagonal matrix. A near-diagonal matrix is symmetric and is of the form: \begin{equation} \mathbf{D} = \begin{bmatrix} D_{11} & 0 & 0 & 0 \\ 0 & D_{22} & 0 & 0 \\ 0 & 0 & 0 & D_{43} \\ 0 & 0 & D_{43} & 0 \end{bmatrix} \end{equation}
The canonical form of a hyperbolic paraboloid listed in Xu (Table 1, pp518) is not strictly in near-diagonal form: \begin{align} 0 &= z - xy \\ &= \mathbf{X}^T \mathbf{D'} \mathbf{X} \\ &= \mathbf{X}^T \begin{bmatrix} 0 & -0.5 & 0 & 0 \\ -0.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.5 \\ 0 & 0 & 0.5 & 0 \end{bmatrix} \mathbf{X} \end{align} but may be re-written in near-diagonal form as: \begin{align} 0 &= -x^2 + y^2 + z \\ &= \mathbf{X}^T \mathbf{D} \mathbf{X} \\ &= \mathbf{X}^T \begin{bmatrix} -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0.5 \\ 0 & 0 & 0.5 & 0 \end{bmatrix} \mathbf{X} \end{align}
The Xu form $\mathbf{D'}$ and the near-diagonal form $\mathbf{D}$ are related via an affine transform: \begin{align} \mathbf{D'} = \mathbf{M}^T \mathbf{D} \mathbf{M} \end{align} where: \begin{equation} \mathbf{M} = \begin{bmatrix} 0.5 & -0.5 & 0 & 0 \\ 0.5 & 0.5 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{equation}
Since the Xu's canonical form can be found from the near-diagonal form, our immediate task is to find the transformation $\mathbf{T}$ that casts a given hyperbolic paraboloid $\mathbf{A}$ into near-diagonal form. Although Dupont recommends determining the transform via Gaussian Reductions to avoid radicals, perhaps the easiest first step is to remove the rotation via an Eigendecomposition as suggested by Levin.
Consider the principle submatrix $\mathbf{A}_u$. Since it is symmetric, of rank 2 and has known signs for the eigenvalues, an eigendecomposition may be found such that \begin{align} \mathbf{A}_u = {\mathbf{R}_u}^T \mathbf{\Lambda} \mathbf{R}_u \end{align} where $\mathbf{R} \in SO(3)$ and $\mathbf{\Lambda}$ is in the form: \begin{equation} \mathbf{\Lambda} = \begin{bmatrix} \lambda_+ & 0 & 0 \\ 0 & \lambda_- & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} \end{equation} where $\lambda_+ > 0$ and $\lambda_- < 0$. Therefore, the first affine transform may be built as: \begin{align} \mathbf{T}_1 = \begin{bmatrix} \mathbf{R}_u & \mathbf{0} \\ \mathbf{0} & 1 \end{bmatrix} \end{align}
The second transform forces the eigenvalues to unit magnitude and it can be verified that: \begin{align} \mathbf{T}_2 = \begin{bmatrix} \frac{1}{\sqrt{\lambda_+}} & 0 & 0 & 0 \\ 0 & \frac{1}{\sqrt{\|\lambda_-\|}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \end{align}
Let $\mathbf{A}_2 = {\mathbf{T}_2}^T {\mathbf{T}_1}^T \mathbf{A} \mathbf{T}_1 \mathbf{T}_2$, which will be of the form: \begin{equation} \mathbf{A}_2 = \begin{bmatrix} 1 & 0 & 0 & a \\ 0 & -1 & 0 & b \\ 0 & 0 & 0 & c \\ a & b & c & d \end{bmatrix} \end{equation}
It can be verified that the transform required to cast into near-diagonal form is: \begin{equation} \mathbf{T}_3 = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ \frac{-a}{c} & \frac{-a}{b} & \frac{-1}{2c} & \frac{-d}{2c} \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation}
Therefore, the complete transformation required to cast the hyperbolic paraboloid into Xu's canonical form is: \begin{align} \mathbf{T} = \mathbf{T}_1 \mathbf{T}_2 \mathbf{T}_3 \mathbf{M} \end{align}