This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 1

Finding the axis of revolution of an algebraic surface of revolution ´ Juan G. Alcazar and Ron Goldman Abstract—We present an algorithm for extracting the axis of revolution from the implicit equation of an algebraic surface of revolution based on three distinct computational methods: factoring the highest order form into quadrics, contracting the tensor of the highest order form, and using univariate resultants and gcds. We compare and contrast the advantages and disadvantages of each of these three techniques and we derive conditions under which each technique is most appropriate. In addition, we provide several necessary conditions for an implicit algebraic equation to represent a surface of revolution. Index Terms—Algebraic surfaces, surfaces of revolution, axis of revolution.

F

1

I NTRODUCTION

A

SURFACE OF REVOLUTION is a surface generated by rotating a space curve, the directrix or profile curve, about a straight line, the axis of revolution. Surfaces of revolution are ubiquitous both in theoretical mathematics and in computational science and engineering. Spheres, cylinders, cones, and tori, are some of the simplest surfaces of revolution, which appear both in elementary mathematics and in common human artifacts. Vases, pipes, wheels, architectural columns, cups and plates can all be constructed from surfaces of revolution. In Geometric Modeling, surfaces of revolution are typically represented either implicitly, i.e. as the set of zeroes of a polynomial F (x, y, z), or parametrically, i.e. by means of a rational parametrization (x(t, s), y(t, s), z(t, s)). Starting from one of these representations, a number of problems concerning surfaces of revolution have been considered in the literature: implicitizing [3], [17], [18], intersecting [8], [12], offsetting [16] and reparametrizing [1]. The problem of finding a surface of revolution approximating a discrete set of points has also been addressed in [14], [15], [22]. Furthermore, some special properties of quadrics of revolution are investigated in [10]. In easy cases, given the implicit representation of an algebraic surface of revolution the axis of revolution may be self-evident. For example, it is not difficult to guess that z − (x2 + y 2 ) = 0 represents a circular paraboloid, resulting from rotating the parabola {z = x2 , y = 0} around the z-axis. However, if the implicit equation is more complicated, determining the axis of revolution is not trivial. • Juan G. Alc´azar is with Departamento de F´ısica y Matem´aticas, Universidad de Alcal´a, E-28871 Madrid, Spain. E-mail: [email protected] • Ron Goldman is with Department of Computer Science, Rice University, Houston, Texas 77005, USA. E-mail: [email protected]

The primary contribution of this paper is to develop an algorithm to extract the axis of revolution from the implicit equation of an algebraic surface of revolution. Our algorithm is based on three different computational methods, each one with advantages and disadvantages: factoring the highest order form into quadrics, contracting the tensor of the highest order form, and using resultants and gcds. Two of these methods (factoring the highest order form into quadrics and contracting the tensor of the highest order form) rely on the crucial observation that the highest order form of the implicit equation of an algebraic surface of revolution is also a surface of revolution whose axis passes through the origin and is parallel to the axis of revolution of the original surface (see Section 2.1). Thus the highest order form typically contains enough information to compute the direction of the axis of revolution of an algebraic surface of revolution. A point on the axis of revolution can then be found by straightforward computational techniques. In some cases, however, for instance when the highest order form is equal to (x2 + y 2 + z 2 )p , the surface represented by the highest order form degenerates to a single point, and as a consequence an alternative approach, provided by the third method, is needed. During our analysis, we shall also derive some interesting and novel algebraic properties of algebraic surfaces of revolution. Some of these properties are in fact necessary conditions for an algebraic equation to represent a surface of revolution. These properties concern the factorization of the highest order form of the implicit equation of the surface, as well as the contraction of the tensor associated with this form. The problem of extracting the axis of revolution of a surface of revolution from the equations of the surface has so far not received too much attention. This problem is implicit in [14], although somewhat hidden because the problem addressed there is to approximate noisy data with a surface of revolution. In [14] the input is a point cloud possibly corresponding to a surface of ¨ revolution, and Plucker coordinates, together with least-

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 2

squares techniques, are used to approximate the best ¨ choice for the axis of revolution. Plucker coordinates are also used in a very recent paper [20], developed independently of our paper, where the authors address the very same problem: finding the axis of revolution from the equation of a surface of revolution. The method ¨ in [20] also uses Plucker coordinates, but with an ingenious device that avoids extending the field where the coefficients of the implicit equation are defined, therefore making the method efficient and fast. Essentially, the method in [20] requires computing five independent normals to the surface and leads to a linear system in ¨ the Plucker coordinates of the axis of revolution. From a computational point of view, the algorithm in [20] is certainly sound, although the development of the algorithm does not shed much light on the structure of the implicit equation of a surface of revolution.

In contrast to [20], the development of the algorithm in our paper does provide information on the structure of the implicit equation. Furthermore, our algorithm computes rapidly and efficiently the direction of the axis of revolution, in a wide range of cases. Indeed, when the implicit equation of the surface has odd degree, we prove that the form of highest order of the implicit equation of the surface has a linear factor (see Section 2.1), whose coefficients provide the direction of the axis of revolution. This linear factor can easily be computed either by factoring the polynomial, or by means of an ad-hoc method presented in Section 2.1.1. Moreover, we also prove that the direction of the axis of revolution can be computed from the tensor defined by the form of highest degree of the implicit equation (see Section 2.2). After contraction, this tensor provides either a vector, when the degree of the surface is odd, or a 3 × 3 matrix, when the degree of the surface is even. If the degree is odd and the vector is not the zero vector, then the vector is parallel to the axis of revolution of the surface. If the degree is even and the matrix is not a multiple of the identity, then the axis of revolution is parallel to an eigenvector associated with the only simple eigenvalue of the matrix. When these two approaches (factoring polynomials and contracting tensors) fail, we resort to a third method using resultants and gcds; although this third method always works, it is certainly less efficient than the other two.

Three sections follow. In Section 2 we present our three methods for finding the axis of revolution from the equation of a surface of revolution. In Section 3, we compare and contrast the strengths and weaknesses of these three methods and we propose a general algorithm to solve the problem of extracting the axis of revolution from the equation of a surface of revolution. In Section 4 we compare the performance of our method to the method in [20].

2

M ETHODS FOR FINDING THE AXIS OF REVO LUTION OF A SURFACE OF REVOLUTION . Here we shall consider three methods for finding the axis of revolution from the implicit equation of a surface of revolution: factoring the highest order form into quadrics, contracting the tensor of the highest order form, and using resultants and gcds. In this section, we consider an algebraic surface S, implicitly defined by an equation F (x, y, z) = 0, where F is a squarefree polynomial. The methods we present here are valid whenever the coefficients of F belong to a field of characteristic 0. Nevertheless, the computations greatly simplify when the coefficients are rational numbers. We exclude the case where F depends on just two variables: if, say, F = F (x, y), the analysis reduces to checking whether F (x, y) = 0, seen as an algebraic curve in the xy-plane, is a product of concentric circles. Notice here that in certain cases one can apply a change of variables to reduce the number of variables of a polynomial [2], thereby simplifying the problem. We also exclude the trivial case where F is a plane. In general we will denote space vectors as x = (x, y, z). 2.1 First approach: factoring the highest order form into quadrics. The method we present next focuses on finding the direction of the axis of rotation. Once the direction of this axis is found, we shall see that straightforward algebraic techniques can be used to find a point on this axis. Let F (x) = 0 implicitly define an algebraic surface of revolution S, and let A denote the axis of revolution of S. Since any rotation around A can be represented by a rotation R about a line through the origin parallel to A followed by an appropriate translation b, there is a collection of transformations T (x) = Rx + b such that the surface S is invariant under T (x). We will write F (x) as F (x) = FN (x) + FN −1 (x) + · · · + F1 (x) + F0 (x), where Fi (x), for i = 0, 1, . . . , N , corresponds to the homogeneous form of degree, or order, i. The following result about the highest order form FN (x) is one of the key insights of this paper. Lemma 1. FN (x) also defines a surface of revolution about an axis of revolution A0 through the origin parallel to A. Proof: Let G(x) = F (T (x)) = F (Rx + b). Since F is invariant under T (x), the expressions F (Rx + b) and F (x) both define the same surface; therefore F (x) = λ(R, b) · G(x) for any rotation T about A. In particular, Fk (x) = λ(R, b) · Gk (x) for k = 0, 1, . . . , N . Since GN (x) = FN (Rx) it follows that FN (x) = λ(R, b) · FN (Rx), so the surface defined by FN (x) is invariant under any rotation T˜(x) = Rx. Lemma 2. If the axis A passes through the origin, then for k = 1, 2, . . . , N , Fk (x) defines a surface of revolution about A.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 3

Proof: If the axis A passes through the origin, then for any rotation T (x) = Rx about A, we have Fk (Rx) = λ(R) · Fk (x) for every k = 1, . . . , N . Therefore Fk (x) is invariant under every rotation around A. Corollary 3. F (x) defines a surface of revolution about an axis passing through the origin if and only if for k = 1, 2, . . . , N , Fk (x) defines a surface of revolution about the same axis passing through the origin. Let us consider first the case N = 2, where F (x) defines a quadric surface. In this case the form of highest degree is also a quadric, and this quadric can be written as F2 (x) = xT QF x, where QF is a 3 × 3 symmetric matrix. We shall now show that quadrics of revolution can easily be detected and their axis of revolution can easily be computed. We begin with a technical lemma. Lemma 4. Let F (x) define a quadric of revolution, and let A be the axis of revolution. Then for any rotation matrix R about the line A0 through the origin parallel to A, QF = R−1 · QF · R. Proof: By Lemma 1, F2 (x) defines a surface of revolution about A0 . So F2 (Rx) = λ(R) · F2 (x). Since F2 (Rx) = (Rx)T · QF · (Rx) = xT (RT · QF · R)x, we conclude that λ(R)·QF = RT ·QF ·R. Furthermore, R is a rotation matrix, so R−1 = RT ; therefore QF and λ(R)·QF are similar matrices. Hence, QF and λ(R) · QF have the same eigenvalues. Since the eigenvalues of λ(R) · QF are the eigenvalues of QF multiplied by λ(R), we conclude that λ(R) = 1. The following result provides a method for computing the direction of the axis of revolution of a (possibly degenerate, and possibly complex) quadric of revolution. Theorem 5. Let F (x) define a quadric of revolution different from a sphere, and let A be the axis of revolution. Then the matrix QF has only one simple eigenvalue, and the associated eigenvector is parallel to the axis of revolution. Proof: Let us see first that the axis of revolution corresponds to an eigenvector of QF . Indeed, let R denote any rotation matrix about a line A0 through the origin parallel to the axis of revolution A. If v is parallel to A0 , then R−1 v = v. Thus by Lemma 4 QF v = (R−1 · QF · R)(R−1 v) = R−1 (QF v). So QF v is also invariant under the rotation matrix R; therefore QF v is also parallel to v. Hence, v is an eigenvector of QF . Now let us see that v corresponds to a simple eigenvector λ. Indeed, λ cannot be a triple eigenvalue because then the quadric is a sphere. So QF must have another eigenvalue µ 6= λ. Let w be an eigenvector associated with µ. Hence, w is not parallel to A, and therefore R−1 w is also not parallel to v. Now since QF (R−1 w) = (R−1 · QF · R)(R−1 w) = R−1 (QF w) = R−1 (µw) = µ(R−1 w),

it follows that R−1 w is also an eigenvector of QF corresponding to the eigenvalue µ. However, R−1 w cannot be parallel to w, because otherwise w would be, in turn, parallel to v. So there are at least two independent eigenvectors associated with µ, which is therefore a double eigenvalue of QF . For N > 2, the following theorem provides a connection with quadrics and is the main result of this subsection. Theorem 6. Let F (x) define a surface of revolution with axis A, and let FN (x) be the form of highest degree of F (x). Then: (1) FN (x) factors completely into a single real plane and several quadrics, possibly with multiplicity in both cases.The real plane may or may not be present in the factorization. (2) The quadrics are either circular, possibly imaginary cones with the same axis and vertex, or degenerate quadrics corresponding to a line or a single point. The lines and the plane, if any, also pass through the vertex of the cones. The point, if any, coincides with the vertex of the cones. (3) The direction of the normal vector to the real plane, if any, the directions of the cones’ axes and the direction of the lines corresponding to degenerate quadrics, are the same, and coincide with the direction of the axis of revolution A. Proof: By Lemma 1, FN (x) defines a surface of revolution about a line A0 through the origin, parallel to A. Since the nature of the factors of FN (x) is not affected by a rigid motion, we can begin by assuming that A0 is the z-axis. Since FN (x) is invariant under any rotation around the z-axis, using cylindrical coordinates (ρ, θ,p z) we can write FN = FN (ρ, z). Furthermore, since ρ = x2 + y 2 and FN (x) is algebraic in x, y, z, we have that FN = FN (ρ2 , z). So FN (x) can be written as FN (x) = a0 (x2 + y 2 )p z q + a1 (x2 + y 2 )p−1 z q+2 + · · · +aj (x2 + y 2 )p−j z q+2j , for some j ∈ {0, 1, . . . , p}. Therefore we deduce that FN (x) = (x2 + y 2 )p−j · z q · [aj z 2j + · · · +a1 (x2 + y 2 )j−1 z 2 + a0 (x2 + y 2 )j ].

(1)

The first factor, of multiplicity p − j, corresponds to a degenerate quadric defining a line, namely the z-axis. The second factor, of multiplicity q, is the xy-plane, normal to the z-axis. As for the last polynomial, we shall seek factors for it of the form x2 + y 2 − rz 2 , which correspond to circular cones with vertex at the origin and axis equal to the z-axis. To determine the values of r for which x2 +y 2 −rz 2 is a factor of (1), we plug x2 +y 2 = rz 2 into the last polynomial in (1) and get aj + · · · + a1 rj−1 + a0 rj = 0. This expression is a univariate polynomial in r, so this polynomial has j solutions, possibly with multiplicity. Therefore, we get j factors all corresponding to circular cones with the same vertex, the origin, and the same axis, the z-axis. Notice that if r = 0 is a root, then the cone collapses to a line, the z-axis, and if r = −1 is a root, then

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 4

the cone collapses to a single point, the origin. Since the nature of the factors is not affected by a rigid motion, the theorem follows. Figure 1 illustrates Theorem 6. The following corollaries can be derived directly from Theorem 6.

Direction of A

Fig. 1: Illustrating Theorem 6: the figure shows the different possibilities for the factors of FN (x), and the relative positions of the surfaces defined by these factors with respect to the line through the origin parallel to A. Corollary 7. If F (x) has odd degree, then FN (x) = 0 contains a single real plane, and the direction of the axis of revolution is the direction of the normal vector to this plane. Corollary 8. If the axis of revolution of F (x) goes through the origin and F1 (x) is not identically zero, then F1 (x) divides every form of higher, odd, degree. Furthermore, the normal vector to the plane surface defined by F1 (x), when it is not identically zero, provides the direction of the axis of revolution. Corollary 9. If FN (x) depends on just one variable, then the axis of F (x) is parallel to the corresponding coordinate axis. Furthermore, if FN (x) depends on just two variables, then the axis of F (x) is parallel to the coordinate axis corresponding to the variable not included in FN (x). If F (x) and the quadratic factors of FN (x) have rational coefficients, then the quadratic factors of FN (x) can be found by applying standard factorization techniques, implemented in most computer algebra systems. However, the quadratic factors of FN (x) do not necessarily have rational coefficients, even when the coefficients of F (x) are rational; in fact, the quadratic factors may not even be real. So in the most general case we may need to factor over the complex numbers, i.e. to compute an absolute factorization of FN (x). This factoring problem has been addressed for instance in [6], [9], [13], but we must say that factoring over the complex numbers is not an easy task. Nevertheless once the factors are known, the direction of the axis of revolution can be found by applying Theorem 5 to any of the quadratic factors. In order to actually find a point on the axis A, one can pick any plane π normal to the direction of A. The

intersection of the surface of revolution S and the plane π must be a product of concentric circles with a common center corresponding to a point on A. By applying a rotation that transforms A into the new z-axis, S ∩ π is transformed onto a planar algebraic curve g(x, y) = 0, which is again a product of concentric circles. So their center can be found, for instance, by computing the leftmost and right-most points of this curve or the highest and lowest points of this curve, and then determining the middle point. If g˜(x, y) denotes the square-free part of g(x, y), the x-coordinates of the left-most and right-most points to the largest and smallest real roots of   correspond ∂g ˜ . Therefore the x-coordinate of the center is Resy g˜, ∂y the arithmetic mean of these two numbers. Similarly, the y-coordinate of the center is the arithmetic mean of the   ∂g ˜ largest and smallest real roots of Resx g˜, ∂x . Another method for finding a point on the axis A will be given in Section 2.3. Example 1. Consider the surface F (x) = (x2 +y 2 )2 +z 4 −z 3 . The form of highest degree is FN (x) = (x2 + y 2 )2 + z 4 , which factors into two imaginary circular cones around the z-axis:   FN (x) = (x2 + y 2 + iz 2 ) · (x2 + y 2 − iz 2 ) . The matrices associated with these two quadratic factors have 1 as a double eigenvalue, and a simple eigenvalue equal to i for the first factor, and −i for the second factor. The eigenvector associated with these simple eigenvalues is parallel to (0, 0, 1). So the axis of revolution is parallel to the z-axis. Now intersecting the surface F (x) = 0 with the plane z = 0 yields the curve {x2 +y 2 = 0, z = 0}, which contains just one real point, the origin. Therefore, the axis of revolution is the z-axis. In addition, FN −1 (x) = −z 3 also defines a surface of revolution about the z-axis, namely the plane z = 0, as predicted by Corollary 3. Notice that if FN (x) = (x2 +y 2 +z 2 )p , which for example happens with cyclides (see Exercise 19.3, page 305 in [11]), Theorem 6 does not provide any information about the direction of A, since in this case FN (x) collapses to a single point. We refer the reader to Section 2.3 for an alternative method for analyzing this case. 2.1.1 The case N odd. If N is odd, then by Corollary 7 the polynomial FN (x) necessarily has one real linear factor. Furthermore, this linear factor provides the direction of the axis of revolution A. By Corollary 2.6 in [20], the axis of revolution is parallel to a vector with rational coefficients. Therefore, one can find the direction of the axis of revolution by just factoring over the rationals the highest order form of the implicit equation, and picking the linear factor. Notice that factoring over the rationals is a very fast and standard technique implemented in the main computer algebra systems. If the coefficients of the implicit equation are not rational and N is odd, we can apply a computational approach to find the linear factor of FN (x)

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 5

without completely factoring FN (x). We will assume that FN (x) is square-free, since otherwise we can easily remove multiple factors by dividing FN (x) by  ∂FN ∂FN N , , . Furthermore, we will focus on gcd FN , ∂F ∂x ∂y ∂z factors that are not of the form Ax + By, since these factors lead to a special case that can be treated similarly. Now FN (x) has a linear factor iff the surface S 0 defined by FN (x) contains a plane π that we will assume is not normal to the xy-plane. Since FN (x) is by assumption square-free, the singular locus of S 0 is at most 1dimensional. Therefore, for almost all lines parallel to the z-axis the intersection with S 0 does not contain any singular point of S 0 . So let {x = x0 , y = y0 } contain only regular points of S 0 ∩ π. This line must surely intersect π. Let Q = (x0 , y0 , β) represent the intersection of the line {x = x0 , y = y0 } with the plane π; notice that β satisfies ξ(β) = FN (x0 , y0 , β) = 0. Since by hypothesis this line does not intersect the singular locus of S 0 , the gradient ∇ (FN ) (x0 , y0 , β) = (a(β), b(β), c(β)) 6= 0. But the gradient is normal to the surface S 0 , so the gradient is normal to the plane π. Therefore π can be written as π : a(β)(x − x0 ) + b(β)(y − y0 ) + c(β)(z − β) = 0. Since π is a component of S 0 , the resultant of π and FN with respect to any variable x, y, z is identically 0. Writing the resultant with respect to, say x, of these P polynomials as Resx (π, FN ) = i,j ηij (β)y i z j , it follows that β is a real root of the greatest common divisor of the ηij (β)’s and ξ(β). Notice that in our case we have just one real linear factor, so β is in fact the real root of the greatest common divisor of the ηij (β)’s and ξ(β). Factors of the type Ax + By can be found in a similar way. √ √ 3 2 y − 8 2z 3 x2 − Example 2. Let G(x) = −16z√4 x − 2 2z √ √ 10z 2 xy 2 +8z 2 x3 +5 2zx2 y 2 − 2zy 4 +4 2x4 z+x3 y 2 +xy 4 , which is a homogeneous polynomial of degree 5. To find the linear factors of G(x), let us pick the line {x0 = −1, y0 = 1}. The intersections of this line with the surface defined by G(x, y, z) are the points (−1, 1, β), with √ √ √ ξ(β) = ( 2β + 1)(2β 2 − 2 2β + 2)(−4 2β + 1) = 0. Furthermore, ∇(G)(−1, 1, β) = (A(β), B(β), C(β)), where: √ √ 2 A(β) = 16 2β 3 − 16β 4 + 14β − 26 2β + 4, √ √ 2 3 B(β) = 20β √ −2 4 2β 3+ 6 2β −√6, C(β) − 30 2β + 64β + 4β + 8 2. Let π(x, y, z) = A(β)(x+1)+B(β)(y−1)+C(β)(z−β) = 0. After computing Resx (G(x, y, z), π(x, y, z)), we find that the gcd of the coefficients of Resx (G(x, y, z), π(x, y, z)) and ξ(β) is √ √ 2 2 2 3 β − β + . 2 2 √









The roots of this polynomial are − 22 , 22 − 22 i, 22 + 22 i. So we get three linear factors of G(x) corresponding to these three β-values, one of them real, the other two complex conjugate. In fact, the two complex conjugate factors give rise√ to a real quadratic factor of G(x). Substituting β = − 22 in

√ (A(β), B(β), √ C(β)), we√get (25, 0, −25 2), which is parallel to (1, 0, − 2), so x − 2z is a linear factor of G(x). 2.1.2 Summary of the algorithm corresponding to the first approach. We finish the subsection with an algorithm summarizing this computational method. Notice that this method is useful either when the degree N is odd, or when N is even but one of the quadratic factors of FN (x) has rational coefficients. The main computational tools that we need are factoring, or, if N is odd, resultants and gcds, which are required to apply the method in Section 2.1.1. Subalgorithm 1 Input: A square free polynomial F (x) depending explicitly on x, y, z, defining a surface of revolution S. Output: The axis of revolution A of the surface S. If FN (x) is equal to (x2 + y 2 + z 2 )p , return Use a different method. Otherwise consider two cases: Case 1: N is odd. If F (x) has rational coefficients, factor FN (x) and pick the linear factor. This factor provides the direction of A. Otherwise, apply the method in Section 2.1.1 to find the direction of A. Compute a point on A. Case 2: N is even: If FN (x) has some quadratic factor with rational coefficients, then find the direction of A by using Theorem 5. Then compute a point on A. Otherwise, Use a different method.

2.2 Second approach: contracting the tensor of the highest order form. In the case when N is even, factoring the highest order form into quadrics may require factoring over the complex numbers -that is, we may require an absolute factorization. We now present an alternative method that also focuses on the highest order form, which avoids factoring altogether. The method of contracting the tensor of the highest order form, like the method of factoring the highest order form into quadrics, concentrates on finding the direction of the axis of revolution. Once this direction is found, we can use the straightforward method of the previous section to find a point on the axis of revolution. Informally, a tensor is an entity depending on several indices, where the nature of each index corresponds to one of two possible categories (contravariant or covariant). A slightly more formal description is that a tensor is a multidimensional array defined relative to a vector space with a fixed basis whose components transform in a certain way when the basis is changed. In our case the ambient vector space is a Euclidean space of dimension

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 6

M . Readers unfamiliar with tensors may find useful the book [19], as well as Chapter 10 in [7]. To get better acquainted with tensors, let us examine first two familiar examples: vectors and matrices. A contravariant vector is an entity  T = ti where i ∈ {1, . . . , M }, such that when a change of basis of the underlying vector space represented by a square matrix Q = (Qri ) is applied, the components ti of T transform by the rule X tˆi = Qir · tr . r=1,...,M

For instance, the coordinates of a vector with respect to a basis provide an example of a contravariant vector; in fact, in the above expression one can recognize the usual formula for the new coordinates of a vector after a change  of basis. Similarly, a covariant vector is an entity T = tj , with j ∈ {1, . . . , M }, whose components tj transform by the rule X tˆj = (Q−1 )sj · ts . s=1,...,M

n Q = (Qri ) is applied, the components tij11,...,i ,...,jm transform by the rule in  −1 i1 n m n ·trs11 ,...,r · · · Q−1 rn ·Qsj11 · · · Qsjm tˆji11,...,i ,...,sm . (4) ,...,jm = Q r1

This expression is certainly cumbersome, but fortunately we do not really need to make use of it, since eventually all of our calculations will lead to either vectors or matrices. An n-contravariant and m-covariant tensor is said to be of order (n, m). If n + m = 1, then the tensor is in fact a vector, contravariant if n = 1, m = 0 and covariant if n = 0, m = 1. If n + m = 2, the tensor can be represented by a matrix. We say that a tensor is symmetric if the components of the tensor are invariant under any permutation of the indices; for instance, symmetric matrices correspond to symmetric tensors of order two. There are several operations that can be carried out with tensors. For us the most useful operation is contraction. Contraction reduces a tensor of order (n, m) to another tensor of order (n − 1, m − 1) by summing over one contravariant and one covariant index. Summing over the last contravariant and covariant indices of a tensor T , the contraction of T produces a new tensor S with components i ,...,i

An example of a covariant vector is the list of the covariant coordinates of a vector in a certain basis, which are the projections of the vector onto the directions of the basis vectors (see Section 10.2 in [7]). Both in the case of contravariant and covariant vectors we have just one index; therefore, these vectors correspond to tensors of order 1. Tensors of order 2 can be represented by matrices. For instance, a 1-contravariant and 1-covariant tensor is  an entity T = tij , where i, j ∈ {1, . . . , M }, whose components tij transform in the following way under a change of basis with change of basis matrix Q: X X tˆij = (Q−1 )sj · Qir · trs . (2) r=1,...M s=1,...,M

In this formula one can recognize the matrix equality ˆ = Q−1 · T · Q, where T is the matrix (ti ), and T ˆ is the T j i matrix (tˆj ). Therefore, an example of a 1-contravariant and 1-covariant tensor is the set of entries of a matrix associated with a linear mapping. In order to simplify the notation of (2), the Einstein summation notation is commonly used. Here, whenever a subindex and superindex agree, it is understood that we must sum over the range of variation of that index. Using this notation, we can rewrite equation (2) as tˆij = (Q−1 )sj · Qir · trs .

(3)

In general, a mixed n-contravariant and m-covariant tensor is a multidimensional array  n T = tij11,...,i ,...,jm where ip , jq ∈ {1, . . . , M } for p = 1, . . . , n, q = 1, . . . , m, such that when a change of basis defined by a matrix

i ,...,i

,1

i ,...,i

,M

= tj11 ,...,jn−1 sj11 ,...,jn−1 + · · · + tj11 ,...,jn−1 . m−1 m−1 ,1 m−1 ,M  i For instance, a tensor T = tj of order (1, 1) can be contracted to a tensor of order 0, the scalar equal to t11 + t22 + t33 , i.e. the trace of the matrix representing the tensor. Observe that when n − m = 0, by repeteadly contracting we eventually arrive at a tensor of order (1,1), which corresponds to a matrix. If n − m = 1 or m − n = 1, by repeatedly contracting we eventually arrive at a (contravariant or covariant) vector. Since we work in a Euclidean space with a Euclidean metric, any index can be considered as either covariant or contravariant. Therefore, one does not need to worry about the contravariant or covariant nature of every index (see page 240 in [7]). The reason why tensors are useful for us is that a homogeneous polynomial f in the variables x, y, z can be represented as a symmetric (n, m) tensor, where n+m is equal to the degree N of the polynomial [5]. We will consider n − m = 0 if N is even, and n − m = 1 if N is odd. More precisely, let x := x1 , y := x2 , z := x3 , and let (x1 , x2 , x3 ) be the contravariant coordinates of (x1 , x2 , x3 ); recall that since the ambient space is R3 , in fact xk = xk for k = 1, 2, 3. Using Einstein’s summation notation j1 jm n f (x1 , x2 , x3 ) = cij11,...,i ,...,jm xi1 · · · xin x · · · x ,

where ip , jq ∈ {1, 2, 3} for any p, q, i1 + · · · + in + j1 + n · · · + jm = N , and cji11,...,i ,...,jm is equal to the coefficient in f of the monomial xi1 · · · xin xj1 · · · xjm , divided by the number of permutations of {i1 , . . . , in , j1 , . . . , jm }.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 7

Example 3. Consider the homogeneous polynomial of order 3 defined by f (x, y, z) = 8x3 + 12x2 y − 24x2 z + 6xy 2 − 24xyz + 24xz 2 +y 3 − 6y 2 z + 12yz 2 − 8z 3 . In this case, we have: 3 c1,1 8, 1 = {coeff. of x } = 2,2 3 1, c2 = {coeff. of y } = 3 c3,3 3 = {coeff. of z } = −8. 2,1 1,2 The components c1,1 2 = c1 = c1 are equal to the coefficient 1,1 2,1 2 of x y in f , divided by 3, i.e. c2 = c1,2 1 = c1 = 12/3 = 4. Similar results hold for other components with two repeated indices. Also, we have that 3,1 2,3 3,1 1,3 2,1 c1,2 3 = c3 = c2 = c2 = c1 = c2 ,

and each of these components is equal to the coefficient of xyz in f , divided by 6, i.e. −24/6 = −4. Similarly, we have

ˆ = λ(R) · C. Therefore this result follows by know that C the same arguments as in Lemma 4 and Theorem 5. ˆ = R · c. On the (2) Since C i1 is a contravariant vector c ˆ = λ(R)·c, so the direction of the vector c 6= other hand, c 0 is unchanged by the rotation defined by R. Therefore, c is parallel to A0 . Example 4. Consider the surface defined by F (x) = 8x3 + 12x2 y − 24x2 z + 6xy 2 − 24xyz + 24xz 2 +y 3 − 6y 2 z + 12yz 2 − 8z 3 − 51x2 − 24xy + 48xz − 33y 2 +24yz − 51z 2 + 36x − 9y − 90z − 54. Notice that the form of F of highest degree is the polynomial in Example 3. Therefore 1,3 1,2 8+2+8= 18 c1 = c1,1 1 + c2 + c3 = 2,1 2,2 2 c = c1 + c2 + c2,3 = 4 + 1 + 4 = 9 3 3,3 3,2 = (−8) + (−2) + (−8) = −18 + c + c c3 = c3,1 3 2 1

Hence c = (18, 9, −18) is parallel to the axis of revolution. In fact, the form of highest degree of F (x) factors as (2x + 1,3 3,1 c1,1 = −8 1,2 2,1 2,2 y − 2z)3 , which is consistent with Theorem 6. 3 = c1 = c1 c2 = c2 = c1 = 2 3,1 3,3 c1,3 = c = c = 8 2,3 3,2 3 3 1 c2,2 = −2 Notice that when N is odd, Theorem 10 provides 3,2 3,3 3 = c2 = c2 c2,3 = 4 3 = c3 = c2 the direction of A whenever c 6= 0. When N is even, k Theorem 10 also provides the direction of A whenever C and f (x1 , x2 , x3 ) = ci,j k · xi xj x . is not a multiple of the identity matrix. Unfortunately, if Now let F (x) define a surface of revolution S of FN (x) = (x2 +y 2 +z 2 )p , where the methods of Section 2.1 degree N about an axis A. From Lemma 1, FN (x) also are not applicable, C is a multiple of the identity matrix. defines a surface of revolution S 0 about an axis A0 Indeed, in this situation unless i = j, the component through the origin parallel to A. Therefore any rotation ci,k...,k j,k,...,k is 0, because this component corresponds to a T (x) = Rx about A0 leaves S 0 invariant, i.e. FN (Rx) = term of FN (x) where the power of some variable is odd. λ(R) · FN (x). Since FN (x) is a homogeneous polynomial Moreover, by symmetry c11 = c22 = c33 . On the other ,...,in , in x, y, z, FN (x) can be represented by a tensor Cji11,...,j m hand, if C is a multiple of the identity matrix, that does where n − m = 0 if N is even, and n − m = 1 if N is odd. not necessarily imply that FN (x) = (x2 + y 2 + z 2 )p . For Similarly FN (Rx) can also be represented by another 3 3 2 3 ,...,in ˆ i1 ,...,in = λ(R) · C i1 ,...,in . By instance, if FN (x) = xy − x y + x yz − y z, then C is tensor, Cˆji11,...,j , where C j1 ,...,jm j1 ,...,jm m the null matrix. repeatedly contracting, we eventually reach either the Once the direction of A is found, a point on A can be ˆ i1 = λ(R) · C i1 , if N is even, or the equality equality C j1 determined as explained in Section 2.1. Another method j1 ˆ i1 = λ(R) · C i1 , if N is odd. Here we represent in for determining a point on A, as well as an alternative C bold face the new tensors generated by contracting the method for the cases when Theorem 10 fails to provide the direction of A, is presented in the next section. original tensors.  Let C denote the matrix representing the tensor C ij11 in the case N even, and let c denote the  vector C i1 in the case N odd. The following theorem is 2.2.1 Summary of the algorithm corresponding to the our main result connecting tensors of the highest order second approach. We finish the subsection with an algorithm summarizing form to the direction of the axis of revolution. this computational method. This method is useful if by ,...,in Theorem 10. Let Cji11,...,j be the tensor representing the m contraction we reach either a nonzero vector c, or a surface of revolution FN (x). Then: matrix C that is not a multiple of the identity matrix. (1) If N is even and the matrix C is not a multiple of the Notice that this method relies only on Linear Algebra, so identity matrix, then the matrix C has only one simple we do not need to use any sophisticated computational eigenvalue, and the associated eigenvector is parallel to tools. the axis of revolution A. (2) If N is odd and c 6= 0, the vector c is parallel to the axis 2.3 Third approach: resultants and gcds. of revolution A. Factoring the highest degree form into quadrics and Proof: (1) Since C ij11 is a (1, 1) tensor, the change of contracting the tensor of the highest degree form of F (x) basis defined by R, where R is any rotation about A0 , revolution when the highest ˆ representing both fail to find the axis of transforms this tensor so that the matrix C 2  i  degree form F (x) = (x + y 2 + z 2 )p . Contracting the N 1 −1 ˆ ˆ C j1 satisfies C = R · C · R. On the other hand, we tensor of the highest degree form also fails when we 1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 8

Subalgorithm 2 Input: A square free polynomial F (x) depending explicitly on x, y, z, defining a surface of revolution S. Output: The axis of revolution A of the surface S. If FN (x) is equal to (x2 + y 2 + z 2 )p , return Use a different method. Otherwise consider two cases. Case 1: N is odd: Contract the tensor associated with FN (x) to reach a vector c. If c = 0 return Use a different method. If c 6= 0, A is parallel to c. Then find a point on A. Case 2: N is even: Contract the tensor associated with FN (x) to reach a matrix C. If C is a multiple of the identity matrix, return Use a different method. If C is not a multiple of the identity, find the eigenvector associated with the only simple eigenvalue; A is parallel to this eigenvector. Then find a point on A.

the axis of revolution A0 of the surface S 0 defined by G(x, λ0 ) passes through the origin.

contract to a matrix that is a multiple of the identity matrix, or to the zero vector. Therefore we seek an alternative way that will work when the other methods fail. Our final approach uses only univariate resultants and gcds. By assumption, F is square-free and therefore the singular locus of the surface S defined by the polynomial F is at most 1-dimensional. Hence, we can always find (x0 , y0 ) such that the line defined by the equations {x = x0 , y = y0 } does not contain any singular point of S. The intersection between this line and the surface S consists of the points (x0 , y0 , z) where F (x0 , y0 , z) = 0. Let P = (x0 , y0 , z0 ) be one of these points. Then P is a regular point of S, and therefore ∇(F )(P ) = (a, b, c) 6= 0 is parallel to the normal line to S at P , which can be written as

where

Lemma 11. For a generic x0 = (x0 , y0 , z0 ) ∈ S, both GN −1 (x, λ) and G1 (x, λ) are not identically zero. Proof: Let us see first that G1 (x, λ) is not identically zero. Indeed, if G1 (x, λ) is identically zero, then in particular G1 (x, 0) is also identically zero. Since G(x, 0) = F (x + x0 ), by the chain rule ∇(G(x, 0))x=0 = ∇(F (x0 )). Therefore since 0 = ∇(G1 (x, 0)) = ∇(G(x, 0))x=0 , it follows that ∇(F (x0 )) = 0, which is impossible because generically x0 is a regular point of S. Hence, G1 (x, λ) is not identically zero. Now let us see that GN −1 (x, λ) is also not identically zero. Let X FN (x) = aijk · xi y j z k . (5) 0 ≤ i, j, k ≤ N i+j+k =N

Then GN −1 (x, λ) =

X 0 ≤ i, j, k ≤ N i+j+k =N

aijk ·ψijk (x, x0 , λ)+FN −1 (x), (6)

x = x0 + λ · a, y = y0 + λ · b, z = z0 + λ · c. Let x0 = (x0 , y0 , z0 ) and a = (a, b, c); then Pλ = x0 + λa represents a generic point on the normal line to the surface S at the point P . Since every normal line to S intersects the axis of revolution A, there exists some λ0 such that the point x0 + λ0 a ∈ A. Consider now the translation x → t(x, λ) = x + (x0 + λa). This translation transforms S into a new surface S 0 which is also a surface of revolution, with an axis of revolution A0 parallel to A. We will denote the polynomial defining S 0 by G(x, λ) = F (t(x, λ)), and will denote also the forms of degree N − 1 and degree 1 of G(x, λ) by GN −1 (x, λ) and G1 (x, λ). Notice that when λ = λ0 ,

ψijk (x, x0 , λ) = i(x0 + λa)xi−1 y j z k + j(y0 + λb)xi y j−1 z k +k(z0 + λc)xi y j z k−1 . If FN depends on just one coordinate, say z, then (6) is simpler and equal to GN −1 (x, λ) = N a0,0,N (z0 + λc)z N −1 + FN −1 (x). Generically c 6= 0, since otherwise F = F (x, y), which is a case excluded in this paper. Therefore, in this particular case GN −1 (x, λ) is not identically zero. If FN depends on two coordinates, say x, y, then the expression for GN −1 (x, λ) is also simpler, and reduces to X GN −1 (x, λ) = ai,j,0 ·ψi,j (x, x0 , λ)+FN −1 (x), 0 ≤ i, j ≤ N i+j =N

(7) where ψi,j (x, x0 , λ) = i(x0 + λa)xi−1 y j + j(y0 + λb)xi y j−1 . Generically a, b cannot both be identically zero, since otherwise the normal to S would be parallel to the zaxis, implying that S is a plane, contrary to assumption. Hence, once again GN −1 (x, λ) is not identically zero. Now let us finally consider the general case, where FN depends on x, y, z. In this situation, since (x0 , y0 , z0 ) is generic (x0 , y0 , z0 ) is regular; therefore a, b, c cannot all be zero. So from the right hand-side of (6) we observe that GN −1 (x, λ) is also not identically zero. From Lemma 11, we get an algorithm to find a point on A, given the direction ~v of A. The idea is the following: after picking a generic point (x0 , y0 , z0 ) ∈ S, we

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 9

have G1 (x, λ0 ) = A(λ0 )x + B(λ0 )y + C(λ0 )z. When λ = λ0 , by Lemma 2 the vector (A(λ0 ), B(λ0 ), C(λ0 )) is either 0 or parallel to the direction of the axis A. In the first case gcd(A(λ), B(λ), C(λ)) 6= 1, and the roots of the gcd provide tentative values for λ0 . In the second case, gcd(A(λ), B(λ), C(λ)) = 1, and imposing (A(λ), B(λ), C(λ)) = α · ~v , with α 6= 0, one can solve for λ0 . Notice that each potential value of λ0 provides a different tentative point on the axis of revolution, namely the point x0 + λ0 a. Therefore, in general we get several candidate lines Ai . One can test each of these lines by checking if S is invariant under a generic rotation around each line Ai . We need to see now how to proceed if the direction of A is not previously known. But before addressing this question, we need the following result. Lemma 12. An algebraic surface F (x), where F is a real polynomial, is a product of (possibly complex) spheres centered at the origin if and only if every nonzero homogeneous form of F (x) can be written as αi · (x2 + y 2 + z 2 )i , where αi ∈ R.

Proof: (⇒) Let F (x) = (x2 + y 2 + z 2 − r1 )(x2 + y 2 + z − r2 ) · · · (x2 + y 2 + z 2 − rn ), where r1 , . . . , rn ∈ C. Then F (x) = g(x2 + y 2 + z 2 ), where g(ω) is the univariate polynomial g(ω) = (ω − r1 ) · · · (ω − rn ). Writing 2

g(ω) = an ω n + an−1 ω n−1 + · · · + a1 ω + a0 , Pn where ai ∈ R for i = 0, . . . , n,Pwe have F (x) = i=0 ai · n 2 2 2 i (x2 +y 2 +z 2 )i . (⇐) If F (x) = i=0 αi ·(x ) , then Pn +y +z 2 2 2 i F (x) = g(x + y + z ) where g(ω) = i=0 αi ω . Writing g(ω) = αn (ω − r1 ) · · · (ω − rn ), where r1 , . . . , rn are the roots of g(ω), we conclude that F (x) is the product of spheres centered at the origin. Let us see now how to use Lemma 11 to find both the direction and a point of A, when the direction of A is not previously known. Recall that G1 (x, λ) = A(λ)x + B(λ)y + C(λ)z. We will distinguish here the cases where either gcd(A(λ), B(λ), C(λ)) = 1 or gcd(A(λ), B(λ), C(λ)) 6= 1. (I) Let gcd(A(λ), B(λ), C(λ)) = 1. By Lemma 11, generically GN −1 (x, λ) and G1 (x, λ) are not identically zero. By Corollary 8 when λ = λ0 , G1 (x, λ0 ) divides GN (x, λ0 ), if N is odd, or GN −1 (x, λ0 ), if N is even. Assume now that the common factor of G1 (x, λ0 ) and Gα (x, λ0 ), where α = N or α = N −1 according to whether N is odd or even, depends on the variable z; we can argue in a similar way if instead the common factor depends on the variable x or the variable y. Then the resultant ξ(x, y, λ) = Resz (G1 (x, λ), Gα (x, λ)) must vanish for λ = λ0 . Now we distinguish two cases:

(I.1) (Generic case) If ξ(x, y, λ) is not identically zero, then λ0 must be a root of the gcd of the coefficients ξij (λ) of ξ(x, y, λ), seen as a polynomial in x, y. Furthermore, if λ0 corresponds to a point on A, then (A(λ0 ), B(λ0 ), C(λ0 )) is parallel to A. Moreover, the point on A corresponding to λ0 is x0 + λ0 a. Therefore, for each tentative value of λ0 we get a candidate line Ai . (I.2) (Special case) If ξ(x, y, λ) is identically zero, then G1 (x, λ) divides Gα (x, λ) as a polynomial in x and λ, and therefore we cannot proceed as in case (I.1). However, we know that the direction of the axis A is parallel to (A(λ), B(λ), C(λ)), for some unknown value λ0 of λ. Now we can find tentative values for λ by imposing the condition that S is invariant under a random rotation about the line Lλ parallel to (A(λ), B(λ), C(λ)) through the point Pλ . In detail, let Tλ,θ denote the rotation by a random angle θ about the line Lλ . Tλ,θ is an affine transformation whose equation can be computed by well-known methods (see [11]). Since Tλ,θ leaves the surface S invariant for λ = λ0 , we have that F (Tλ0 ,θ (x)) and F (x) define the same surface; therefore, the coefficients of F (Tλ0 ,θ (x)) and F (x) are proportional. In particular, the polynomials F (Tλ0 ,θ (x)) and F (x) share a factor, and therefore the resultant of F (Tλ,θ (x)) and F (x) with respect to any variable, x, y or z, must vanish identically. Solving for the gcd of the coefficients of these resultants provides tentative values for λ0 , and therefore candidate lines Ai . (II) Let gcd(A(λ), B(λ), C(λ)) 6= 1. By proceeding as in case (I), we can get some candidate lines Ai . However, now it can happen that λ0 is a common root of A(λ), B(λ), C(λ). In this case (A(λ0 ), B(λ0 ), C(λ0 )) = 0, and even though we still have a point on the axis of revolution, namely the point x0 + λ0 a, the direction of the axis of revolution is not known. So let us see how to proceed with the values λ1 , . . . , λn that are roots of gcd(A(λ), B(λ), C(λ)). For each λj , consider G(x, λj ) = F (t(x, λj )). (II.1) If G(x, λj ) has some nonzero form Gi (x, λj ) of odd degree, then we can find the direction of A using the method in Section 2.1.1. (II.2) Assume that all the forms of odd degree of G(x, λj ) are identically zero. If G(x, λj ) defines a surface of revolution, then by construction the axis of revolution A0 , parallel to A, passes through the origin. Then by Lemma 2, every form Gi (x, λj ) defines a surface of revolution about A0 . If every nonzero form can be written as Gi (x, λj ) = αi · (x2 + y 2 + z 2 )i with αi ∈ R, then by Lemma 12 G(x, λj ) is a product of (possibly complex) spheres centered at the origin.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 10

Hence, S is also a product of (possibly complex) spheres, with a common center. So assume that we are not in this case. Since we already know that the axis of revolution A0 passes through the origin, we just need another point of A0 , different from the origin, to compute the direction of A0 , which is also the direction of A. Since G(x, λj ) is not a product of spheres centered at the origin, we can find a new point on A0 by applying again the technique of this subsection on the surface defined by G(x, λj ). In turn, we get tentative values for the direction of A, and therefore we get several candidate lines Ai . Notice that in both cases (II.1) and (II.2), in general we have several candidate lines Ai for the axis of revolution A of the surface S. Hence, in order to test each Ai , one checks whether or not S is invariant under a generic rotation about Ai . Furthermore, we conjecture that case (I) is generic. Observe that while the approaches in Section 2.1 and Section 2.2 are not applicable when FN (x) = (x2 + y 2 + z 2 )p , the approach described in this section is certainly applicable in this case. Notice also that although x0 , y0 can be chosen to be rational, z0 need not be rational. In this case we work symbolically with z0 , representing z0 as a number β that is a root of the polynomial F (x0 , y0 , β) = 0. Simplifications of expressions containing such a β, and gcd’s of polymomials whose coefficients are polynomials in β can be efficiently computed. For instance, one can find algorithms for performing these computations in the computer algebra system Maple. Example 5. Consider the surface implicitly defined by F (x) = (x2 + y 2 + z 2 )4 + 4x3 − 9x2 z + 4xy 2 +4xz 2 − 9y 2 z + 6x2 − 18xz + 2y 2 + 2z 2 + 4x − 9z + 2. This surface is constructed in the following way: first we consider the surface of revolution generated by rotating the curve {g(y, z) = (y 2 + z 2 )4 − 9y 2 z + 1, x = 0} around the z-axis, and then we apply the translation t˜(x) = (x, y, z) → (x + 1, y, z). In particular, notice that the form of highest degree of F (x) is of the type (x2 + y 2 + z 2 )p . By following the ideas in this section, we first consider the point x0 = 1, y0 = −1 on the xy-plane. The line {x = x0 , y = y0 } intersects the surface S at the point (1, −1, β), where β 4 + 10β 2 − 45β + 26 = 0. Furthermore, ∇(F )(1, −1, β) = (a, b, c), where 2

2

a = 8β − 36β + 40, b = −4β + 18β − 20, c = 4β(β 2 + 5) − 45. By substituting x → t(x, λ) = x + (1 + λ · a, −1 + λ · b, β + λ · c),

we get a polynomial G(x, λ), where the form of degree 1 is G1 (x, λ) = A(λ)x + B(λ)y + C(λ)z, with A(λ) = −233280β 3 λ3 − 6912β 3 λ2 − 1054944β 2 λ3 −432β 3 λ − 10896β 2 λ2 + 4265712βλ3 + 840β 2 λ +60696βλ2 − 1793120λ3 + 8β 2 − 864βλ + 6840λ2 −36β + 2388λ + 40, B(λ) = 116640β 3 λ3 + 3456β 3 λ2 + 527472β 2 λ3 +216β 3 λ + 5448β 2 λ2 − 2132856βλ3 − 420β 2 λ −30348βλ2 + 896560λ3 − 4β 2 + 432βλ − 3420λ2 +18β − 1194λ − 20, C(λ) = −244944β 3 λ3 + 13344β 3 λ2 + 663408β 2 λ3 +96β 3 λ + 108β 2 λ2 + 49520βλ3 + 4β 3 + 108β 2 λ −49716βλ2 − 420948λ3 + 2004βλ + 12096λ2 +20β − 2700λ − 45. One can check (for instance using the computer algebra system Maple) that gcd(A(λ), B(λ), C(λ)) = 1. On the other hand, the form of degree 3 of G(x, λ) is G3 (x, λ) = (32β 2 λ − 144βλ + 160λ + 8)x3 +(16β 3 λ + 80βλ + 4β − 180λ − 9)x2 z +(−16β 2 λ + 72βλ − 80λ − 4)x2 y +(32β 2 λ − 144βλ + 160λ + 8)xz 2 +(32β 2 λ − 144βλ + 160λ + 8)xy 2 +(16β 3 λ + 80βλ + 4β − 180λ)z 3 +(−16β 2 λ + 72βλ − 80λ − 4)yz 2 +(16β 3 λ + 80βλ + 4β − 180λ − 9)y 2 z +(−16β 2 λ + 72βλ − 80λ − 4)y 3 . Resz (G1 (x, λ), G3 (x, λ)) is a polynomial in x, y with coefficients that are polynomials in λ. The gcd of these coefficients is 397 2 2340 1563 369 3 β − β − β+λ+ . − 3598 1799 1799 1799 So by setting this gcd equal to zero we get 369 3 397 2 2340 1563 β + β + β− . 3598 1799 1799 1799 Substituting λ0 into G1 (x, λ) yields G1 (x, λ0 ) = C(β)z, where C(β) = 429722219374065 β 3 + 945950869303227 β2 831755057 831755057 λ0 =

+ 11040214482261369 β− 1663510114

3318489360041355 . 831755057

Therefore the direction of the axis of revolution is parallel to (0, 0, C(β)), or equivalently parallel to (0, 0, 1). Finally, since 1 + λ0 a = −1 −1 + λ0 b = 0, 3 β + λ0 c = 14499 3598 β +

16389 2 1799 β

+

95499 1799 β



51147 1799 ,

the point   14499 3 16389 2 95499 51147 −1, 0, β + β + β− 3598 1799 1799 1799 belongs to the axis of revolution A. Therefore, A is the line {x = −1, y = 0}.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 11

2.3.1 Summary of the algorithm corresponding to the third approach. We finish this subsection with an algorithm summarizing this computational method. The main computational tools that we need here are resultants and gcds. Subalgorithm 3 Input: A square free polynomial F (x) depending explicitly on x, y, z, defining a surface of revolution S. Output: The axis of revolution A of the surface S. 1: Let x0 = (x0 , y0 , β), with F (x0 , y0 , β) = 0, such that the line {x = x0 , y = y0 } does not contain any singularity of F (x) = 0. 2: Compute a = ∇F (x0 ). 3: Apply the translation x → t(a, λ) = x + x0 + λa to S. Let S 0 be the new surface, and let G(x) be its implicit equation. 4: Let G1 (x) = A(λ)x + B(λ)y + C(λ), GN −1 (x) be the homogeneous forms of G(x) of degrees 1, N − 1. If either one is identically zero, restart the algorithm. 5: Generic case: gcd(A(λ), B(λ), C(λ)) = 1, and the resultant Resz (G1 , GN −1 ) is not identically zero. 6: Compute the gcd h(λ) of the coefficients of the resultant Resz (G1 , GN −1 ). Each root λi of h(λ) provides a tentative point Pi = x0 + λi a ∈ Ai , and a direction vi = (A(λi ), B(λi ), C(λi )) parallel to Ai . 7: Test each tentative axis Ai defined by Pi and vi . 8: Special Case 1: gcd(A(λ), B(λ), C(λ)) = 1 but the resultant Resz (G1 , GN −1 ) is identically zero. 9: Proceed as in (I.2) 10: Special Case 2: gcd(A(λ), B(λ), C(λ)) 6= 1. 11: Proceed as in (II).

3

T HE

GENERAL ALGORITHM .

We have presented three different computational methods for finding the axis of revolution of a surface of revolution: (1) Factoring the highest order form into quadrics; (2) Contracting the tensor of the highest order form; (3) Resultants and gcds. Each of these methods has strengths and weaknesses, which are summarized in Table 1 jointly with the computational tools required in each case.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics

Method 1 2 3

Advantages Efficient if FN (x) has factors over Q[x, y, z] or if N is odd. Efficient, fast if C 6= k · Id or c 6= 0 Works in the general case. Applicable if FN (x) = (x2 + y 2 + z 2 )p

Disadvantages We may need to factor over the complex numbers. Not applicable if FN (x) = (x2 + y 2 + z 2 )p . Not applicable if C = k · Id or c = 0. In particular, not applicable if FN (x) = (x2 + y 2 + z 2 )p . Works under a generic choice of a point P ∈ S. We may need to test several solutions.

Computational tools Factoring or resultants and gcds (ad-hoc method of Section 2.1.1) Linear Algebra Resultants and gcds

TABLE 1: Advantages, disadvantages and computational tools required for each of the three different computational methods for finding the axis of revolution of a surface of revolution. Method 1=“Factoring”; Method 2=“Tensors”; Method 3=”Resultants and gcds”.

12

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 13

Based on Table 1, we propose the following algorithm to determine the axis of revolution A of a surface of revolution S defined by an implicit algebraic equation F (x) = 0 with highest degree form FN (x). Algorithm Axis of Revolution Input: A square free polynomial F (x) depending explicitly on x, y, z, defining a surface of revolution S. Output: The axis of revolution A of the surface S. If FN (x) 6= (x2 + y 2 + z 2 )p (in particular, if N odd), then: Apply Subalgorithm 2 (Tensors). If Subalgorithm 2 fails, then apply Subalgorithm 1 (Factoring into quadrics). If FN (x) = (x2 +y 2 +z 2 )p or the previous steps did not solve the problem, apply Subalgorithm 3 (Resultants and gcds).

4

C OMPARISON METHOD.

WITH

V RSEK

AND

L AVICKA ’ S

The method in [20] reduces to solving a homogeneous ¨ linear system with six unknowns (the Plucker coordinates of the axis of revolution), and directly provides both the direction of the axis of revolution, and a point on the axis of revolution. In our case, only Subalgorithm 3 simultaneously provides both the direction and a point of the axis of revolution. When Subalgorithm 2 or Subalgorithm 1 are invoked, we first compute the direction of the axis of revolution, and then we separately compute a point on the axis of revolution. From a computational point of view, and although our method exhibits a very good performance as well, in a generic situation we observe better timings for the method in [20]; especially when Subalgorithm 3 needs to be called. This comparison is understandable, since solving linear systems is very easy and fast. Compared to this approach, our algorithm can perform better only when the direction of the axis of revolution is derived by means of Subalgorithm 2 or, in some cases, Subalgorithm 1, and a point on the axis does not need to be computed. For instance, if the polynomial defining the surface is homogeneous, then the axis of revolution, if any, must pass through the origin, so a point on the axis of revolution is known from the beginning. Our algorithm is also especially advantageous when the form of highest order is simple enough to allow us to quickly derive the direction of the axis of revolution, even if afterwards we need to compute a point on the axis. We illustrate the above statements with four somewhat typical examples. The examples have been computed with the computer algebra system Maple 18, run in a personal computer revving up to 2.90 GHz, with 8 Gb of RAM. (1) Let F (x, y, z) = z 15 + z 5 (x2 + y 2 ) − z 3 . According to Corollary 3 and Theorem 6, F (x, y, z) = 0 defines a

surface of revolution about the z-axis. By applying the orthogonal change of coordinates x → A · x + b, where   1/3 2/3 2/3 A = 2/3 −2/3 1/3  2/3 1/3 −2/3 and b = [1, −1, 2]t , we get a new surface of revolution. This new surface is defined by an implicit equation G(x, y, z) = 0 where G(x, y, z) is a dense polynomial of degree 15 and 816 terms, whose infinity norm (i.e. the absolute value of the maximum of the coefficients of G) is close to 1014 . Our method correctly computes the axis of revolution in 0.234 seconds. The algorithm in [20] also correctly computes the axis of revolution in 0.110 seconds. (2) Let F (x, y, z) = z 9 (x2 + y 2 )2 (x2 + y 2 − z 2 )4 . This expression represents a homogeneous polynomial of degree 21, which according to Theorem 6 defines a surface of revolution about the z-axis. Since F is homogeneous, we know that the axis of revolution passes through the origin. Hence, we need only compute the direction of the axis of revolution. After expanding the polynomial (which has 25 terms and an infinity norm of 60), our algorithm finds the axis of revolution after 0.016 seconds. The algorithm in [20] also finds the axis of revolution after 0.047 seconds. (3) Let F (x, y, z) be the polynomial in Example 5 (see Section 2.3). Unlike the two previous examples, in this case we need to invoke Subalgorithm 3. The execution of our algorithm takes 0.281 seconds. The method in [20] takes 0.016 seconds. (4) Let F (x, y, z) = z 17 + z 7 (x2 + y 2 ) − z 3 . By Corollary 3 and Theorem 6, F (x, y, z) = 0 defines a surface of revolution about the z-axis. After applying the orthogonal change of coordinates x → A · x + b, where √  √ √2/2 √2/2 0 A =  2/2 − 2/2 0 0 0 1 and b = [1, −1, 2]t , we get a new surface defined by an implicit equation G(x, y, z) = 0. The polynomial G(x, y, z) has 42 terms and degree 17, and its infinity norm is close to 107 . The highest order form of G(x, y, z) is z 17 , so one immediately has that the axis of revolution is parallel to the z-axis. Furthermore, the intersection of G(x, y, z) = 0 with the plane z = 1 yields √ 2187x2 + 2187(y + 2)2 + 129140136 = 0, which corresponds √ to a (complex) circle centered at the point (0, − 2, 1)√lying on the plane z = 1. Hence, the point (0, − 2, 1) belongs to the axis of revolution. The whole computation takes 0.016 seconds. The method in [20] provides the same result after 0.047 seconds.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 14

ACKNOWLEDGMENTS We would like to thank Bernard Mourrain for discussing some aspects of the problem with us. We would also like to thank Jens Gravesen for some insightful conversations concerning tensors. Juan G. Alc´azar is supported by the Spanish Ministerio de Econom´ıa y Competitividad and by the European Regional Development Fund (ERDF), under the project MTM2014-54141-P, and is a member of the Research Group ASYNACS (Ref. CCEE 2011/ R 34). Ron Goldman is partially supported by a Giner de los R´ıos grant from the Universidad de Alcal´a.

[21] Walker R.J. (1978), Algebraic curves, Springer-Verlag. [22] Wong K.Y., Mendona P., Cipolla R. (2002), Reconstruction of surfaces of revolution from single uncalibrated views, Image and Vision Computing Vol. 22, Issue 10, pp. 829–836.

R EFERENCES [1]

[2] [3] [4] [5] [6]

[7] [8]

[9] [10] [11] [12] [13] [14] [15]

[16] [17] [18] [19] [20]

Andradas C., Recio T., Sendra J.R., Tabera L.F., Villarino C. (2014), Reparametrizing Swung Surfaces over the Reals, Applicable Algebra in Engineering, Communication and Computing, Vol. 25, pp. 3965. Carlini E. (2006), Reducing the number of variables of a polynomial, Algebraic Geometry and Geometric Modeling, Mathematics and Visualization 2006, pp. 237-247, Springer-Verlag. Chionh E.W. (2008), Shifting planes always implicitize a surface of revolution, Computer Aided Geometric Design Vol. 26, Issue 4, pp. 369-377. Cox D., Little J., O’Shea (1992), Ideals, varieties and algorithms, Springer. Comon P., Golub G., Lim L., Mourrain B. (2008), Symmetric tensors and symmetric tensor rank, Siam J. Matrix Anal. App., Vol. 30, No. 3, pp. 1254–1279. Corless R., Galligo A., Kotsireas I., Watt S. (2002), A GeometricNumeric Algorithm for Absolute Factorization of Multivariate Polynomials, Proceedings ISSAC 2002, pp. 37–45, ACM New York, NY, USA. Farouki R.T. (2008), Pythagorean-Hodograph Curves, SpringerVerlag. Fioravanti M., Gonzalez-Vega L., Necula I. (2006), On the intersection with revolution and canal surfaces. Algebraic Geometry and Geometric Modeling, Mathematics and Visualization, pp. 169-184, Springer-Verlag. Gao S. (2003), Factoring multivariate polynomials via partial differential equations, Mathematics of Computation Vol. 72, pp. 801–822. Goldman R. (1983), Quadrics of Revolution, IEEE Computer Graphics and Applications, Vol. 3, No. 3, pp. 68-76. Goldman R. (2009), An Integrated Introduction to Computer Graphics and Geometric Modeling, CRC Press, Taylor and Francis, New York. Jinyuan J., Baciub G., Kwokb K.W. (2004), Quadric decomposition for computing the intersections of surfaces of revolution, Graphical Models, Vol. 66, Issue 5, pp. 303-330. Kaltofen E., May J.P., Yang Z., Zhic L. (2008), Approximate factorization of multivariate polynomials using singular value decomposition, Journal of Symbolic Computation, Vol. 43, Issue 5, pp. 359–376. Pottmann H., Randrup T. (1998), Rotational and Helical Surface Approximation for Reverse Engineering, Computing Vol. 60, pp. 307– 322. Qian X., Huang X. (2004), Reconstruction of surfaces of revolution with partial sampling, Journal of Computational and Applied Mathematics Vol. 163, Issue 1, pp. 211-217. Proceedings of the International Symposium on Computational Mathematics and Applications San Segundo F., Sendra J.R. (2011), Offsetting Revolution Surfaces, Automated Deduction in Geometry, Lecture Notes in Computer Science Vol. 6301, pp. 179-188. ¨ Shalaby M., Juttler B. (2008), Approximate Implicitization of Space Curves and of Surfaces of Revolution, Geometric Modeling and Algebraic Geometry, Springer, pp. 215–227. Shi X., Goldman R. (2012), Implicitizing rational surfaces of revolution using µ-bases, Computer Aided Geometric Design, Vol. 29, Issue 6, pp. 348-362 Simmonds, J.G. (1994), A brief on tensor analysis, Undergraduate texts in Mathematics, Springer-Verlag. Vrsek J., Lavicka M. (2015), Determining surfaces of revolution from their implicit equations., Journal of Computational and Applied Mathematics Vol. 290, pp. 125–135.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TVCG.2015.2498602, IEEE Transactions on Visualization and Computer Graphics 15

´ Juan G. Alcazar is Interim Associate Professor at the Department of Physics and Mathematics ´ in Madrid (Spain). of the Universidad de Alcala, He received his B.S. in Mathematics from the Universidad Complutense de Madrid in 1995 and his PhD in Mathematics from the Universidad de Alcala´ in 2007. He joined the faculty at Universidad de Alcala´ in 2007. His interests include computational algebraic geometry, curves and surfaces and applications in Computer Aided Geometric Design. He has worked on offset curves and surfaces, on topological questions on curves, surfaces and families of curves and surfaces, from a computational point of view, and more recently, on symmetry and similarity detection in curves and surfaces.

Ron Goldman is a Professor of Computer Science at Rice University in Houston, Texas. Professor Goldman received his B.S. in Mathematics from the Massachusetts Institute of Technology in 1968 and his M.A. and Ph.D. in Mathematics from Johns Hopkins University in 1973. Professor Goldman’s current research interests lie in the mathematical representation, manipulation, and analysis of shape using computers. His work includes research in computer aided geometric design, solid modeling, computer graphics, polynomials and splines. He is particularly interested in algorithms for polynomial and piecewise polynomial curves and surfaces, as well as in applications of algebraic and differential geometry to geometric modeling. His most recent focus is on the uses of quaternions and Clifford algebras in computer graphics. Dr. Goldman has published over a hundred articles in journals, books, and conference proceedings on these and related topics. He has also published three books on Computer Graphics and Geometric Modeling: Pyramid Algorithms: A Dynamic Programming Approach to Curves and Surfaces for Geometric Modeling, An Integrated Introduction to Computer Graphics and Geometric Modeling, and Rethinking Quaternions: Theory and Computation. Dr. Goldman is currently an Associate Editor of Computer Aided Geometric Design. Before returning to academia, Dr. Goldman worked for ten years in industry solving problems in computer graphics, geometric modeling, and computer aided design. He served as a Mathematician at Manufacturing Data Systems Inc., where he helped to implement one of the first industrial solid modeling systems. Later he worked as a Senior Design Engineer at Ford Motor Company, enhancing the capabilities of their corporate graphics and computer aided design software. From Ford he moved on to Control Data Corporation, where he was a Principal Consultant for the development group devoted to computer aided design and manufacture. His responsibilities included data base design, algorithms, education, acquisitions, and research. Dr. Goldman left Control Data Corporation in 1987 to become an Associate Professor of Computer Science at the University of Waterloo in Ontario, Canada. He joined the faculty at Rice University in Houston, Texas as a Professor of Computer Science in July 1990.

1077-2626 (c) 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

Finding the Axis of Revolution of an Algebraic Surface of Revolution.

We present an algorithm for extracting the axis of revolution from the implicit equation of an algebraic surface of revolution based on three distinct...
566B Sizes 1 Downloads 8 Views