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$.
$\newcommand{\Reals}{\mathbf{R}}$Edit of January 12, 2017, to address clarifications from the comments.
If $S \subset \Reals^{3}$ is a smooth surface invariant under some one-parameter group $\Gamma$ of ambient (Euclidean) isometries, then "the geometry of $S$ is $\Gamma$-invariant". Specifically, scalar geometric quantities such as principal curvatures (or mean curvature, or Gaussian curvature) are constant along the orbits of $\Gamma$. Assuming at least one principal curvature is non-constant (i.e., that $S$ is not a portion of a sphere, plane, or cylinder), level curves of the principal curvature functions are $\Gamma$-invariant, i.e., lines, circles, or right circular helices.
To detect numerically whether or not a surface $S$ approximated by a discrete surface mesh is "$\Gamma$-invariant", one might calculate some scalar geometric quantity $K$, approximate some level curve of $K$, and then determine whether or not this level curve has constant curvature and torsion:
If yes, the curvature and torsion can be used to recover the (approximate) axis and pitch of the helical motion. If necessary, these data can be used to confirm the surface is "helical" with the stated axis and pitch.
If no, the surface $S$ is not helical.
The prospects of numerical success presumably depend on the fineness of the surface mesh and how closely the mesh approximates any helical structure the surface may have. It's entirely possible a more accurate and/or efficient strategy can be found with more analytic work.
(Original answer of December 24, 2016, with a few minor edits for attribution and notational consistency.)
I don't know of published literature, much less a standard name for these surfaces, but here are some observations:
Extrinsically, a self-sliding surface $S$ in Euclidean three-space is invariant under a one-parameter group $\Gamma$ of Euclidean isometries.
If $\Gamma$ consists of translations, $S$ is a (generalized) cylinder.
If $\Gamma$ consists of rotations, $S$ is a surface of rotation.
Otherwise, the orbits of $\Gamma$ are helices. In the (unpublished) past, a student (J. M. Antonio) and I called these "pasta surfaces".
Intrinsically, a self-sliding surface admits local Clairaut coordinates $(u, v)$ in which the metric has the form
$$
g = E(u)\, du^{2} + G(u)\, dv^{2},
$$
with the group action given by translation in $v$, and with $u$ a parameter for a "profile" (a curve transverse to the group orbits). (The group action foliates $S$ by curves. The field of orthogonals is one-dimensional, hence integrable: Through each point $p$, there exists a profile through $p$ that is everywhere perpendicular to the orbits. Let $u$ be a parameter for this curve, and let $v$ be induced by the group action.)
Thanks to the one-parameter family of isometries (which manifests analytically as $v$-independence of the metric components), the geometry of $S$ is "constant in $v$". Particularly, under an equivariant embedding (i.e., an embedding $i:S \hookrightarrow \mathbf{R}^{3}$ for which a one-parameter group of Euclidean isometries preserves the image and acts by intrinsic coordinate translation), the principal curvatures (and consequently the Gaussian and mean curvatures) are constant along $v$-coordinate curves.
It's also not difficult to show that one-parameter families of pasta surfaces of constant Gaussian curvature exist, analogous to the one-parameter families of surfaces of rotation having constant Gaussian curvature. (For studying Gaussian curvature, it's particularly convenient to arrange that $EG \equiv 1$, i.e., that $E(u) = 1/G(u)$. The Gaussian curvature in these "momentum" coordinates is linear in the metric components: $K = -\frac{1}{2}G''(u)$.) The only complete, embedded examples are spheres, cylinders, and planes, but examples based on the pseudosphere give isometric embeddings of arbitrarily large hyperbolic disks.
(I've never looked at constant mean curvature, but don't expect obstructions to existence.)
Edit of December 26, 2016: Up to a rigid motion of Euclidean three-space, every self-sliding surface (other than a generalized cylinder) arises as follows.
Fix a parametric curve $\gamma(u) = (r(u), 0, z(u))$ whose image lies in the half-plane $x > 0$, $y = 0$ and whose velocity is non-vanishing; a real number $k$ (the "pitch" or "slope" of a helix); and an arbitrary smooth function $\psi$ (the "gauge"); and put
$$
X(u, v) = \left[\begin{array}{@{}c@{}}
r(u) \cos(v + \psi(u)) \\
r(u) \sin(v + \psi(u)) \\
z(u) + k(v + \psi(u)) \\
\end{array}\right].
$$
The image of $X$ is the result of sliding the image of $\gamma$ (the green curve below) under the action
$$
(x, y, z) \mapsto (x\cos t - y\sin t, x\sin t + y\cos t, kt)
$$
of the additive group of real numbers. This action represents rotation about the $z$-axis if $k = 0$ and helical motion about the $z$-axis if $k \neq 0$. (Formally, translation along the $z$-axis arises by letting $k \to \infty$, but I assume you're mostly interested in helical motions.)
(The gauge function $\psi$ does not change the surface geometrically, but "deforms" the $u$-coordinate curves "internally" by sliding them along the $v$-coordinates curves (circles or helices). If $\psi \equiv 0$, the $u$-coordinates curves ($v$ constant) are plane curves, the images of $\gamma$ under the group action. For some purposes it's more convenient to choose $\psi$ so that the $u$- and $v$-coordinate curves are pointwise orthogonal, as in the diagram below.)
Best Answer
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.