Polygonal decompositions of quadrilateral subdivision meshes

I.P. Ivrissimtzis,
MPI – Informatik, Stuhlsatzenhausweg 85, 66123, Saarbrucken, Germany
ivrissim@mpi-sb.mpg.de

R. Zayer,
MPI – Informatik, Stuhlsatzenhausweg 85, 66123, Saarbrucken, Germany
zayer@mpi-sb.mpg.de

H.-P. Seidel
MPI – Informatik, Stuhlsatzenhausweg 85, 66123, Saarbrucken, Germany
hpseidel@mpi-sb.mpg.de

Abstract

We study a polygonal decomposition of the 1-ring neighborhood of a quadrilateral mesh. This decomposition corresponds to the eigenvectors of a matrix with circulant blocks, thus, it is suitable for the study of subdivision schemes. First, we calculate the extent of the local mesh area we have to consider in order to get a geometrically meaningful decomposition. Then we concentrate on the Catmull-Clark scheme and decompose the 1-ring neighborhood into 2n planar 2n-gons, which under subdivision scheme transform into 4n planar n-gons coming in pairs of coplanar polygons and quadruples of parallel polygons. We calculate the eigenvalues and eigenvectors of the transformations of these configurations showing their relation with the tangent plane and the curvature properties of the subdivision surface. Using direct computations on circulant-block matrices we show how the same eigenvalues can be analytically deduced from the subdivision matrix.

Keywords: evolving polygons; subdivision; Catmull-Clark; circulant-block matrices.

Introduction

The Catmull-Clark [5,7] and the Doo-Sabin [5,7] were the first bivariate subdivision schemes to appear in the literature. They are both based on quadrilateral meshes. The former, based on the bicubic B-spline, produces $C^2$ continuous surfaces in the regular part of the mesh. Around the extraordinary vertices it generally produces $C^1$ continuous surfaces. Because of its various practical applications, a large part of the research effort in subdivision has been devoted to this particular scheme, e.g. [10,18,13].

In the literature of subdivision there are two main tools for the analysis of the smoothness properties of a subdivision scheme, i.e. the generating functions [8] and the spectral analysis of the subdivision matrix. The latter was first proposed in [7] and it was used in [2] to obtain necessary conditions for tangent plane continuity. In [16] the characteristic map was introduced and necessary conditions for $C^1$ continuity were obtained. The technique of the characteristic map has been expanded and generalized in several directions [14,15]. For a survey of the main methods proposed for studying subdivision surfaces see e.g. [20,19].

Motivation

The motivation behind the present paper lies in the observation that the literature on spectral methods for subdivision focuses on the eigenanalysis of the subdivision matrix. At the same time, only little attention has been given to the spectral analysis of the initial input mesh.

In [11] working on 1-ring neighborhoods of triangle meshes, we carried out the details of their decomposition into planar regular polygons, gaining geometric insights on how a subdivision scheme works. For example, we showed that many local properties of a subdivision surface depend on the initial control polygon only, and not on the particular subdivision scheme. We also found that the normals of the subdivision surface can be very sensitive to small perturbations of the initial control points.

In this paper we apply similar techniques to the 1-ring neighborhood of a quadrilateral mesh. The main difference is that we are now able to study geometric properties of the subdivision surface related to second order smoothness. We believe that such a geometric approach helps to make the study of subdivision simultaneously more intuitive and more rigorous. It allows us to explicitly construct the natural quadratic configurations introduced in [17]. Furthermore, it helps us show that in a block of the subdivision matrix corresponding to a frequency, the eigenvectors play a role as significant as the eigenvalues. As in [11], our geometric approach emphasizes the influence of the initial coarse mesh on the properties of the subdivision surface.

Overview of the method

The main obstacle in a direct generalization of the polygonal decomposition in [4,1,11] to the 1-ring neighborhood of a quadrilateral mesh is the non-trivial connectivity. Fig. 1 shows the 1-ring of a 5-valent vertex $O$ and we notice that some vertices are connected with $O$ while others are not. If we decompose the pentagons

$$A=(a_0,a_1,a_2,a_3,a_4),\ \ \ B=(b_0,b_1,b_2,b_3,b_4)$$

separately into regular eigen-pentagons, then the corresponding components will generally lie on intersecting planes and their transformation under a subdivision matrix will soon become intractable from a geometric point of view. Instead, we can see the initial 1-ring as a decagon and decompose it into ten regular eigen-decagons. After one iteration of a subdivision scheme each eigen-decagon further splits into two pentagons. This happens because in any reasonable subdivision scheme the vertex $O$ has different influence on the polygons $A,B$, see Fig. 2. The two pentagons created this way, will generally lie on different but, crucially, parallel planes, allowing us to split the problem into two parts. The first one deals with the transformations along these parallel planes and the second deals with the transformations in the direction perpendicular to the planes.

Figure 1: The 1-ring neighborhood of $O$.

Figure 2: A planar decagon splits into two planar parallel pentagons.

The rest of this paper is organized as follows. In Section 2 we give a brief background on circulant matrices and polygonal decompositions. In Section 3 we discuss how large the decomposed neighborhood should be to give a geometrically meaningful result. In Section 4 we carry out the polygonal decomposition of the 1-ring neighborhood of a quadrilateral mesh for the purpose of studying the Catmull-Clark scheme. The geometric aspects of this polygonal decomposition are further explored in Section 5. In Section 6 we study the subdivision matrix, tailoring its eigenanalysis to fit with the polygonal decomposition.

Background

We define an $n$-gon $P$ as an $n$-tuple of points of ${\bf R}^3$.

$$P=(p_0,p_1,\dots,p_{n-1})$$ (1)

Notice that there is no restriction on the position of the vertices of the polygon which can intersect itself, degenerate, or even collapse to a single point if all the vertices are in the same position. If $P$ is planar we can also write it as an $n$-tuple of complex numbers, i.e. as a vector in ${\bf C}^n$.

Consider the matrix

$$F_{n}= \begin{pmatrix}1&1&\dots&1\cr 1&\omega_{n}&\dots&\omega_{n... ...cr \vdots&\vdots&&\vdots\cr 1&\omega_{n}^{n-1}&\dots&\omega_{n}\cr\end{pmatrix}$$ (2)

where

$$\omega_{n}=e^{(i2\pi /n)}$$ (3)

is the $n$-th primitive root of unity. According to our definition of a polygon the rows of $F_n$

$$\vec{v_0},\vec{v_1},\dots,\vec{v_{n-1}}$$ (4)

can also be seen as regular planar $n$-gons. Fig. 3 shows these polygons for $n=5,6$.

Figure 3: Up: the planar regular pentagons. Down: The regular planar hexagons.

Moreover, as the rows of $F_n$ are linearly independent, they form a basis for ${\bf C}^n$ and thus, each planar $n$-gon can be written as unique linear combination of the vectors in Eq.(4).

Moreover, following [4], every $n$-gon $P$ in ${\bf R}^3$ has at least one decomposition into $n$ planar regular $n$-gons

$$P=z_0\vec{v_0}+z_1\vec{v_1}+\cdots+z_{n-1}\vec{v_{n-1}}$$ (5)

where $z_j$ are complex numbers with the additional property that the pairs of polygons

$$\{ z_j \vec{v_j}, z_{n-j} \vec{v_{n-j}} \},\ \ \ j=1,2,\dots,\lfloor \frac{n-1}{2} \rfloor$$ (6)

are coplanar. The multiplications in Eq.(5) can be seen as taking place in the complex plane, while the addition is in ${\bf R}^3$. Notice that two different pairs of polygons are not necessarily coplanar.

The matrix $F_n$ and the vectors in Eq.(4) are closely related to the circulant matrices. Recall that a circulant matrix is defined as

$$circ(c_0,c_1,\dots,c_{n-1}) = \begin{pmatrix}c_0&c_1&\dots&c_{n-1... ...&c_{n-2}\cr &&\dots&\cr c_2&c_3&\dots&c_1\cr c_1&c_2&\dots&c_0\cr \end{pmatrix}$$ (7)

The following proposition shows that the linear combination in Eq.(5) has the eigenvectors of the circulant matrices as components.

Proposition 2.1   The eigenvectors of any $n\times n$ circulant matrix $C$ are the rows of $F_n$. The corresponding eigenvalues are y

$$\lambda_k=\sum_{j=0}^{n-1} c_{j} e^{i(2\pi jk/n)}$$ (8)

that is, the values of the generating polynomial of $C$

$$p(z)=c_0+c_1z+\cdots+c_{n-1}z^{n-1}$$ (9)

at $\omega_n^k$, $k=0,1,\dots,n-1$

A proof of the proposition and a deep study of the circulant matrices can be found in [6].

Suitable decomposition neighborhoods

In this section we determine the decomposition neighborhood of a control vertex $P$. We are interested in computing local geometric invariants of the subdivision surface at $P$, e.g. the tangent plane and the three natural quadratic configurations [17]. As these geometric invariants describe continuous derivatives of the subdivision surface at $P$, the aim is to find a neighborhood of $P$ that will contain sufficient information to compute these derivatives.

The area of the initial mesh sufficient to determine the continuous derivatives of the subdivision surface at $P$ depends on the support of the subdivision scheme. That is, it depends on the area of the subdivision surface that will be affected by a change in the position of $P$.

To study the support in more detail we separate the mesh vertices into three distinct classes with respect to $P$, i.e. the points inside the support of $P$ (including $P$ itself), the points on the support boundary, and the points outside of the support. For any nicely behaving subdivision scheme, the relationship between $P$ and another vertex $P'$ of the control mesh, is symmetric with respect to the support. That is, $P'$ is inside, on the boundary, or outside the support of $P$, if and only if $P$ is inside, on the boundary, or outside the support of $P'$, corresponding. See Fig. 4 for the regular case of the Catmull-Clark scheme.

If $P'$ is outside the support of $P$, then it has no influence on the position or any of the continuous derivatives of $P$. Interestingly, this is also the case when $P'$ is on the boundary of the support of $P$, see Fig. 4 (right). To see this, notice that by symmetry $P$ is also on the boundary $B'$ of the support of $P'$, and any disk-like neighborhood of $P$ is divided by $B'$ into two semi-disk like neighborhoods $A_{in}$ and $A_{ex}$, one inside and the other outside the support of $P'$. Any change in the position of $P'$ can affect only $A_{in}$. Thus, it can not change any partial derivative at $P$, because it would cause a derivative discontinuity.

Figure 4: The point $P$ is shown in black. The points inside, on the boundary and outside the support of $P$ are shown in red, blue and white, respectively.

With a similar argument, we also have to consider the vertices on the boundary of the support to prove derivative continuity at $P$. As Fig. 4 shows, the number of these boundary points can be significantly larger than the number of interior points, see also [21]. For example, in the regular case of the Catmull-Clark scheme we need 25 points to prove continuity, compared to nine only points for computing continuous derivatives. Moreover, several methods for proving derivative continuity, like the one in [16], need a semi-regular mesh around $P$ and thus require one or two subdivision steps as an initialization process. This way, the number of vertices involved in a proof may be even larger.

In the regular case, computations of the support and its boundary (which can be either polygonal or fractal) can be found in [12]. On the other hand, it is not straightforward to make the same computations for the irregular case. Nevertheless, here we are not interested in describing the support but only in deciding whether a vertex of the initial control mesh is inside, on the boundary, or outside the support of $P$.

For the case of the Catmull-Clark scheme, using arguments similar to [12], we find that the vertices inside the support of $P$ are in its 1-ring neighborhood, while the vertices on the boundary of the support are the boundary vertices of the 2-ring. Therefore, in the following we use a 1-ring decomposition for computing the geometric invariants of the Catmull-Clark scheme.

1-ring decomposition for the Catmull-Clark scheme

As mentioned in earlier, the 1-ring neighborhood of a quad mesh can be written as a $2n$-gon

$$P = (a_0,b_0,a_1,b_1,\dots,a_{n-1},b_{n-1})$$ (10)

An example for $n=5$ is shown in Fig. 1. This $2n$-gon configuration can be decomposed of into $2n$ regular planar $2n$-gons

$$P=z_0\vec{v_0}+z_1\vec{v_1}+\cdots+z_{2n-1}\vec{v_{2n-1}}$$ (11)

The explicit computation of this decomposition can be found in [4].

The main idea for choosing the decomposition in $2n$-gons instead of decomposing separately two $n$-gons is the observation that if

$$z_j\vec{v_j} = (a_0^j,b_0^j,a_1^j,b_1^j, \dots ,a_{n-1}^j,b_{n-1}^j)$$ (12)

with $j=0,1,\dots,2n-1$, is one component of the sum in Eq.(11), then the $n$-gons

$$\begin{tabular}{rcl} $A^j$\ & = & $(a_0^j,a_1^j, \dots ,a_{n-1}^j)$\ \\ $B^j$\ & = & $(b_0^j,b_1^j, \dots ,b_{n-1}^j)$\ \nonumber\\ \end{tabular}$$

are also planar regular, and moreover, the four $n$-gons

$$A^j, B^j, A^{2n-j}, B^{2n-j}\ \ \ \ \ j=1,\dots,n-1$$ (13)

which are initially coplanar, remain planar and parallel throughout the subdivision process.

The latter property becomes clear when we decompose the subdivision process into several stages. The Catmull-Clark rules for the introduction of new points give

$$a_i\rightarrow\frac{3}{8}(a_i+O) + \frac{1}{16}(a_{i-1}+b_{i-1}+b_i+a_{i+1})$$ (14)

$$b_i\rightarrow\frac{1}{4}(a_i+b_i+a_{i+1}+O)$$ (15)

$$O\rightarrow (1-\alpha-\beta)O + \alpha\sum_{i=0}^{n-1}a_i + \beta\sum_{i=0}^{n-1}b_i$$ (16)

where $\alpha,\beta$ are real variables that can be used to fine-tune the performance of the scheme, see Fig. 5.

In Eq.(14,15,16) we use the substitution symbol $\rightarrow$ to avoid heavy use of indices. In all the three equations the left hand side corresponds to points at step $k$ and the right hand side to points at step $k+1$. Eq.(14) can be decomposed in three steps

$$a_i\rightarrow \frac{3}{4}a_i +\frac{1}{8}(a_{i-1}+a_{i+1})$$ (17)

$$a_i\rightarrow \frac{8}{10}a_i +\frac{1}{10}(b_{i-1}+b_i)$$ (18)

$$a_i\rightarrow \frac{5}{8}a_i +\frac{3}{8}O$$ (19)

while Eq.(15) decomposes into

$$b_i\rightarrow \frac{1}{3}(b_i+a_i+a_{i+1})$$ (20)

$$b_i\rightarrow \frac{3}{4}b_i + \frac{1}{4}O$$ (21)

Figure 5: The insertion of new points.

Then we use the linearity of the subdivision process to further split each of the above substitutions into equations where each point has the same superscript, i.e. corresponds to the same component of Eq.(11). For example Eq.(18) can be written

$$\sum_{j=0}^{2n-1}a_i^j \rightarrow \frac{8}{10}\sum_{j=0}^{2n-1}a_i^j +\frac{1}{10}(\sum_{j=0}^{2n-1}b_{i-1}^j + \sum_{j=0}^{2n-1}b_i^j)$$ (22)

and then decomposed into the substitutions

$$a_i^j \rightarrow \frac{8}{10}a_i^j + \frac{1}{10}(b_{i-1}^j+b_i^j)$$ (23)

with $i=0,\dots,n-1,\ j=0,\dots,2n-1$.

With this decomposition of the subdivision step we can prove

Proposition 4.1   $A^j$ and $B^j$ remain planar polygons on parallel planes throughout the subdivision process, for $j=0,1,\dots,2n-1$.

Proof: First we notice that as $A^j,B^j$ are initially coplanar and thus planar and on parallel planes. Then we check that all the five substitutions in Eq.(17,18,19,20,21) retain this property. The Eq.(17) is on the plane of $A^j$. Eq.(18) sends $A^j$ on a plane parallel to $A^j,B^j$ placed at a ratio $\frac{8}{10}:\frac{2}{10}$ between these planes. Similarly, Eq.(20) sends $B^j$ to a plane parallel to $A^j,B^j$ at ratio $\frac{2}{3}:\frac{1}{3}$. Finally, Eq.(19,21) keep these property because they are similarities. $\Box$

Notice that $A^j,B^j$ do not remain coplanar generally. Indeed, we can see that if $O$ is not on the initial plane of $A^j,B^j$, then after one iteration the $A^j,B^j$ are separated into two different but parallel planes. The reason is that the influence of $O$ on $A^j,B^j$ can be described as a similarity with center $O$ and ratio $5/8,3/4$ respectively, and the different ratios send the initially coplanar polygons into two different but parallel planes, see Fig. 2.

Furthermore, with similar arguments we can show that

Proposition 4.2   $A^j$ and $A^{2n-j}$ remain coplanar throughout the subdivision process, for $j=0,1,\dots,2n-1$.

and

Proposition 4.3   The vertex $O$ and the centers of the polygons $A^j,B^j$ remain collinear throughout the subdivision process.

Proof: Initially they are collinear because $A^j,B^j$ have the same center. By checking all the transformations we verify that they remain collinear. $\Box$

We also notice that the relation between the initial $A^j,B^j$ is given by

$$B^j = \omega_{2n}^j A^j,\ \ \ j=1,\dots,n-1$$ (24)

Evolution of the eigencomponents

Proposition 4.1 allows us to split the study of the evolution of the 1-ring neighborhood under the Catmull-Clark scheme into two parts. One dealing with the transformation of the planar regular polygons and the other dealing with changes in the distance between the parallel planes they lie on. We will show that the former is responsible for first order smoothness properties and the latter for quadratic properties of the limit subdivision surface. We carry out the computations in detail, even though many of the formulas are familiar from the standard eigenanalysis of the subdivision matrix.

Tangent plane properties.

In order to proceed,we calculate the exact transformation of each planar polygon at one subdivision step of the Catmull-Clark scheme. From Eq.(17,18,20) we see that there are three circulant transformations involved

$$C_\alpha = circ(\frac{3}{4},\frac{1}{8},0,\dots,0,\frac{1}{8})$$ (25)

$$C_\beta = circ(\frac{1}{2},0,0,\dots,0,\frac{1}{2})$$ (26)

$$C_\gamma = circ(\frac{1}{2},\frac{1}{2},0,\dots,0,0)$$ (27)

Notice that the difference between the last two circulant transformations is rather subtle, being the result of the enumeration of vertices around the polygon.

Let

$$\lambda^\alpha_0,\lambda^\alpha_1,\dots,\lambda^\alpha_{n-1}$$ (28)

$$\lambda^\beta_0, \lambda^\beta_1, \dots, \lambda^\beta_{n-1}$$ (29)

$$\lambda^\gamma_0, \lambda^\gamma_1, \dots, \lambda^\gamma_{n-1}$$ (30)

be the corresponding eigenvalues of $C_\alpha, C_\beta, C_\gamma$. Combining the three substitutions on $A^j$ we get

$$A^j \rightarrow \frac{5}{8}(\frac{4}{5}\lambda^\alpha_j A^j + \frac{1}{5}\lambda^\beta_j B^j)$$ (31)

and similarly

$$B^j \rightarrow \frac{3}{4}(\frac{2}{3}\lambda^\gamma_j A^j + \frac{1}{3}B^j)$$ (32)

Eq.(31,32) combined give a recursive sequence with matrix of coefficients

$$M_j = \begin{pmatrix}\frac{3}{8}+\frac{1}{16}(\omega_{n}^j+\omega... ...}{16}(1+\omega_n^{-j})\cr \frac{1}{4}(1+\omega_n^j) & \frac{1}{4} \end{pmatrix}$$ (33)

The two eigenvalues of $M_j$ are

$$\lambda_{1,2}=\frac{1}{4}+\frac{1}{8}\cos^2j\theta \pm \frac{1}{4}\cos j\theta\sqrt{1+\frac{1}{4}\cos^2 j\theta}$$ (34)

with $\theta=\frac{\pi}{n}$. To find the limit behavior of an initial vector $(A^j,B^j)^T$ we write it in the base of the two eigenvectors of $M_j$. Then, up to a scaling factor, the component corresponding to the largest eigenvalue (if it exists) will give the limit vector while the eigenvalue itself will give its limit speed of shrinkage.

To find the corresponding eigenvectors we have to solve

$$(M_j -\lambda_{1,2} I) \begin{pmatrix}X \cr Y \end{pmatrix} = \begin{pmatrix}0 \cr 0 \end{pmatrix}$$ (35)

where $\lambda_{1,2}$ are the two eigenvalues given in Eq.(34) and $I$ is the $2\times 2$ identity matrix. The matrices $(M_j -\lambda_{1,2} I)$ are singular and we can obtain the solutions by solving one of the two equations. Solving the equation coming from the second row we get

$$\frac{1}{4}(1+\omega_n^j)X + (\frac{1}{4}-\lambda_{1,2})Y = 0$$ (36)

The subspace of each eigenvector is generated by

$$(1, \frac{\frac{1}{4}(1+\omega_n^j)}{\lambda_{1,2}-\frac{1}{4}})^T$$ (37)

To find the initial input vector we notice that $A^j=z_j \vec{v_j}$, and by Eq.(24) $B^j=\omega_{2n}^j z_j\vec{v_j}$, where $\vec{v_j}$ is the $j$th regular $n$-gon. Thus the initial vector is equal to $(1,\omega_{2n}^j)^T$, up to a linear factor. We see that the initial vector $(1,\omega_{2n}^j)^T$ is not on the subspace of the second eigenvector because

$$\left\vert\begin{array}{cc} 1 & \frac{\frac{1}{4}(1+\omega_n^j)}{\lambda_{1,2}-\frac{1}{4}} \\ 1 & \omega_{2n}^j \\ \end{array}\right\vert \neq 0$$ (38)

Therefore the initial vector always has a component on the dominant eigenvector of $M_j$.

From Eq.(34) we can immediately verify

Proposition 5.1   The regular convex components $A^1,B^1,A^{2n-1},B^{2n-1}$ have the largest eigenvalues.

That means that the Catmull-Clark scheme has nice properties regarding tangent plane and $C^1$ continuity. It also means that any possible singularities can be studied as in [11].

Remark 5.1   The components $A^{n-1},B^{n-1},A^{n+1},B^{n+1}$ are also convex regular but they are filtered out by the subdivision scheme in the same way as the non-convex components.

Application: In the limit, the ratio $r$ of $A^j,B^j$ given by $A^j=r B^j$ can be found by dividing the two components of the largest eigenvector of $M_j$. Using Eq.(37,34) this is equal to

$$\frac{\frac{1}{4}(1+\omega_n^j)}{\frac{1}{4}\cos j\theta\sqrt{1+\frac{1}{4}\cos^2 j\theta} + \frac{1}{8}\cos^2 j\theta}$$ (39)

For $n=4$ and $j=1$ we get a ratio equal to $1+i$, showing that $B^1$ is obtained from $A^1$ by scaling by $\sqrt2$ and rotating by $\frac{\pi}{4}$. For large $n$ and $j=1$ the ratio tends towards $4/(\sqrt5 +1)$.

These ratios are of key importance for the behavior of the subdivision scheme. They describe the $n$-valent 1-ring neighborhoods that are invariant under the Catmull-Clark scheme. Finding such invariant neighborhoods is trivial for $n=4$ because of the symmetry properties of the scheme in the regular case. But the answer is not straightforward for the other valences, see Fig. 6.

Figure 6: The ratio between $A^1$ and $B^1$ for $n=4$ (left) and $n=12$ (right).

Elliptic properties

We will refer to the properties of the subdivision scheme related to the distances between $O,A^j,B^j$ as elliptic properties. The existence of a natural parabolic configuration has been pointed out in [17], and its relation to the subdivision matrix has been studied in [3].

Consider the initial decomposition of the 1-ring neighborhood of $O$. As shown in Section 5.1 the components determining the tangent plane at $O$ are the $A^1,A^{2n-1},B^1,B^{2n-1}$. As $A^1,A^{2n-1}$ are convex regular, their sum $A^1+A^{2n-1}$ is the affine image of a regular convex polygon, see for example [9]. In particular it is inscribed on an ellipse $E_a$ with axes given by

$$l_1=\vert w_1 \vert + \vert w_{2n-1} \vert,\ \ \ \ \ l_2=\bigl\vert \vert w_1 \vert - \vert w_{2n-1} \vert\bigr\vert$$ (40)

with $w_1,w_{2n-1}$ given by

$$A^1=w_1\vec{v_1},\ \ \ A^{2n-1}=w_{2n-1}\vec{v_{n-1}}$$ (41)

where $\vec{v_1},\vec{v_{n-1}}$ are the convex regular $n$-gons. Similarly, the sum of $B^1+B^{2n-1}$ is inscribed on an ellipse $E_b$. In fact, initially $E_a=E_b$ as both ellipses are the same with the ellipse inscribing the $2n$-gonal regular convex component of Eq.(11).

As we have seen, after one iteration of the subdivision scheme the two polygons $A^1+A^{2n-1}$ and $B^1+B^{2n-1}$ will be separated into two parallel planes. The following proposition shows that the Catmull-Clark scheme does not rotate any of the two initial inscribing ellipses $E_a,E_b$ but only scales them.

Proposition 5.2   At each subdivision step $A^j \rightarrow r_A^j A^j$ and $B^j \rightarrow r_B^j B^j$ where $r_A^j,r_B^j \in {\bf R}$. In other words, the components $A^j,B^j$ are not rotated.

Proof: The transformation of $A^j,B^j$ is given by the matrix $M_j$ in Eq.(33) while the initial input is $(1,\omega_{2n}^j)^T$. Assume as an inductive hypothesis that the initial input is $(r_1,r_2\omega_{2n}^j)^T$ with $r_1,r_2$ real numbers. Direct computation shows that

$$M_j(r_1,r_2\omega_{2n}^j)^T = (r'_1,r'_2\omega_{2n}^j)^T$$ (42)

where $r'_1,r'_2$ are real numbers. The proof follows by induction. $\Box$

We will say that a subdivision scheme has nice elliptic properties if in the limit the two ellipses $E_a,E_b$ fit into an elliptic paraboloid with center $O$, see Fig. 7.

Figure 7: Left: The ellipses $E_a,E_b$ after several iteration of the subdivision scheme. Right: The elliptic paraboloid fitting the natural elliptic configuration [17].

If $d_a,d_b$ are the distances between $O$ and the planes of $E_a,E_b$, respectively, the elliptic properties of a scheme depend on the speed of shrinkage $d_a,d_b$ and their ratio in the limit. These depend on subdominant eigenvalue and eigenvector of the matrix of coefficients of the linear recursive sequence

$$\begin{pmatrix}1-\alpha-\beta & \alpha & \beta \cr \frac{3}{8} & ... ...c{1}{2} & \frac{1}{8} \cr \frac{1}{4} & \frac{1}{2} & \frac{1}{4} \end{pmatrix}$$ (43)

with $\alpha,\beta$ as in Eq.(16), provided that the initial input, written in the base of the eigenvectors of the matrix, has a component in the subspace of the subdominant eigenvector. In particular the subdominant eigenvalue of the above matrix should be equal to $\lambda_1^2$, with $\lambda_1$ as in Eq.(34), so that the natural elliptic configuration shrinks in the $z$-direction with quadratic speed compared to the its shrinkage along the tangent plane. Then, if $(u_O,u_a,u_b)^T$ is the subdominant eigenvector, in the limit we have

$$\frac{u_b-u_O}{u_a-u_O}=\frac{d_b}{d_a}$$ (44)

and this ratio should also be equal to the square of the scale ratio between the two ellipses which is given in Eq.(39).

Remark 5.2   Notice that there are only two essentially different cases for the initial input, namely (1,1,1) and (0,1,1), depending on whether $O$ is on the plane of $A^1$ or not. In the first case there is no elliptic component, while in the second case it generally exists.

Example: As it has been pointed out by several researchers, the Catmull-Clark scheme does not satisfy all these elliptic properties for all the valences $n$. Nevertheless in the regular case $n=4$ the matrix in Eq.(43) becomes

$$\begin{pmatrix}\frac{9}{16} & \frac{3}{8} & \frac{1}{16} \cr \fra... ...c{1}{2} & \frac{1}{8} \cr \frac{1}{4} & \frac{1}{2} & \frac{1}{4} \end{pmatrix}$$ (45)

with subdominant eigenvalue equal to 1/4 and corresponding eigenvector $(-2,1,4)^T$. Then, the ratio $d_b/d_a$ is equal to 2, while the scale ratio between the ellipses is $\sqrt2$ as we saw in Section 5.1. Thus, we can verify the nice elliptic properties of the Catmull-Clark scheme in the regular case.

Hyperbolic properties

We will refer to the properties of the subdivision scheme related to the components $A^2,A^{2n-2},B^2,B^{2n-2}$ as hyperbolic. Their eigenvalues are smaller than those related to tangent plane continuity, therefore we are interested in the part of these components that is perpendicular to the tangent plane.

Similarly to the elliptic case, the hyperbolic components have to shrink with quadratic speed compared to the tangent plane. Fig. 8 shows the effect of a hyperbolic component on the shape of the subdivision surface.

Figure 8: Left: The eigencomponent $A^2$. Right: If $A^2$ is the only component outside the tangent plane we get a hyperbolic paraboloid.

Application: We notice that the coarse mesh of the basis function has only elliptic component near $O$. This happens because the 1-ring neighborhood of $O$ is planar and thus, the components perpendicular to the tangent plane are zero. Therefore, the behavior of the basis function near $O$ can be misleading regarding quadratic properties of the subdivision scheme.

To draw a simple configuration with only hyperbolic components, we first place $O$ on the plane of $A^1,A^{2n-1},B^1,B^{2n-1}$. This way the elliptic component becomes zero. Then we add to the 1-ring neighborhood a linear combination of $\vec{v_2},\vec{v_{2n-2}}$ lying on a plane perpendicular to the 1-ring neighborhood.

Exact calculation on CB-matrices

The standard approach to the spectral analysis of subdivision schemes is based on block circulant matrices with each block corresponding to a frequency,see e.g. [14,22]. In order to better capture the geometry, we use circulant block matrices with each circulant block corresponding to a pair of polygons. This suits the setting of the polygonal decomposition.

As we have seen in Section 5, the critical eigenvalues and eigenvectors of the subdivision scheme can be computed from small square matrices of dimension 2 or 3. In the literature these eigenvalues are usually calculated from the full subdivision matrix of dimension $2n+1$. The relations between the eigenvalues of the full matrix and the eigenvalues of its block-diagonal components have been studied in [3].

These relations are due to the circulant block structure of a $2n\times 2n$ block of the full subdivision matrix. Here we show analytically how both eigenvalues and eigenvectors of such matrices relate to the eigenvalues and eigenvectors of the blocks. The results are used to show how the addition of constant term blocks affects the eigenvalues and eigenvectors of the full matrix. Furthermore we show the implications of this on the parameter fine tuning.

First, we establish some properties of circulant matrices and square matrices with circulant blocks (CB-matrices).

Proposition 6.1   Let $C$ be an $n\times n$ matrix $C$ such that $c_{ij}=c$ for $0<i,j<n-1$. The only nonzero eigenvalue of $C$ is $n\times c$ with corresponding eigenvector $v$, where all components of $v$ are equal to $1$.

Proof: Direct consequence of Proposition 2.1. $\Box$

Recall that by Proposition 2.1 all the circulant matrices of a certain order $n$ have the same system of eigenvectors independently of the components of the matrix. This fact is of paramount importance in the following proposition.

Proposition 6.2   Let $B$ be a CB-matrix composed of $n\times n$ circulant blocks $C_{ij}$ of size $m\times m$ where $0<i,j<n-1$. If $v$ is a an eigenvector of the circulant matrices of order $m$ then there exist a vector $\alpha^T=(\alpha_0,..,\alpha_{n-1})$ such that $V^T=(\alpha_{0}v,\alpha_{1}v,..,\alpha_{n-1}v)$ is an eigenvector of $B$. Furthermore the vector $\alpha$ is the eigenvector of a certain matrix $E$. Each element $e_{ij}$ of $E$ is the eigenvalues of $C_{ij}$ associated with the vector $v$.

Proof: It suffices to find the vector of coefficients $\alpha$. Suppose that $V^T=(\alpha_{0}v,\alpha_{1}v,..,\alpha_{n-1}v)$ is an eigenvector of B, then there exist $\lambda$ such that $BV=\lambda V$. On the other hand we have

$$BV=B \begin{pmatrix}\alpha_{0}v\cr \vdots\cr \alpha_{n-1}v\cr\end{pmatrix}$$ (46)

This can be rewritten as

$$BV= \begin{pmatrix}C_{0,0}v&\dots&C_{1,n-1}v\cr \vdots&\vdots&\vd... ...nd{pmatrix} \begin{pmatrix}\alpha_{0}\cr \vdots\cr \alpha_{n-1}\cr\end{pmatrix}$$ (47)

Using the previously stated remark that $v$ is an eigenvector of all circulant matrices of order $m$, the above equation turns out to be

$$BV= \begin{pmatrix}e_{0,0}v&\dots&e_{1,n-1}v\cr \vdots&\vdots&\vd... ...nd{pmatrix} \begin{pmatrix}\alpha_{0}\cr \vdots\cr \alpha_{n-1}\cr\end{pmatrix}$$ (48)

In the following we use the dot operation as a shorthand to avoid expanding the hole system. It should be understood as multiplying the blocks of a vector by the coefficient of another vector, i.e.

$$BV=<tex2html_comment_mark>3 E \begin{pmatrix}\alpha_{0}\cr \vdots... ...\alpha_{n-1}\cr \end{pmatrix}. \begin{pmatrix}v\cr \vdots\cr v\cr \end{pmatrix}$$ (49)

From the above equation it is clear that if the vector $\alpha$ is an eigenvector of E then $V$ is an eigenvector of $B$ . Literally $E\alpha=\lambda\alpha$. This leads to

$$BV=\lambda \begin{pmatrix}\alpha_{0}\cr \vdots\cr \alpha_{n-1}\cr \end{pmatrix} . \begin{pmatrix}v\cr \vdots\cr v\cr \end{pmatrix} =\lambda V.$$ (50)

$\Box$

Remark 6.1   From the above proposition, it is worthwhile to note that the eigensystem of the matrix B can be fully determined if all the matrices E are diagonalizable.

Proposition 6.3   Given two block matrices $B$ and $A$ composed of $n\times n$ blocks $B_{ij}$ and $A_{ij}$ of size $m\times m$ where $0<i,j<n-1$ respectively, such that the $B_{ij}$'s are circulant and the $A_{ij}$'s are constant i.e. ( $a_{lk}^{ij}=a_{ij}$ for all $0<l,k<m-1$).

Let $C=B+A$. If remark 6.1 holds, then the eigenvalues of $C$ and the eigenvalues of $B$ are all identical except for $n$ eigenvalues. If $E_{1}$ is the matrix of eigenvalues of the $B_{ij}$'s corresponding to the eigenvector of ones i.e $(1,1,..,1)$ then the non-matching eigenvalues can be found as the eigenvectors of $E_{1}+T$ where the components of the matrix $T$ are given by $t_{ij}=na_{ij}$.

Proof: The result follows immediately from propositions 6.1 , 6.2 and remark 6.1. $\Box$

Application: The full subdivision matrix for the Catmull-Clark scheme has the form

$$S= \begin{array}{\vert c\vert ccc\vert ccc\vert} \hline * &... ... * & & C_3 & & & C_4 & \\ & & & & & & \\ \hline \end{array}$$

where $C_1,C_2,C_3,C_4$ are circulant matrices. Its eigenvalues are the roots of its characteristic polynomial

$$det(S-xI)$$ (51)

To take out the trivial eigenvalue 1, we add all the columns to the first column, making it a constant column with all the entries equal to $1-x$. Then we subtract the first row from all the other rows. Then, the characteristic polynomial of $S$ is the same as the characteristic polynomial of

$$S'= \begin{array}{\vert c\vert ccc\vert ccc\vert} \hline 1 ... ...C_3 - A & & & C_4 - B & \\ & & & & & & \\ \hline \end{array}$$

where $A,B$ are constant blocks i.e., $a_{ij}=\alpha$ and $b_{ij}=\beta$.

From Proposition 6.3 and the form of the matrix $S'$, we see that the fine tuning parameters $\alpha,\beta$ in Eq.(16) affect only two eigenvalues of the subdivision matrix. Hence, they affect only elliptic properties of the Catmull-Clark without interfering with tangent plane or hyperbolic properties.

Furthermore, this property holds regardless of any a priori knowledge of the underlying algebraic or geometric multiplicity of the subdivision matrix. In the case where the matrix is diagonalizable, our method outlined in proposition 5.3 allows for a direct computation of the eigen basis. In our elaboration throughout the paper, we restricted the range of our decomposition analysis to a one ring neighborhood. We showed that this is sufficient for the computation of all the continuous derivatives. Moreover, increasing the range might only increase the rank deficiency of the subdivision matrix by introducing null columns.

Conclusion - Future Work

We constructed an explicit decomposition of the 1-ring neighborhood of a quadrilateral mesh into planar regular polygonal components which are related to the eigenproperties of the subdivision matrix. We believe that this technique leads to a simultaneously more intuitive and more rigorous study of the subdivision surfaces and we highlighted many aspects of subdivision which are obscured in the literature because of the lack of such a geometric construction.

In the future we plan to use the technique presented here to study more systematically artifacts on subdivision surfaces. We also plan to develop tools for the design of initial control polyhedra that will give subdivision surfaces with prescribed properties.

Bibliography

1

F. Bachmann and E. Schmidt.
n-Ecke.
B.I.-Hochschultaschenbücher, 1970.

2

A.A. Ball and D.J.T. Storry.
Conditions for tangent plane continuity over recursively generated B-spline surfaces.
ACM Transactions on Graphics, 7(2):83-102, 1988.

3

L. Barthe, C. Gerot, M. Sabin, and L. Kobbelt.
Simple computation of the eigencomponents of a subdivision matrix in the fourier domain.
In N. Dodgson, M. Floater, and M. Sabin, editors, Proc. of MINGLE Workshop 2003, page (to appear). Springer, 2004.

4

E.R. Berlekamp, E.N. Gilbert, and F.W. Sinden.
A polygon problem.
Am. Math. Mon., 72:233-241, 1965.

5

E. Catmull and J. Clark.
Recursively generated B-spline surfaces on arbitrary topological meshes.
Computer-Aided Design, 10:350-355, 1978.

6

P.J. Davis.
Circulant matrices.
Wiley-Interscience., 1979.

7

D. Doo and M. Sabin.
Behaviour of recursive division surfaces near extraordinary points.
Computer-Aided Design, 10:356-360, 1978.

8

N. Dyn.
Subdivision schemes in computer-aided geometric design.
In W. Light, editor, Advances in numerical analysis, volume 2, pages 36-104. Clarendon Press, 1992.

9

J.C. Fisher, D. Ruoff, and J. Shilleto.
Perpendicular polygons.
Am. Math. Mon., 92:23-37, 1985.

10

Mark Halstead, Michael Kass, and Tony DeRose.
Efficient, fair interpolation using Catmull-Clark surfaces.
In SIGGRAPH 93, Conference Proceedings, pages 35-44, 1993.

11

I. Ivrissimtzis and H.-P. Seidel.
Evolutions of polygons in the study of subdivision surfaces.
Computing, page (to appear), 2003.

12

Ioannis Ivrissimtzis, Malcolm Sabin, and Neil Dodgson.
On the support of recursive subdivision.
ACM Transactions on Graphics, 23(4):1043-1060, 2004.

13

Ahmad H. Nasri and A. Abbas.
Designing Catmull-Clark subdivision surfaces with curve interpolation constraints.
Computers and Graphics, 26(3):393-400, 2002.

14

J. Peters and U. Reif.
Analysis of algorithms generalizing $B$-spline subdivision.
SIAM Journal on Numerical Analysis, 35(2):728-748, 1998.

15

H. Prautzsch and U. Reif.
Necessary conditions for subdivision surfaces.
Advances in Computational Mathematics, 10:209-217, 1999.

16

U. Reif.
A unified approach to subdivision algorithms near extraordinary vertices.
Computer Aided Geometric Design, 12(2):153-174, 1995.

17

M. Sabin.
Eigenanalysis and artifacts of subdivision curves and surfaces.
In A. Iske, E. Quak, and M. Floater, editors, Tutorials on Multiresolution in Geometric Modelling, pages 69-92. Springer, 2002.

18

Jos Stam.
Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values.
In SIGGRAPH 98, Conference Proceedings, pages 395-404, 1998.

19

J. Warren and H. Weimer.
Subdivision methods for geometric design : a constructive approach.
Kaufmann, 2002.

20

D. Zorin.
$C^k$ Continuity of Subdivision Surfaces.
PhD thesis, California Institute of Technology, Department of Computer Science, 1996.

21

D. Zorin and P. Schröder.
Siggraph 00 course notes, subdivision for modeling and animation., 2000.

22

Denis Zorin.
A method for analysis of ${C^1}$-continuity of subdivision surfaces.
SIAM Journal on Numerical Analysis, 37(5):1677-1708, October 2000.