How to determine that a surface is symmetric

Here is an answer to this question in the case of compact surfaces (without boundary); maybe these ideas can be used in the general case as well. Let $F$ be a compact surface in $R^3$, bounding a solid $S$. In the setting you are interested in, one will probably have $F=\{x: f(x)=c\}$ and $S=\{x: f(x)\le c\}$. I will also assume that $f$ is a polynomial (I do not think it is really necessary, but I use this assumption to detect surfaces of revolution).

The solid $S$ then has a Loewner-Jones ellipsoid, which is the unique ellipsoid $E$ of the least volume containing $F$ (in fact, every compact in $R^n$ with not contained in a hyperplane has such ellipsoid). This ellipsoid is also known as Jones ellipsoid and least volume ellipsoid. It was discovered by Loewner and then rediscovered by Jones, if I remember correctly. People also use the largest volume ellipsoid contained in $S$, that one will also work.

The key is that, in view of uniqueness of $E$, every symmetry of $F$ is also a symmetry of $E$. There is a vast literature in the computational math community describing various algorithms for finding $E$, here is just a random paper on this subject which I found by googling. Now, suppose you are lucky and your ellipsoid has distinct (beyond the margin of error in the algorithm you would be using) principal axes. Then $E$ has exactly 3 planes of symmetry which pass through its center $C$ and pairs of axes. Then you test if symmetries in these planes preserve $S$.

Now, suppose you are unlucky and $E$ has a rotational symmetry in its axis $A$ (but is not the sphere). By intersecting $S$ with the plane $P$ orthogonal to $A$, we now get a 2-dimensional problem: Given a (closed and bounded) curve $\Gamma$ in the plane and a center $C$, determine if $\Gamma$ has a reflection symmetry in a line $L$ (through $C$). This is relatively easy since each point of intersection of $\Gamma$ with $L$ is a critical point of the distance function from $C$, so you can find such points using Lagrange multipliers and then check if any of the corresponding symmetries preserve $F$. This method will fail if $\Gamma$ is a circle. Instead of intersecting $F$ with planes passing through $C$, we can use other planes $P_k$ orthogonal to $L$. Now, if too many (comparing to the degree of $f$) of these intersections with $F$ are circles, then $F$ is a surface of revolution with the axis $A$. (I think, it should follow from Bezout theorem, but this has to be checked.) In this case, of course, any plane through $A$ will work.

The remaining case when $E$ is a sphere is similar to the circular case in one dimension lower, you would be then looking for critical points of the distance function from the origin.

Remark: The argument with Lowener-Jones ellipsoid should also handle your second question, about surfaces of revolution, in the bounded case.

Here's a stab at it. You can translate your question into whether or not a system of certain polynomial equations has a real solution.

A generic rotation can be described by choosing an axis using two angles in spherical coordinates, and then choosing a third angle to rotate by. The matrix for such a rotation has entries that are quartic polynomials in six variables that represent the sines and cosines of the three angles.

If your surface has a plane of symmetry, then for the right choice of these three angles, followed by a translation in the $z$-axis direction, you should be left with a polynomial that has only even powers of $z$.

If you give names to the seven unknowns (six trig values and one translation value), and apply the generic rotation and translation to the defining equation of the surface, then we can isolate the coefficients of $z,xz,yz,x^2z,z^3,z^5,$ etc. and see if we can solve for them all simultaneously being $0$. We already have three quadratic relations on the six trig variables. If we get equations in the seven unknowns corresponding to $z,xz,yz,x^2z,z^3,z^5$, etc (there will only be finitely many) then we have a system of polynomial equations in seven unknowns.

Such a system is consistent/inconsistent in $\mathbb{C}$ according to whether certain algebraic conditions are met on the equations by Hilbert's nullstellensatz. I think a CAS would be able to determine the nullstellensatz criterion, but I would consult an expert. I'm not sure what can be done about $\mathbb{R}$.

In 2D with your example, it would go like this, using $c=\cos\theta$, $s=\sin\theta$, and $a$ for the translation variable in the $y$-direction:

$$27 x^3 + 108 x^2 y + 144 x y^2 + 64 y^3 - 80 x^2 + 120 x y - 45 y^2 - 200 x + 150 y - 125 = 0$$

Apply rotation, followed by translation:

$$27(xc+(y+a)s)^3 +108 (xc+(y+a)s)^2(xs-(y+a)c) +144(xc+(y+a)s)(xs-(y+a)c)^2 +64(xs-(y+a)c)^3 -80(xc+(y+a)s)^2 +120(xc+(y+a)s)(xs-(y+a)c) -45(xs-(y+a)c)^2 -200(xc+(y+a)s) +150(xs-(y+a)c)-125 = 0$$

I cheat and have a CAS expand:

$$27 a^3 s^3+81 a^2 c s^2 x+81 a^2 s^3 y+81 a c^2 s x^2+162 a c s^2 x y+81 a s^3 y^2+27 c^3 x^3+81 c^2 s x^2 y+81 c s^2 x y^2+27 s^3 y^3 -108 a^3 c s^2-216 a^2 c^2 s x-324 a^2 c s^2 y+108 a^2 s^3 x-108 a c^3 x^2-432 a c^2 s x y+216 a c s^2 x^2-324 a c s^2 y^2+216 a s^3 x y-108 c^3 x^2 y+108 c^2 s x^3-216 c^2 s x y^2+216 c s^2 x^2 y-108 c s^2 y^3+108 s^3 x y^2 +144 a^3 c^2 s+144 a^2 c^3 x+432 a^2 c^2 s y-288 a^2 c s^2 x+288 a c^3 x y-288 a c^2 s x^2+432 a c^2 s y^2-576 a c s^2 x y+144 a s^3 x^2+144 c^3 x y^2-288 c^2 s x^2 y+144 c^2 s y^3+144 c s^2 x^3-288 c s^2 x y^2+144 s^3 x^2 y -64 a^3 c^3-192 a^2 c^3 y+192 a^2 c^2 s x-192 a c^3 y^2+384 a c^2 s x y-192 a c s^2 x^2-64 c^3 y^3+192 c^2 s x y^2-192 c s^2 x^2 y+64 s^3 x^3 -80 a^2 s^2-160 a c s x-160 a s^2 y-80 c^2 x^2-160 c s x y-80 s^2 y^2 -120 a^2 c s-120 a c^2 x-240 a c s y+120 a s^2 x-120 c^2 x y+120 c s x^2-120 c s y^2+120 s^2 x y -45 a^2 c^2-90 a c^2 y+90 a c s x-45 c^2 y^2+90 c s x y-45 s^2 x^2 -200 a s-200 c x-200 s y +150xs-150yc-150ac-125=0 $$

Ignoring terms with even powers of $y$:

$$81 a^2 s^3 y+162 a c s^2 x y+81 c^2 s x^2 y+27 s^3 y^3 -324 a^2 c s^2 y-432 a c^2 s x y+216 a s^3 x y-108 c^3 x^2 y+216 c s^2 x^2 y-108 c s^2 y^3 +432 a^2 c^2 s y+288 a c^3 x y-576 a c s^2 x y-288 c^2 s x^2 y+144 c^2 s y^3+144 s^3 x^2 y -192 a^2 c^3 y+384 a c^2 s x y-64 c^3 y^3-192 c s^2 x^2 y -160 a s^2 y-160 c s x y -240 a c s y-120 c^2 x y+120 s^2 x y -90 a c^2 y+90 c s x y -200 s y -150yc=0 $$

Starting with the relation between $c$ and $s$, and then grouping coefficients of $y$, $xy$, etc setting them equal to $0$:

$$\begin{align} c^2+s^2&=1\\ 81 a^2 s^3 -324 a^2 c s^2 +432 a^2 c^2 s-192 a^2 c^3-160 a s^2-240 a c s-90 a c^2-200 s -150c&=0&(y)\\ 162 a c s^2 -432 a c^2 s +216 a s^3 +288 a c^3-576 a c s^2+384 a c^2 s-160 c s-120 c^2+120 s^2+90 c s&=0&(xy)\\ 81 c^2 s -108 c^3 +216 c s^2-288 c^2 s+144 s^3-192 c s^2&=0&(x^2 y)\\ 27 s^3-108 c s^2+144 c^2 s-64 c^3&=0&( y^3) \end{align}$$

And now the question is equivalent to asking if this system of polynomials in three variables has a solution. Again, I think this is something that a CAS can do, but I'd consult an expert first.