Detecting isometries and symmetries of implicit algebraic surfaces

1.   Introduction

Computing equivalences and symmetries of geometric objects in a specific geometric setup, e.g., projective, affine, or Euclidean, plays a significant role in computer vision, computer graphics, pattern recognition, and computer-aided geometric design. Specifically, the studies addressing the problems regarding algebraic curves and surfaces have recently gained particular interest. Some of these studies are ^{[1,2,3,5,7,9,10,11,12,15,17]} for both rational and implicit curves, and ^{[4,6,8,9,18,19]} for both rational and implicit surfaces.

In this paper, similar to the above studies, we aim to address the problem of computing isometries between two implicit surfaces. An isometry is a transformation that preserves the metric properties of the underlying space. We refer the reader to ^[13] for an extensive account regarding isometries. We call two surfaces isometric or congruent if one of the surfaces is the image of the other under an isometry. Also, we call a surface symmetric if there exists an isometry, other than the identity, which leaves the surface invariant.

For rational surfaces, there are certain approaches ^[4,8,18] addressing the detection of isometries. However, to the best of our knowledge, there is only one algorithm ^[9] regarding the isometries of implicit surfaces. In ^[9], the authors present a comprehensive approach to compute projective and affine equivalences, starting from projective equivalences between finite sets of points, possibly given as roots of univariate polynomials. However, the approach in ^[9] leads to solving multivariate polynomial equations, something that we avoid here. Also, in a very recent paper ^[6], a novel approach concerning symmetries of implicit surfaces is provided. In ^[6], the authors present a complete algorithm to compute four types of notable symmetries, namely, rotational, axial, reflectional, and central symmetries. Nonetheless, the method given in ^[6] is not global, i.e., it requires a separate computation for each type of symmetry. On the contrary, when applied to just one surface, in which case we determine its symmetries, our method can determine all symmetries of a surface in one go.

The two main features of the method presented in this paper are the reduction to isometries fixing the origin, which allows us to transfer the problem to polynomials of smaller degree, namely some homogeneous forms of the resulting polynomials, and the use of polynomial factoring and gcd computation as an alternative to polynomial system solving. The reduction step employs a previous idea already provided in ^[6,7], where the problem is reduced to computing transformations fixing the origin. In some nongeneric bad cases, it is not possible to perform the reduction step; for those cases we provide a backup method inspired by ^[9]. The main algorithm is implemented in Maple and we provide extensive tests, all of which can be found publicly in the first author's website ^[16].

The paper is organized as follows. In Section 2, we recall some basic notions regarding isometries of surfaces, and present the statement of the problem. In Section 3, which comprises three subsections, we first explain the general strategy of the method, show how to reduce the problem to computing isometries fixing the origin, and provide a full algorithm. In Section 4, we provide results concerning extensive experimentation performed on the algorithm. If the reduction method fails, we provide a backup method in Section 5. Finally, we conclude the paper in Section 6.

2.   Preliminaries

Let ${\mathcal S}_f, {\mathcal S}_g\subset {\Bbb R}^3$ be two algebraic surfaces, implicitly defined by $f(x, y, z) = 0$, $g(x, y, z) = 0$. Writing ${{\boldsymbol{x}}} = (x, y, z)^\intercal$, an isometry of the space is a mapping ${\mathcal T}:{\Bbb R}^3\to {\Bbb R}^3$,

where $A$ is a $3\times 3$ orthogonal matrix, i.e., $AA^\intercal = I$, and ${{\boldsymbol{b}}}\in\mathbb{R}^3$. We will use the notation

so

We say that the surfaces ${\mathcal S}_f, {\mathcal S}_g$ are congruent, if there exists an isometry $\mathcal{T}$ such that ${\mathcal T}({\mathcal S}_f) = {\mathcal S}_g$. Under the hypothesis that the surfaces ${\mathcal S}_f, {\mathcal S}_g$ are irreducible, this definition can be transferred to the polynomials $f, g$ defining the surfaces, so that ${\mathcal S}_f, {\mathcal S}_g$ are congruent if, and only if,

where $\lambda\in {\Bbb R}-\{0\}$.

As a special case, we say that the surface ${\mathcal S}_f$ is symmetric if there is an isometry $\mathcal{T}$ other than the identity, leaving $\mathcal{S}_f$ invariant, i.e., $\mathcal{T}(\mathcal{S}_f) = \mathcal{S}_f$. This can be stated in terms of the polynomial $f$ as $f(\mathcal{T}({{\boldsymbol{x}}})) = \lambda f({{\boldsymbol{x}}})$. In this case, it can be shown that $\lambda = \pm 1$ ^[6].

In the rest of the paper, we will assume that none of the surfaces we are dealing with have multiple components and that they are neither cylindrical nor surfaces of revolution. Cylindrical surfaces and surfaces of revolution are the only surfaces admitting infinitely many symmetries ^[8], so under this last hypothesis, we can be sure that the number of symmetries between ${\mathcal S}_f, {\mathcal S}_g$ is finite: Indeed, if $\mathcal{T}_1$ and $\mathcal{T}_2$ are two isometries between ${\mathcal S}_f, {\mathcal S}_g$, then $\mathcal{T}_1\circ\mathcal{T}_2^{-1}$ is a symmetry of $\mathcal{S}_g$ and $\mathcal{T}_2^{-1}\circ\mathcal{T}_1$ is a symmetry of $\mathcal{S}_f$. Thus, the existence of infinitely many isometries between ${\mathcal S}_f, {\mathcal S}_g$ implies that ${\mathcal S}_f, {\mathcal S}_g$ have infinitely many symmetries themselves.

Since isometries preserve the degree, ${\mathcal S}_f$ and ${\mathcal S}_g$ must have the same degree. In our case, we will assume that $\mbox{deg}({\mathcal S}_f) = \mbox{deg}({\mathcal S}_g) = n\geq 3$. Notice that for $n = 1$ where the surfaces are planes and for $n = 2$ where the surfaces are quadrics, the problem can be solved easily using linear algebra. Thus, we can summarize the problem in the following way.

Problem: Given two implicit algebraic surfaces ${\mathcal S}_f, {\mathcal S}_g$, implicitly defined by polynomials $f(x, y, z), g(x, y, z)$ of the same degree $n\geq 3$ without multiple components, neither of them are cylindrical or a surface of revolution. Find an algorithm to check whether or not ${\mathcal S}_f, {\mathcal S}_g$ are congruent and to compute the isometries, if any, between them in the affirmative case.

3.   Detecting isometries between implicit surfaces

3.1.   General strategy

Our approach consists of the following main steps.

$(1)$ Reducing the problem to isometries fixing the origin. Isometries fixing the origin have the form $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$, so compared to Eq (2.1), they have three less parameters. In order to reduce our problem to this case, we need to locate two points ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$ such that $\mathcal{T}({{\boldsymbol{P}}}) = {{\boldsymbol{Q}}}$. By performing translations in each surface so that both ${{\boldsymbol{P}}}, {{\boldsymbol{Q}}}$ are taken to the origin, we move to $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$.

$(2)$ Eliminating $\lambda$. If $f(x, y, z)$ implicitly defines $\mathcal{S}_f$, where the degree of $f(x, y, z)$ is $n$, then we can write

where $f_d(x, y, z)$ are the {homogeneous forms} of $f$ of degree $d$, with $d = 0, 1, \ldots, n$; similarly for $g(x, y, z)$ and $\mathcal{S}_g$. If $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$ is a congruence between $\mathcal{S}_f, \mathcal{S}_g$, we have that $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$ is also a congruence between the surfaces defined by $f_d, g_d$, for all $d$, so

for all $d$. $\lambda$ can be written in terms of ${{\boldsymbol{x}}}$ and the entries of the matrix $A$ and, therefore, it is eliminated.

$(3)$ Elimination and factorization. We first determine a polynomial system in the original variables $x, y, z$, and the components $u, v, w$ of the congruence $\mathcal{T}$ using (2.3). After computing a Gröbner basis to eliminate, say, $u, v$, we can recover the functions $w = w(x, y, z)$ as the linear factors in $w$ of the resulting polynomial. For $u, v$, we proceed in a similar way.

In the first and third steps, we need to make use of two notions from Vector calculus, namely, the Laplacian and Gradient, which behave nicely under isometries. Recall that the laplacian operator $\Delta: = \dfrac{\partial^2}{\partial x^2}+\dfrac{\partial^2}{\partial y^2}+\dfrac{\partial^2}{\partial z^2}$ satisfies ^[2]

and the gradient operator $\nabla: = (\dfrac{\partial}{\partial x}, \dfrac{\partial}{\partial y}, \dfrac{\partial}{\partial z})$ satisfies ^[6]

where $\Vert \cdot\Vert$ denotes the Euclidean norm of a vector.

3.2.   Reduction of the problem to isometries fixing the origin

We start by reducing the problem to finding isometries fixing the origin. In order to do this, we form two other surfaces from the original surface, namely, the surface itself and the surfaces obtained by taking the Laplacian $\Delta f$ and the square of the norm $\Vert \nabla f\Vert^2$ of the gradient of the original surface. Let us denote the surfaces defined by $\Delta f$, $\Vert \nabla f\Vert^2$ and by $\mathcal{L}_f$, $\mathcal{N}_f$. Note that since $f(x, y, z)$ is, by hypothesis, not a plane, $\Vert \nabla f\Vert^2$ is not a constant. However, $\Delta f$ might be a constant, in which case we replace $\Delta f$ by the {Hessian} of $f$, $Hf$. An account regarding the invariance of the Hessian under isometries can be found in ^[7].

Now, using Eqs (3.3) and (3.4), we get that if two surfaces, $\mathcal{S}_f$ and $\mathcal{S}_g$, are congruent by an isometry $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}+{{\boldsymbol{b}}}$, then the surfaces $\mathcal{L}_f$, $\mathcal{N}_f$ and $\mathcal{L}_g$, $\mathcal{N}_g$ are also mutually congruent by the same isometry $\mathcal{T}$. Using this fact, we conclude that the intersection points of the systems $\{\mathcal{S}_f, \mathcal{L}_f, \mathcal{N}_f\}$ and $\{\mathcal{S}_g, \mathcal{L}_g, \mathcal{N}_g\}$ are congruent by $\mathcal{T}$. Assume that the set of intersection points of the system $\{\mathcal{S}_f, \mathcal{L}_f, \mathcal{N}_f\}$ is $\mathcal{P}: = \{P_i\}_{i = 1}^m$ and that the set of intersection points of the system $\{\mathcal{S}_g, \mathcal{L}_g, \mathcal{N}_g\}$ is $\mathcal{Q}: = \{Q_i\}_{i = 1}^m$. Then, since $\mathcal{T}$ is an isometry and $\mathcal{T}(\mathcal{P}) = \mathcal{Q}$, the barycenters of $\mathcal{P}$ and $\mathcal{Q}$ are also congruent by $\mathcal{T}$. Denote the barycenters of the intersection points of the above systems by ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$, respectively. Thus, we have the following result.

Lemma 3.1. If ${\mathcal S}_f, {\mathcal S}_g$ are congruent by ${\mathcal T}$, then ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$ are also congruent by ${\mathcal T}$.

Let us provide a method to determine ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$ efficiently. We explain the method only for ${{\boldsymbol{P}}}$, as the computation for ${{\boldsymbol{Q}}}$ is the same. First, we need to compute a Gröbner basis for the system $\{f, \Delta f, \Vert \nabla f\Vert^2\}$ using a lexicographic order with $z > y > x$ in order to eliminate the variables $y, z$. We denote the first element of the basis by $p_1(x) = \sum_{i = 0}^{m_1}\alpha_ix^i$. Since the first coordinates of the intersection points are the roots of $p_1(x)$, using Cardano-Vieta's formulae, we can determine the first coordinate of ${{\boldsymbol{P}}}$ as $-\dfrac{\alpha_{m_1-1}}{m_1\alpha_{m_1}}$. We then choose a bivariate polynomial $q(x, y)$ in $x, y$ in the basis and set $p_2(y) = Res_x(p_1(x), q(x, y)): = \sum_{i = 0}^{m_2}\beta_iy^i$. Similarly, the second coordinate of ${{\boldsymbol{P}}}$ can be found as $-\dfrac{\beta_{m_2-1}}{m_2\beta_{m_2}}$ and the last one can be found by a Gröbner basis of $\{f, \Delta f, \Vert \nabla f\Vert^2, p_1(x), p_2(y)\}$ using a lexicographic order with $x > y > z$. We denote the first element of the last basis by $p_3(z): = \sum_{i = 0}^{m_3}\gamma_iz^i$, then the third coordinate of ${{\boldsymbol{P}}}$ is $-\dfrac{\gamma_{m_3-1}}{m_3\gamma_{m_3}}$. Finally, we get

Remark 3.1. According to the extension Theorem ^[14], one needs to check whether or not the roots obtained by eliminating y and z can be extended. To do this, we need to extend these solutions to x, y and then to x, y, z. The key question is whether or not the leading coefficients vanish. If not, then the solution can be extended.

Thus, whenever the intersection of the systems $\{\mathcal{S}_f, \mathcal{L}_f, \mathcal{N}_f\}$ and $\{\mathcal{S}_g, \mathcal{L}_g, \mathcal{N}_g\}$ is nonempty and finite, we can determine the barycenters ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$ efficiently. If the intersection is either empty or not finite, we say that we fall in a $\texttt{FAIL}$ case. We will provide an alternative method for these bad cases in Section 5.

Now let us assume that we are not in a fail case. Once we find ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$, we apply translations ${{\boldsymbol{x}}}+{{\boldsymbol{P}}}$ and ${{\boldsymbol{x}}}+{{\boldsymbol{Q}}}$ to the surfaces $\mathcal{S}_f$ and $\mathcal{S}_g$ so that ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$ are moved to the origin. These translated surfaces are congruent by $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$. We will use the following notation for the translated implicit equations of the new surfaces:

Thus, we aim to find the isometries satisfying

where $(u, v, w) = (a_{11}x+a_{12}y+a_{13}z, a_{21}x+a_{22}y+a_{23}z, a_{31}x+a_{32}y+a_{33}z)$.

3.3.   Determining isometries fixing the origin

The polynomials $F, G$ in Eq (3.6) satisfy $F(\mathcal{T}({{\boldsymbol{x}}})) = \lambda G({{\boldsymbol{x}}})$, with $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$. Denoting $A{{\boldsymbol{x}}} = (u(x, y, z), v(x, y, z), w(x, y, z))$, we get $F(u, v, w) = \lambda G(x, y, z)$. We can write $F$ and $G$ in terms of their homogeneous forms as $F(x, y, z) = \sum_{i = 0}^nF_d(x, y, z)$ and $G(x, y, z) = \sum_{i = 0}^nG_d(x, y, z)$, where $F_d, G_d$ are homogeneous of degree $i$. In turn, we get $F_d(u, v, w) = \lambda G_d(x, y, z)$. Let $d_0$ be the smallest integer satisfying $0\leq d_0 \leq n$ such that $F_{d_0}$ and $G_{d_0}$ are not zero. Isolating $\lambda$ for $d_0$, we obtain $\lambda = \dfrac{F_{d_0}(u, v, w)}{G_{d_0}(x, y, z)}$. Substituting the latter in Eq (3.2) and clearing the denominators, we have equations like

where $0\leq i\leq n$. Let ${\mathcal R}$ be the polynomial system consisting of the equations of the type in Eq (3.8) with the smallest $i$, where the corresponding equations in Eq (3.8) are nonzero. We compute a Gröbner basis for ${\mathcal R}$ using the lexicographic order with $w > v > u$. Let us denote the first element of the basis by $\Phi_1$, which can be considered as a polynomial in $u$ with coefficients in $\mathbb{R}[x, y, z]$, then we have the following result, which essentially will follow from Eq (3.7) and the Elimination Theory.

Proposition 3.1. Let $\mathcal{S}_F$ and $\mathcal{S}_G$ be congruent by $(u, v, w) = (a_{11}x+a_{12}y+a_{13}z, a_{21}x+a_{22}y+a_{23}z, a_{31}x+a_{32}y+a_{33}z)$, then $\eta_1(u): = u-(a_{11}x+a_{12}y+a_{13}z)$ is a linear factor of $\Phi_1$.

Let $\Psi_1, \Psi_2$ be the result of substituting $u: = (a_{11}x+a_{12}y+a_{13}z)$ into two polynomials of ${\mathcal R}$; notice that $\Psi_1, \Psi_2$ are polynomials in $v, w, x, y, z$. Let

then we get the following result, which, again, follows from resultant properties and Eq (3.7).

Proposition 3.2. Let $\mathcal{S}_F$ and $\mathcal{S}_G$ be congruent by $(u, v, w) = (a_{11}x+a_{12}y+a_{13}z, a_{21}x+a_{22}y+a_{23}z, a_{31}x+a_{32}y+a_{33}z)$, then $\eta_2(v): = v-(a_{21}x+a_{22}y+a_{23}z)$ is a linear factor of $R(v)$.

Now, substituting the zeros of the linear factors $\eta_1$ and $\eta_2$ in $\Psi_1$ and $\Psi_2$, we get polynomials $\tilde{\Psi}_1, \tilde{\Psi}_2$ in $w, x, y, z$.

Proposition 3.3. Let $\mathcal{S}_F$ and $\mathcal{S}_G$ be congruent by $(u, v, w) = (a_{11}x+a_{12}y+a_{13}z, a_{21}x+a_{22}y+a_{23}z, a_{31}x+a_{32}y+a_{33}z)$, then $\eta_3(w): = w-(a_{31}x+a_{32}y+a_{33}z)$ is a linear factor of the common gcd of $\tilde{\Psi}_1$ and $\tilde{\Psi}_2$.

The following example illustrates these results.

Example 3.1. Consider the surfaces implicitly defined by the following polynomials.

The first surface is called a Cassini surface, whose plotting can be seen in Figure 1. The second surface is obtained by composing the first one with $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}+{{\boldsymbol{b}}}$ and multiplying the resulting polynomial by $9$, where

Applying our algorithm, we find the barycenters ${{\boldsymbol{P}}} = (0, 0, 0)$ and ${{\boldsymbol{Q}}} = (-\frac{1}{3}, \frac{4}{3}, \frac{1}{3})$. Translating $f$ and $g$ by ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$, we get

Forming the system ${\mathcal R}$ for $F$ and $G$ and computing a Gröbner basis, we have that

which admits two linear factors

Substituting the above factors into ${\mathcal R}$, we get $\Psi_1, \Psi_2$. Taking the resultant for each linear factor, we obtain $\eta_{1i}$, $i\in\{1, 2\}$

Substituting these factors into $\Psi_1, \Psi_2$, we get $\tilde{\Psi}_1, \tilde{\Psi}_2$. Taking their gcd, we obtain

All the above solutions together yields eight general solutions with $\lambda = \frac{1}{9}$

Note that while the computation of the barycenters took $0.063$ seconds, it took $0.125$ seconds to compute the whole example.

4.   Experimentation.

The algorithm IsoSurf was implemented in the computer algebra system Maple^{TM [20]} and was tested on a PC with a $2.4$ GHz Intel Core i$5$ processor and $8$ GB RAM. Technical details, examples and source codes of the procedures are provided in the first author's personal website ^[16].

In all of our tests and examples, in order to generate congruent surfaces, we take a surface $\mathcal{S}_f$ defined by the implicit equation $f(x, y, z) = 0$ and apply to $\mathcal{S}_f$ the isometry in Eq (3.10) (see Section 3.3), so that $\mathcal{S}_g$ is the surface implicitly defined by $g(x, y, z) = 0$, with $g = f\circ\mathcal{T}$.

We start by surfaces whose implicit equations are provided in Table 1.

Plots of the surfaces given in Table 1 are presented in Figure 2.

In Table 2, we provide the computation times for computing all isometries and symmetries of the surfaces given in Table 1. The table reveals that computing the symmetries is quite faster than computing the isometries. The reason is that the isometries make the first polynomial denser so that the algorithm deals with a more complicated system. Still, the algorithm handles isometries in a few seconds, as seen in the table. In addition to this, notice that we choose examples so that the surfaces admit lots of isometries/symmetries, including surfaces with $48$ isometries/symmetries.

We continue with computing symmetries of surfaces implicitly defined by random dense polynomials. To guarantee that the randomly generated surfaces have several symmetries, we apply the change of variables $(x, y, z)\to (x^2, y^2, z^2)$ to the random dense polynomial $f(x, y, z)$ defining the surface. The surfaces generated by this method have exactly $8$ symmetries

for all $i, j, k\in\{0, 1\}$, where $\mathcal{T}_{000}$ is the trivial symmetry. We test these surfaces according to various degrees and coefficient bitsizes. Recalling that the bitsize $\tau$ of an integer $m$ is defined to be the integer $\tau = \lceil log_2k\rceil+1$, we call the maximum bitsize of the coefficients of the polynomial defining a surface as the bitsize of the surface. In Table 3, we present the timings for computing symmetries of the randomly generated implicit surfaces with degrees ranging from $4$ to $12$ and bitsizes ranging from $4$ to $64$. In all of the cases, our algorithm can compute all $8$ symmetries in a couple of seconds.

5.   Backup method

In this section, we provide a method similar to the one presented in ^[9] for computing isometries whenever the method presented in Section 3 fails; this happens when the solution set of the systems $\{\mathcal{S}_f, \mathcal{L}_f, \mathcal{N}_f\}$ and $\{\mathcal{S}_g, \mathcal{L}_g, \mathcal{N}_g\}$ is either empty or not finite. In that case, we exploit the relationship between the homogeneous forms of maximum degree of the polynomials defining the surfaces we want to compare. We will consider a projective setup.

Let $\mathcal{S}_f$ and $\mathcal{S}_g$ be two surfaces implicitly defined by $f(x, y, z)$ and $g(x, y, z)$. Let $\tilde{\mathcal{S}}_F$ and $\tilde{\mathcal{S}}_G$ be their projective closures, implicitly defined by $F(x, y, z, t)$ and $G(x, y, z, t)$, where $t$ is the homogenization variable. If $\tilde{\mathcal{S}}_F$ and $\tilde{\mathcal{S}}_G$ are congruent, then

where $\bar{{{\boldsymbol{x}}}} = (x, y, z, t)^\intercal$ and

with ${{\boldsymbol{b}}}\in\mathbb{R}^3$ and $A$ is orthogonal. Collecting $F$ and $G$, with respect to the homogenization variable, we get

where $f_d$ and $g_d$ correspond to the homogeneous forms of degree $d$ of $f$ and $g$. Let us denote

Thus, if Eq (5.1) holds, then we have for all $0\leq d\leq n$,

For $d = n$, it is clear that $\tilde{f}_n(x, y, z) = \lambda g_n(x, y, z)$ implies that $f_n(A{{\boldsymbol{x}}}) = \lambda g_n({{\boldsymbol{x}}})$, so $f_n(x, y, z)$ and $g_n(x, y, z)$ are congruent by $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$, where ${{\boldsymbol{x}}} = (x, y, z)^\intercal$. Since $\mathcal{T}$ is an isometry fixing the origin, we can develop a method similar to the method given in Section 3.3. Denoting $\mathcal{T}({{\boldsymbol{x}}}) = (u, v, w)$, we get a system

Here, we find another equality in $\lambda$. Let us denote

for a polynomial $p$. Consider the polynomial $\Delta^k\Vert \nabla f_n\Vert^2$. Here, since $\Vert \nabla f_n\Vert^2$ is a polynomial of degree $2n-2$, we get, for all $k\geq n$, $\Delta^k\Vert \nabla f_n\Vert^2 = 0$. Let $k_0$ be the greatest integer satisfying $1\leq k_0\leq n-1$ such that $\Delta^{k_0}\Vert \nabla f_n\Vert^2$ is not a zero polynomial. Using the second equality in system (5.6), we get that

Substituting the above equality into the system (5.6) and clearing denominators, we get a system free from the parameter $\lambda$.

Finally, we proceed as in Section 3.3 but with replacing the system ${\mathcal R}$ by the above system. Once we determine the isometry $\mathcal{T}({{\boldsymbol{x}}}) = A{{\boldsymbol{x}}}$, we substitute the matrix $A$ in $M$ defined in Eq (5.2). Therefore, it remains to determine ${{\boldsymbol{b}}}$. To do this, we use the polynomial equality corresponding $d = n-1$ in Eq 5.5, which is a linear system. Let us provide an example below to illustrate the idea.

Example 5.1. Consider the surfaces implicitly defined by the polynomials

For these surfaces the intersections of the systems $\{\mathcal{S}_f, \mathcal{L}_f, \mathcal{N}_f\}$ and $\{\mathcal{S}_g, \mathcal{L}_g, \mathcal{N}_g\}$ are empty. Thus, the algorithm IsoSurf fails to reduce the problem using barycenters ${{\boldsymbol{P}}}$ and ${{\boldsymbol{Q}}}$. Applying the method given in this section gives exactly $8$ isometries between $\mathcal{S}_f$ and $\mathcal{S}_g$. These isometries are, with $\lambda = \frac{1}{9}$,

The whole computation to determine all isometries between $\mathcal{S}_f$ and $\mathcal{S}_g$ using the direct method took $0.187$ seconds.

6.   Conclusions

We have provided a new and computationally efficient algorithm for computing the isometries between two implicit algebraic surfaces, as well as computing symmetries when the input surfaces coincide. We have implemented the algorithm in the computer algebra system Maple to test it on various aspects. The results of the tests demonstrate that the algorithm works efficiently. Even for highly symmetric surfaces with higher degrees, the algorithm can compute all the isometries$/$symmetries in a few seconds.

The algorithm consisted of two main steps. The first step reduced the problem to finding isometries fixing the origin. Once we reduced the problem, the second step computed the isometries by determing their components as linear factors from polynomials originated from a Gröbner elimination process. If the reduction step failed, we used a backup alternative method, inspired by ^[9].

The reason why the method was efficient is that the reduction step allowed us to replace the original polynomials by polynomials of much smaller degrees, which are in fact homogeneous forms of the polynomials obtained after appropriately translating the original polynomials. The main computational tools that we employed were Gröbner basis computations, polynomial factoring and gcds. In particular, we completely avoided multivariate polynomial system solving.

One can wonder whether our ideas can be extended to detecting other equivalences between implicit surfaces, e.g., similarities, affine equivalences, or projective equivalences. For similarities, the behavior of the invariants used in this paper, namely, the gradient and the laplacian, were also good, so there might be some hope; however, for similarities, we have a new unknown, the similarity ratio, which needs to be considered. Affine or projective equivalences are more challenging, since the behavior of the invariants used in this paper were not good anymore.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank the reviewers for their valuable comments that helped us improve the earlier version of the paper. U$\breve{\text{g}}$ur Gözütok is supported by the grant $121C421$, in the scope of $2218$-National Postdoctoral Research Fellowship Program, from TUBITAK (The Scientific and Technological Research Council of Türkiye).

Conflict of interest

The authors declare that they have no conflict of interest.