# A Modified de Casteljau Subdivision
that Supports Smooth Stitching with
Hierarchically Organized Bicubic Bézier Patches

Saied Zarrinmehr
Computer Science and Engineering Department,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Ergun Akleman
Visual Computing & Computational Media,
Texas A&M University, College Station, TX, 77831
<EMAIL_ADDRESS>
Jianer Chen
Computer Science and Engineering,
College of Engineering
Texas A&M University,
College Station, TX, 77831
<EMAIL_ADDRESS>
Joint with Computer Science and Engineering Department.

(October 5, 2023)

###### Abstract

One of the theoretically intriguing problems in computer-aided geometric
modeling comes from the stitching of the tensor product Bézier patches. When
they share an extraordinary vertex, it is not possible to obtain continuity
$C^{1}$ or $G^{1}$ along the edges emanating from that extraordinary vertex.
Unfortunately, this stitching problem cannot be solved by using higher degree
or rational polynomials. In this paper, we present a modified de Casteljau
subdivision algorithm that can provide a solution to this problem. The main
advantage of the modified subdivision is that the continuity $C^{1}$ on a
given boundary edge does not depend on the positions of the control points on
other boundary edges. The modified subdivision allows us to obtain the desired
$C^{1}$ continuity along the edges emanating from the extraordinary vertices
along with the desired $G^{1}$ continuity in the extraordinary vertices (see
Figures 2 and 3). Our modified de Casteljau subdivision, when combined with
topological modeling, provides a framework for interactive real-time modeling
of piecewise smooth manifold meshes with arbitrary topology. Based on this
approach, we have extended an interactive modeling system to topologically
update a smooth surface in real-time by opening or closing holes, creating
handles, and combining and disconnecting surfaces (see Figure 1 for an
example). The smooth surfaces resulting from any given orientable two-manifold
polynomial mesh have the following properties: (1) These surfaces consist of
hierarchically organized $C^{2}$ continuous bicubic Bézier patches. (2) They
provide visual smoothness by providing geometric $G^{1}$ continuity in the
extraordinary vertices and $C^{1}$ continuity along the edges emanating from
the extraordinary vertices. (3) They can be updated interactively in real time
by directly evaluating the underlying Bézier patches.

##  1 Discussion

(a)  Smoothed cube with our method.

(b)  Opening a hole by inserting an edge.

Figure 1: Two images from our interactive modeling session with modified de
Casteljau subdivision. Users can visually experience that inserting an edge
between two opposite faces of the cube opens a $G^{1}$ smooth hole by creating
a 10-sided face.

Tensor product parametric polynomials and rational parametric polynomials are
some of the simplest representations to define surfaces. Using basis functions
such as Bézier or B-spline, parametric polynomial representations of the
tensor product become intuitive to describe the shapes of the final surfaces
Bartels et al. (1987). The coefficients of the Bézier or B-spline functions
are simply considered to be the control points of the final shape of the
surface. Bézier patches are popular since they are particularly intuitive De
Casteljau (1963); Bezier et al. (1974). The resulting surface interpolates
four of the control points, which are the corners of the patch. The
derivatives and second derivatives in the corners of the patch are also
intuitively defined by a set of neighboring control points.

Despite their simplicity and intuitiveness, there are significant problems
with the forms of tensor products in terms of their representational power.
Since their domain is simply a 2D square, they can only represent a very
limited set of topologically distinct two-manifold surfaces. They can
represent only those that can be given as maps from the 2D square domain,
namely genus-0 or genus-1 surfaces. Moreover, with tensor product forms, it is
not possible to include non-quadrilateral faces or non-four-valent vertices,
which are usually called extraordinary vertices. This limitation makes it
harder to model even surfaces of genus-0 or genus-1 Akleman and Chen (2006).

(a)  Detail using method in Akleman et al. (2017a).

(b)  Detail using our method.

Figure 2: Comparison our method with AS&C-method Akleman et al. (2017a) around
a three-valent extraordinary vertex. As it can be seen in this example,
AS&C-method creates broken lines that cause $C^{1}$ discontinuity along the
edges emanating from that extraordinary vertex while our method removes
$C^{1}$ discontinuities by using hierarchically organized Bézier patches that
are constructed by our modified de Casteljau subdivision.

To overcome this problem, one approach has always been to smoothly stitch
quadrilateral patches. Martin Newell’s famous teapot consists of smoothly
stitched bicubic Bézier patches Newell (1975). Smoothly stitched bicubic
Bézier patches were also common during the late 1980s and early 1990s.
Geometric modeling practitioners in the special effects and animation industry
used smooth stitched quadrilateral patches to create complicated
surfaces111Private communications and observations of professional modelers
who worked in special effects and computer animation companies during the
1990s.. Visual smoothness along the patch boundaries was achieved primarily by
manually moving the control vertices.

It is fascinating that professional modelers still manually obtained nice-
looking shapes in that time with this approach despite this being a laborious
and time-consuming process that is prone to errors. Because it is not
theoretically possible to obtain derivative continuity around the
extraordinary vertices, the success of early modelers poses an intriguing and
captivating question. A recent paper by Akleman, Srinivasan, and Chen (AS&C)
Akleman et al. (2017a) has contextualized the ground for this discussion and
provided insight into this question Peters (2017); Akleman et al. (2017b).

AS&C demonstrated that it is easy to construct bicubic Bézier patch control
polyhedra that appear to be stitched with $G^{1}$ continuity regardless of the
underlying mesh topology Akleman et al. (2017a). AS&C-method provided better
visual smoothness than the Catmull-Clark subdivision Catmull and Clark (1978)
for extraordinary vertices of very high valence by providing $G^{1}$
continuity at these points. However, subsequent discussion about this issue
revealed that regardless of how we position control vertices, it is still not
possible to obtain $C^{1}$ or $G^{1}$ continuity along the edges emanating
from that extraordinary vertex Peters (2017); Akleman et al. (2017b). On the
other hand, it turned out that their process causes only subtle
discontinuities. Moreover, these subtle discontinuities of $G^{1}$ occur on
continuous curves $G^{2}$, although it is not often even possible to see the
problem from the silhouette of the shape.

(a)  A detailed view around a 10-valent extraordinary vertex obtained with the
method in Akleman et al. (2017a).

(b)  The same view of the same neighborhood was obtained with our method.

(c)  Another detailed view around the same extraordinary vertex obtained with
AS&C method Akleman et al. (2017a).

(d)  The same view of the same neighborhood was obtained with our method.

Figure 3: Comparison our method with AS&C-method Akleman et al. (2017a) around
the complicated neighborhood of an extraordinary vertex with 10-valent. This
extraordinary vertex comes from a ten-sided face that is obtained by combining
two quadrilateral faces by inserting a single edge shown in Figure 1(b). The
insert edge operation creates the hole shown here. As can be seen in this
example, our method removes the discontinuity $C^{1}$ along the boundary edges
emanating from the extraordinary vertex.

These observations provide a possible explanation for the early popularity of
stitched bicubic Bézier patches among practitioners. Most likely, the process
was not as laborious and time consuming for professional modelers. They might
have obtained visually smooth-looking stitched surfaces with relative ease.
But once modeling tools for Catmull-Clark subdivision Catmull and Clark (1978)
were introduced DeRose et al. (1998) and an exact evaluation method for
Catmull-Clark subdivision was developed for any arbitrary point Stam (1998),
most practitioners could have converted to subdivision surfaces. Moreover, in
modeling practice, people do not deal with very high valence vertices, which
is where smoothness problems with the Catmull-Clark subdivision become more
apparent Akleman et al. (2017a).

##  2 Introduction and Motivation

This discussion establishes the need for a solution to obtain continuity
$C^{1}$ or $G^{1}$ along the edges emanating from that extraordinary vertex
for smooth stitched patches. Although AS&C method Akleman et al. (2017a)
provides a simple solution for interactive applications, the resulting
surfaces cannot be used in applications in which at least $C^{1}$ continuity
should be guaranteed everywhere. Since tensor-product polynomials or rational
polynomials are proven not to have a solution, there is a need to replace
polynomials with analytic functions around extraordinary vertices. Since
analytic functions can be represented as an infinite series of polynomials,
the need for analytic functions suggests that we need a subdivision approach
to solve this problem.

In this paper, we present a new subdivision approach based on a minor
modification of the original de Casteljau subdivision De Casteljau (1963). A
notable property of this new subdivision is that each patch can be subdivided
independently from other patches once the initial control polyhedra are
correctly established. The only problem with this method is that the number of
patches exponentially increases with the number of iterative applications of
the de Casteljau subdivision.

An important property of our modified subdivision is that it allows the number
of patches to be reduced. Assume that after one application of the modified de
Casteljau subdivision, we can evaluate each of the four resulting patches with
the standard bicubic Bézier formulation. Although the four resulting patches
do not form the same shape as the original bicubic Bézier, the whole structure
is still $C^{2}$ continuous everywhere. This suggests that if none of the
corners of a patch is extraordinary, we can simply evaluate it with bicubic
Bézier formulation.

Using this property, we have developed a hierarchical evaluation process that
significantly reduces the number of patches by providing high resolution
around extraordinary vertices, as shown in Figures 2 and 3). This hierarchical
evaluation process is similar to Jos Stam’s exact evaluation of the Catmull-
Clark subdivision at any arbitrary point Stam (1998).

To explain our hierarchical evaluation process, we need to provide a few
definitions. Let a Bézier patch that does not share any extraordinary vertex
with any other Bézier patch be called a regular patch and let the rest be
called extraordinary patches. If a Bézier patch is regular, we can simply
render it using standard Bézier evaluation or original de Casteljau
subdivision.

On the other hand, if the patch is not regular, we apply a modified de
Casteljau subdivision that preserves $C^{1}$ continuity in the edge
boundaries. The modified de Casteljau algorithm is essentially similar to the
original de Casteljau algorithm. It subdivides a $4\times 4$ control
polyhedron of the extraordinary patch into four new $4\times 4$ control
polyhedra by producing four new bicubic Bézier patches that stitch together
with $C^{2}$ continuity. Moreover, like the original subdivision, the modified
subdivision preserves the desired continuity $G^{1}$ in the boundary vertices,
which is given by a position and a tangent plane for each boundary vertex.

The main advantage of the modified subdivision is that the continuity $C^{1}$
in a given boundary edge does not depend on the positions of the control
points on other boundary edges. This allows us to obtain the desired
continuity $C^{1}$ to properly stitch two patches that share an extraordinary
point.

In each iteration of this process, some of the new patches become regular
patches, and we can simply evaluate them using the standard B’ezier
evaluation. For the rest of the extraordinary patches, we continue to apply
the modified de Casteljau subdivision. This process produces hierarchically
organized bicubic Bézier patches that are stitched with $G^{1}$ continuity in
extraordinary vertices and $C^{1}$ continuity along the edges emanating from
extraordinary vertices.

###  2.1 Practical Contributions

There are several practical contributions that can originate from this
theoretical framework. We have identified three main practical contributions,
such as fast computation, greater control of the shapes, and two-manifold mesh
modeling.

#### Fast Computation:

The modified subdivision operates only on a single Bézier control polyhedron.
Therefore, each patch can be evaluated independently of the other patches.
Once one of the subdivided patches is turned into a regular patch, it can be
evaluated directly using the Bézier formulation, allowing real-time rendering.

#### Greater Control of the Shape:

We can provide greater control of the final shape, since with our algorithm it
is possible to provide control on the local region around any extraordinary
vertex using a planar polygon. The number of sides of the planar polygon must
be the same as the valence of the extraordinary vertex. The center and normal
of the planar polygon define the positions and tangents of extraordinary
vertices, similar to controlling the B’ezier surfaces.

#### Smooth two-manifold Modeling:

Practically the key attribute of this method is that it can be used to convert
any quad mesh into a piecewise smooth surface by converting every quad into a
Bézier patch. Since any manifold mesh can be converted to a quadrilateral mesh
using a single application of vertex insertion, which is the remeshing
algorithm of Catmull-Clark subdivision Catmull and Clark (1978), or $\sqrt{2}$
subdivision Li et al. (2004), which is the dual of the simplest subdivision
Peters and Reif (1997), this process allows interactive smooth shape modeling
along with topology changes. In this paper, for interactive smooth surface
modeling, we used vertex insertion in our implementation.

###  2.2 Organization of the Rest of the Paper

The rest of the paper is organized as follows. We first give a subdivision-
based analysis of the root cause of the $C^{1}$ or $G^{1}$ discontinuity
problem along the edges emanating from the extraordinary vertices in Section
3. We later provide a modified de Casteljau subdivision that solves the
discontinuity problem in $C^{1}$ or $G^{1}$ along the edges emanating from the
corresponding extraordinary vertices in Section 4. In Section 5, we provide a
high-level summary of implementation. Finally, in Section 6, we discuss
implementation details and possible extensions.

##  3 Preliminaries

To provide a simple explanation for why Bézier patches cannot provide the
continuity $C^{1}$ or $G^{1}$ along the edges, we revisit both the topology
and geometry of the de Casteljau algorithm. Although most of the results
covered in this section are known, the analysis of the problem as a
preservation of properties through subdivision, as far as we know, is new.
This new method is critical for the formal development of the modified de
Casteljau algorithm. For other related analysis methods, please refer to
Peters (2017).

The control polyhedron of any tensor product Bézier surface is always a
$N\times M$ polyhedral grid, where $N$ and $M$ can be any integer. Figure 4
shows what happens topologically with an application of the de Casteljau
subdivision - without loss of generality - on the $4\times 4$ polyhedral grid.
The classical de Casteljau algorithm always creates four copies of the
original $N\times M$ polyhedral grid. Each of the new polyhedral grids defines
the same Bézier equation. Therefore, we can use the de Casteljau subdivision
to analyze the behavior of the Bézier surfaces.

(a)  Elements of $4\times 4$ control polyhedron.

(b)  Application of the De-Casteljau subdivision.

Figure 4: The de Casteljau subdivision splits the original control polyhedron
into four topologically identical copies. Here, the control points of the
original $4\times 4$ control polyhedron are shown in blue and the control
points of one of the copies are shown with yellow control points. The
positions of the control points of each of these copies produce the original
Bézier formula.

To perform the analysis, we will provide geometric effects of the de Casteljau
subdivision in an unusual form, as shown in Figure 5. In this case, instead of
using the classical recursive definition of de Casteljau, we provided kernels
to produce each new control point for one of the patches. Since in practice we
are interested in bicubic Bézier surfaces, we provide kernels for the $4\times
4$ grid.

(a)  New position of $p_{00}$.

(b)  New position of $p_{01}$.

(c)  New position of $p_{02}$.

(d)  New position of $p_{03}$.

(e)  New position of $p_{10}$.

(f)  New position of $p_{11}$.

(g)  New position of $p_{12}$.

(h)  New position of $p_{13}$.

(i)  New position of $p_{20}$.

(j)  New position of $p_{21}$.

(k)  New position of $p_{22}$.

(l)  New position of $p_{23}$.

(m)  New position of $p_{30}$.

(n)  New position of $p_{31}$.

(o)  New position of $p_{32}$.

(p)  New position of $p_{33}$.

Figure 5: This figure shows the De-Casteljau kernels to produce each control.

Since the algorithm is symmetric, the kernels for other points can be obtained
directly by simply renaming control points. The following equation shows how
the kernel is used to compute the new $P^{\prime}_{mn}$’s (that is, the
positions of yellow points shown in Figure 4(b)) from the old
$P^{\prime}_{ij}$ (that is, the positions of blue points shown in Figure
4(a)).

$P^{\prime}_{mn}=\displaystyle\sum_{i=0}^{3}\sum_{j=0}^{3}\left(\frac{w_{mnij}}{\sum\sum
w_{mnij}}\right)P_{ij},$

where $w_{mnij}$’s are the weights given in Figure 5.

Now, let two Bézier patches be given by a set of points $P^{0}_{ij}$ and
$P^{1}_{ij}$. Using classical Bézier formulation, we can show that these two
Bézier patches that share a boundary edge and have the same $C^{1}$ continuity
in the boundary if the following two conditions hold for all $i$’s:

$\displaystyle P^{0}_{i0}$ $\displaystyle=$ $\displaystyle P^{1}_{i0}$ (1)
$\displaystyle P^{0}_{i1}-P^{0}_{i0}$ $\displaystyle=$ $\displaystyle
P^{1}_{i0}-P^{1}_{i1}$ (2)

These two conditions are easy to obtain if there are only two patches. We can
randomly select all $P^{0}_{i1}$ and $P^{1}_{i1}$ and compute
$P^{0}_{i0}=P^{1}_{i0}$ for all $i$’s as follows:

$P^{0}_{i0}=P^{1}_{i0}=\frac{P^{0}_{i1}+P^{1}_{i1}}{2}$

###  3.1 C1 Condition in Extraordinary Vertex

Now, let $K$ number of Bézier patches that are defined by a set of control
points, $P^{k}_{ij}$ where $0\leq k<K$, share a vertex with continuity
$C^{1}$. In this case, $P^{k}_{0i}$’s also become boundary vertices, and we
cannot freely choose them. The notation we use also becomes cumbersome for
this extraordinary case. Therefore, in order to do the analysis, we first need
to simplify our notation. Let $P^{k}_{00}$ be the same for all $k$’s and the
origin, and let $\vec{V}_{k}=P^{k}_{11}-P^{k}_{00}$. These vectors form a
polygon with $K$ sides that is not necessarily planar (see Figure 6(a)), which
is defined by $K$ vectors such as
$\vec{V}_{0},\vec{V}_{1},\ldots\vec{V}_{K-1}$ that are denoted by yellow
circles in Figure 6(a). Based on condition 2, for all consecutive vectors
$\vec{V}_{k}$ and $\vec{V}_{k+1}$, a new average vector is calculated as
$\frac{(\vec{V}_{k}+\vec{V}_{k+1})}{2}$, which are denoted by blue circles in
Figure 6(a). These vectors give us the blue points that are control points on
the boundaries between patches (see the red lines in Figure 6(a)). Since each
blue point must be in the middle of the line that connects two yellow points,
the control points still visually appear like the original polygon with $K$
sides as shown in Figure 6(a).

(a)  Planar convex polygon around an extraordinary vertex.

(b)  Application of the first iteration of De-Casteljau subdivision.

Figure 6: An example that shows five patches that share a 5-valent
extraordinary point. Although five initial quadrilaterals that are defined by
20 control points form a convex planar pentagon, the application of the de
Casteljau subdivision turns this convex planar pentagon into a star-shaped
planar decagon (10-sided polygon) by breaking line segments.

According to the condition 2, the average of
$\frac{(\vec{V}_{k}+\vec{V}_{k+1})}{2}$ and
$\frac{(\vec{V}_{k+2}+\vec{V}_{k+3})}{2}$ must be the origin. This gives us
the following condition for all $k$’s:

$\vec{V}_{k}+\vec{V}_{k+1}+\vec{V}_{k+2}+\vec{V}_{k+3}=0$

This set of equations can have a non-zero solution only when $K=4$. This
observation suggests that $C^{1}$ continuity cannot be achieved at any
extraordinary point. Although it is not possible to satisfy conditions 1 and 2
for $K\neq 4$, if four patches share a vertex satisfying these conditions, it
becomes trivial. We can freely choose positions of $P^{0}_{11}$, $P^{1}_{11}$,
$P^{2}_{11}$, and $P^{3}_{11}$. Then $P^{k}_{00}$ is simply chosen as their
average222Choosing these control points freely can cause self-intersections,
but, we can ignore that issue.. This well-known result is practically useless
since this arrangement of control points again produces a $N\times M$
polyhedral grid as control polyhedra for four Bézier patches.

###  3.2 G1 Condition in Extraordinary Vertex

As a second option, we can identify the conditions on the positions of
$P^{k}_{0i}$’s for obtaining $G^{1}$ continuity at extraordinary points when
$K\neq 4$. Using de Casteljau for analysis consists of two steps. We first
need to identify an initial configuration of $P^{k}_{0i}$’s to satisfy the
desired property. Second, we need to demonstrate that this property is
preserved under the de Casteljau subdivision. For the conditions $G^{1}$, the
control points of the neighboring quadrilaterals around an extraordinary point
must belong to the same implicit surface. After the application of the de
Casteljau subdivision, the transformed control points of the neighboring
quadrilaterals around an extraordinary point must belong to the same - but not
necessarily to the original - implicit surface.

In the general case, we need to deal with implicit versions of bilinear
surfaces formed by neighboring quadrilaterals, which is not that easy. A
simple approach is to ensure that they share the same plane. In this case, it
is easy to satisfy the same affine equation by simply moving all points in the
neighboring quadrilaterals to the same plane. In the simple vector notation,
we used around extraordinary vertices, all $K$ vectors as
$\vec{V}_{0},\vec{V}_{1},\ldots\vec{V}_{K-1}$ must be perpendicular to the
same vector. To avoid self-intersections, the origin and the vectors must form
a star shape, that is, any line emanating from the origin must intersect with
the polygon that is formed by $\vec{V}_{0},\vec{V}_{1},\ldots\vec{V}_{K-1}$
only once.

To further simplify the analysis, let
$\vec{V}_{0},\vec{V}_{1},\ldots\vec{V}_{K-1}$ form a convex polygon and let
$P^{k}_{00}$ be the origin (i.e.
$\vec{V}_{0}+\vec{V}_{1}+\ldots+\vec{V}_{K-1}=0$). This can simply be achieved
by providing a planar convex polygon and using its corners as
$P^{k}_{11}$’s.We can still enforce the condition 2 and for all consecutive
vectors $\vec{V}_{k}$ and $\vec{V}_{k+1}$, we can compute a new average vector
as $\frac{(\vec{V}_{k}+\vec{V}_{k+1})}{2},$ which are again denoted by blue
circles in Figure 6(a). This is basically similar to the previous case. The
only difference is that the $K$-sided polygon is now planar.

(a)  $n=3$.

(b)  $n=4$.

(c)  $n=8$.

Figure 7: Examples demonstrate unbroken line property can only be preserved
for $n=4$. This is because the initial structure provides the same bilinear
configuration for four quadrilaterals. Since the de Casteljau algorithm
consists of bilinear transformations, the lines do not break for $n=4$. The
structures of the original control points that form convex polygons $n$ are
shown in lighter color in (b).

This configuration along with the de Casteljau subdivision makes it easy to
evaluate $G^{1}$ continuity around an extraordinary vertex. It also provides
an intuitive understanding of what is happening locally around the
extraordinary vertex. All we need to check is whether the de Casteljau
subdivision preserves the planarity condition of this polygon to confirm the
preservation of the $G^{1}$ condition in an extraordinary vertex. Note that
any linear combination of the vectors
$\vec{V}_{0},\vec{V}_{1},\ldots\vec{V}_{K-1}$ still provides another vector in
the same vector space. de Casteljau subdivision just provides such a linear
combination of these vectors, therefore, the transformed vectors still define
the same vector space. In other words, the planar property is preserved under
the de Casteljau subdivision. The fact that the de Casteljau subdivision is
simply a linear combination of control vertices suffices to demonstrate that
the $G^{1}$ condition in extraordinary vertices is always preserved.

$G^{1}$ continuity is one of the reasons why AS&C method Akleman et al.
(2017a) provides visually smoother results than the Catmull-Clark subdivision
Catmull and Clark (1978). The second reason is that the AS&C method uses
polygons formed by the Doo-Sabin subdivision Doo and Sabin (1978), which is
known to provide planar regular and convex polygons in limit. As we discuss
later, having more regular and convex polygon is useful to obtain visually
desirable results. However, the Doo-Sabin stage may not necessarily be needed,
since regular and convex polygons can still be obtained by methods other than
Doo-Sabin.

###  3.3 C1 or G1 Condition in Boundaries

The next question is to check if $C^{1}$ or $G^{1}$ conditions are preserved
in the boundaries. To preserve the condition $C^{1}$ in the boundaries, the
boundary points must always be the average of two inner points after
subdivision. To provide the $G^{1}$ condition in non-planar boundary regions,
the necessary and sufficient condition is that the line that connects one
boundary point and two inner points must be straight, i.e. the boundary point
must be a Barycentric average of the two inner points Bartels et al. (1987).

In this section, we demonstrate that conditions $C^{1}$ and $G^{1}$ in the
boundaries can never be preserved when we use planar convex polygons to
produce control points. Note that the $C^{1}$ condition requires the blue
points to be in the middle of the line connecting two yellow points. Let us
call this property the unbroken line property. If the property of the unbroken
line is not preserved by the de Casteljau subdivision, we cannot obtain
continuity of $C^{1}$ or $G^{1}$. Unfortunately, the lines start to break in
the planar region, as shown in Figure 6(b). These broken lines propagate along
the edges when we apply the de Casteljau subdivision shown in Figure 5 to
other boundary regions.

(a)  Two patches that share an extraordinary point. This image shows Bézier
control points on both sides of the boundary. Due to the extraordinary point,
$P^{0}_{01}$, $P^{0}_{00}$, $P^{1}_{00}$, and $P^{1}_{01}$ are not on the same
line.

(b)  Control points after the application of de Casteljau subdivision. Note
that de Casteljau splits the original patches into two in the boundary.

(c)  Separated patches to demonstrate the structure of broken lines.

Figure 8: An example demonstrates the change in boundary control points after
one application of the original de Casteljau subdivision around a boundary
shared by two patches.

To demonstrate this problem, we can now apply actual de Casteljau subdivision
kernels to these vectors. The vertex at the origin, of course, never moves. In
the first step of de Casteljau, the original blue points at
$\frac{(\vec{V}_{k}+\vec{V}_{k+1})}{2}$ move to
$\frac{(\vec{V}_{k}+\vec{V}_{k+1})}{4}$. On the other hand, the original
yellow points at $\vec{V}_{k}$ move to
$\frac{\vec{V}_{k}}{2}+\frac{(\vec{V}_{k-1}+v_{k+1})}{8}$ since the new point
is the calculated average of four control points. It is easy to show that the
property of the unbroken line is preserved if and only if $n=4$ (see Figure
7). For all other cases, the process breaks the lines and the initial planar
polygon with $K$ sides becomes a planar polygon with $2K$ sides (see Figures
9, 7(a) and 7(c) ). In other words, the straight line segments will be broken
originally if $K\neq 4$.

The problem is that these broken lines are not contained in the planar region
(see Figure 8). They propagate along the boundary-edge. As can be verified by
using two de Casteljau kernels in Figures 5(j) and 5(n) that are used to
calculate $P^{\prime}_{21}$ and $P^{\prime}_{31}$, all boundary lines in newly
created patches that share the original extraordinary vertex are broken (see
Figure 8(c)). An additional problem is that the opposite end of the newly
created patches that share the 4-valent vertex now also has a broken line
caused by the computation of $P^{\prime}_{31}$’s. As a result, the broken
lines travel further towards the 4-valent vertex in every iteration.

The immediate neighborhood of an extraordinary point is still planar, since
the first set of broken line segments is still on the original plane. However,
the rest of the broken lines define two inconsistent bilinear surfaces with
neighboring control points along the boundary that separates two consecutive
patches. These broken lines cause discontinuity $C^{1}$ and $G^{1}$ along the
boundary between two patches. The problem of broken lines cannot be solved by
using higher-degree or rational polynomials, since the kernels for boundary
points will still be in the same basic form for both higher-degree and
rational cases.

In summary, the broken line segments that are by-products of extraordinary
vertices cannot be contained in the planar region around extraordinary points.
They create nonplanar bilinear regions by traveling along the boundary in
every subsequent step of the de Casteljau algorithm while causing $C^{1}$ and
$G^{1}$ discontinuity along boundary edges. Although the problem is expected
to be worse when $K$, the valence of extraordinary vertex increases,
discontinuity along boundary edges $C^{1}$ and $G^{1}$ is not visible to a
casual viewer during interactive modeling. This is also not unexpected, as
these discontinuities occur along $C^{2}$ continuous curves. This effect is
subtle, since the broken lines gradually become more and more straightened by
a series of interpolations.

(a)  The original de Casteljau around a planar extraordinary region.

(b)  The modified de Casteljau around the planar extraordinary region.

Figure 9: An example that demonstrates the difference between the first steps
of original and modified de Casteljau algorithms visually.

For interactive applications, therefore, AS&C method Akleman et al. (2017a)
can still be used successfully. This algorithm still provides stitched
surfaces that appear to be visually $G^{1}$ continuous. On the other hand, if
the resulting smooth shapes are to be used for simulation, manufacturing, or
3D printing, there is a need to remove $C^{1}$ discontinuities that are
important for the shape quality.

##  4 Modified de Casteljau Subdivision

Although it is not possible to obtain $C^{1}$ continuity along the edges
emanating from extraordinary vertices, we observe that this problem can be
resolved by a variety of minor modifications to the de Casteljau subdivision
that can guarantee the avoidance of broken lines. These are essentially
subdivision algorithms, and we cannot have a closed polynomial form anymore.

(a)  $n=3$.

(b)  $n=4$.

(c)  $n=8$.

Figure 10: Examples demonstrating the first step of modified de Casteljau
subdivision.

Among these solutions, one of them, in particular, appears to be much simpler
and more powerful than others. Notice that the problem of broken lines stems
from the additional term $\frac{(\vec{V}_{k-1}+\vec{V}_{k+1})}{8}$, we simply
remove it and move the original yellow points at $\vec{V}_{k}$, to
$\frac{\vec{V}_{k}}{2}$ (see 9(b)). This simple change can replace the
original bilinear transformation in Figure 5(f) with a uniform scale around
the origin in Figure 11(f) to calculate $P^{\prime}_{11}$.

Some other de Casteljau kernels must also change to avoid line breaks. As
discussed above, the kernels in Figures 5(j) and 5(n), which are used to
calculate $P^{\prime}_{21}$ and $P^{\prime}_{31}$, also produce broken lines.
We want to point out that the problems in these two kernels come from the
interpolation of the bilinear transformation that is used to calculate
$P^{\prime}_{11}$. If we replace the original bilinear term with scaling, we
obtain new kernels for $P^{\prime}_{21}$ and $P^{\prime}_{31}$ as shown in
Figures 11(j) and 11(n). Due to the symmetry, we have to use the same kernels
for $P^{\prime}_{12}$ and $P^{\prime}_{13}$ as shown in Figures 11(g) and
11(h).

We want the boundary curves to stay the same. Therefore, the kernels for
computing $P^{\prime}_{00}$, $P^{\prime}_{01}$, $P^{\prime}_{02}$,
$P^{\prime}_{03}$, $P^{\prime}_{10}$, $P^{\prime}_{20}$ and $P^{\prime}_{30}$
will remain the same. The kernels for $P^{\prime}_{22}$, $P^{\prime}_{23}$,
$P^{\prime}_{32}$, and $P^{\prime}_{33}$ can also remain the same. However,
since we do not use the original de Casteljau subdivision, we do not create
bicubic Bézier Patches anymore. In other words, we no longer get $C^{2}$
continuity. If we keep the original kernels for these four points, the newly
created patches will no longer be merged with continuity $C^{2}$, since the
positions of $P^{\prime}_{12}$, $P^{\prime}_{21}$, $P^{\prime}_{31}$, and
$P^{\prime}_{13}$ are not the same as before. Therefore, to obtain $C^{2}$
again between newly created patches, we need to readjust the kernels
corresponding to $P^{\prime}_{22}$, $P^{\prime}_{23}$, $P^{\prime}_{32}$, and
$P^{\prime}_{33}$.

To compute new kernels for these four points, we again used the original de
Casteljau process, which provides these kernels through an iterative process.
When we replace the boundary kernels in an iterative process, we obtain the
kernels given in Figures 11(k), 11(l), 11(o), and 11(p). All these kernels can
be obtained simply using the iterative process of de Casteljau just by
replacing the original bilinear kernels of $P^{\prime}_{11}$, $P"_{13}$,
$P^{\prime\prime\prime}_{31}$, and $P^{\prime\prime\prime\prime}_{33}$ where
the prime notation refers to four different subpatches of de Casteljau.

Since we keep the boundary kernels the same, this modification does not change
the shapes of the boundary curves. It only changes the partial derivative at
the extraordinary point. Therefore, the difference in positions is really
negligible except for a significant visual improvement in mesh quality, as
shown in Figures 2 and 3. Another useful property of this particular
modification is that we do not have to consider neighborhood patches. Each
patch can still be evaluated independently.

(a)  New position of $p_{00}$.

(b)  New position of $p_{01}$.

(c)  New position of $p_{02}$.

(d)  New position of $p_{03}$.

(e)  New position of $p_{10}$.

(f)  New position of $p_{11}$.

(g)  New position of $p_{12}$.

(h)  New position of $p_{13}$.

(i)  New position of $p_{20}$.

(j)  New position of $p_{21}$.

(k)  New position of $p_{22}$.

(l)  New position of $p_{23}$.

(m)  New position of $p_{30}$.

(n)  New position of $p_{31}$.

(o)  New position of $p_{32}$.

(p)  New position of $p_{33}$.

Figure 11: This figure shows the Modified De-Casteljau kernels to produce each
control point.

##  5 Implementation, Discussion and Results

In this paper, we introduce an extremely simple modification of the de
Casteljau subdivision that can guarantee the production of continuity $C^{1}$
and $G^{1}$ for any extraordinary vertex. However, the modified de Casteljau
subdivision only guarantees that the limit surface is free of discontinuities
$G^{1}$ or $C^{1}$. In any resolution, if we compute Bézier surface with
direct evaluation using higher resolution control points obtained from the de
Casteljau subdivision, there will still be $G^{1}$ discontinuity that could be
visible in resolutions higher than current resolutions. Therefore, it is
better to obtain the desired resolution directly with the modified de
Casteljau subdivision if there is an extraordinary vertex.

The modified de Casteljau can be applied to $n=4$ cases even when the control
points are not in the same plane. However, in this case, the resulting surface
is not a polynomial. Therefore, it is better to apply original de Casteljau
algorithms for regular patches whose corners are shared only by $4$ other
patches. With each application of the de Casteljau subdivision, each of the
original patches will be subdivided into four patches. This process will
populate surfaces with regular bicubic Bézier patches that do not need to be
further subdivided. We have to continue subdivision only around extraordinary
vertices. This process creates hierarchically organized bicubic Bézier patches
as shown in Figures 2 and 3.

(a)  Start with an orientable two-manifold polygonal mesh.

(b)  Assign K-sided polygons to every face, vertex, and edge of the initial
mesh.

(c)  Create Bézier control polyhedra using K-sided polygons.

(d)  Obtain a smooth surface that consist of bicubic Bézier patches.

Figure 12: The steps of our process. In this case, the initial two-manifold
polygonal mesh is a regular tetrahedron.

(a)  Initial two-manifold mesh.

(b)  Assign a polygon to each face, vertex, and edge.

(c)  Create Bézier control polyhedra.

(d)  Obtain a smooth surface with modified de Casteljau.

Figure 13: One more example demonstrating the steps of our process. Note that
in the polygonal rendering, these two models appear the same. Once we turn
them into smooth models, the difference between them becomes visible, and we
can appropriately visualize the edge insertion operation that creates a
10-sided face.

(a)  Initial two-manifold mesh.

(b)  Assign a polygon to each face, vertex, and edge.

(c)  Create Bézier control polyhedra.

(d)  Obtain a smooth surface with modified de Casteljau.

Figure 14: Another example demonstrating the five steps of our process. Note
that in the polygonal rendering, these two models appear the same. Once we
turn them into smooth models, the difference between them becomes visible, and
we can appropriately visualize the edge insertion operation that creates a
10-sided face.

(a) The polygon mesh obtained from surface in Figure 12(d).

(b) The polygon mesh obtained from surface in Figure 13(d).

(c) The polygon mesh obtained from surface in Figure 14(d).

Figure 15: These examples demonstrate the smoothed polygonal mesh obtained
from two-manifold meshes shown in Figures 12, 13, and 13.

###  5.1 High-Level Summary of Our Implementation

Using this theoretical approach, we have developed a simple process to
interactively model complicated smooth surfaces. Our process produces
hierarchically organized bicubic Bézier patches that are stitched with $G^{1}$
continuity in extraordinary vertices and $C^{1}$ continuity along the edges
emanating from extraordinary vertices. A high-level summary of our
implementation approach to producing smooth stitched bicubic Bézier patches
from any given two-manifold mesh can be given as a five-step process. Figures
12, 13, and 13 provide examples that demonstrate our process. It also shows
the superiority of smooth modeling over polygon modeling: Although there is no
visual difference between two polygonal models, since the additional edge is
not visible, smooth models clearly show the hole resulting from the insert
edge operation.

  1. 1.

Start with an orientable two-manifold mesh (see Figure 12(a)).

  2. 2.

Assign a planar convex polygon with $K=n$ sides to every face with $n$, a
planar convex polygon with $K=m$ sides to every $m$ valent vertex, and a
polygon with $K=4$ sides to every edge of the initial two-manifold mesh.
Polygons with the $4$ sides, i.e. quadrilaterals, do not have to be planar.
(see Figure 12(b))

  3. 3.

Using the vertex insertion scheme, which is the remeshing scheme of Catmull-
Clark subdivision (Catmull and Clark, 1978), obtain a quadrilateral
decomposition of the initial two-manifold mesh. Using this quadrilateral
decomposition, construct a $4\times 4$ bicubic Bézier control points in such a
way that every resulting control polyhedron shares the same boundary
conditions with its four neighboring control polyhedra (see Figure 12(c)).

  4. 4.

If the patch is regular, simply render it by evaluating the corresponding
Bézier polynomial. If the patch is not regular, subdivide it into two four
patches by applying the modified de Casteljau subdivision that preserves the
boundary conditions. Repeat the process until we reach the desired resolution
(see Figure 12(d)).

  5. 5.

If desired, convert it back to a two-manifold polygonal mesh (see Figure
15(a)).

(a)  With a star polygon.

(b)  With a convex polygon.

Figure 16: An example of the effect of the shape of the $K$-sided polygons
over the visual quality of extraordinary neighborhoods: A 10-valent
neighborhood obtained by using a $10$-sided (a) star polygon, and (b) convex
and more regular polygon.

We included some additional images that can provide more information about the
features of our system. Uploaded videos provide more information about the
system. These are not essential for this paper. Therefore, we include these
images to provide additional visual information for a clarification of the
advantages of smooth modeling over direct polygonal modeling when dealing with
an unusual topology.

##  6 Discussion, Conclusion and Future Work

The theoretical analysis demonstrated that we can freely choose the sizes,
positions, and shapes of $K\neq 4$-sided polygons as long as they are planar
star polygons. Their shapes and sizes play a role similar to Kochaneck-Bartels
splines, and the decisions about their sizes and shapes belong to shape
designers Kochanek and Bartels (1984). On the other hand, it is still possible
to give some intuitive guidelines similar to Catmull-Rom splines Catmull and
Rom (1974); Yuksel et al. (2011).

We observe that it is better to have an almost-regular convex polyhedron. In
other words, we want
$\vec{V}_{k}\cdot\vec{V}_{k}\approx|\vec{V}_{k}||\vec{V}_{k+1}|\cos(2\pi/K)$
and $|\vec{V}_{k}|\approx|\vec{V}_{k+1}|$ for every $k$. Figure 16 shows the
superiority of such almost regular convex polygons over star polygons. Such
regular polygons can usually be obtained with Doo-Sabin refinement rules. Doo-
Sabin refinement is known to be powerful in obtaining convex, almost regular,
and almost planar polygons with a few iterations Akleman et al. (2017a).
However, we observed that for some of the unusual configurations, it takes
more steps to obtain a reasonably regular polygon. Therefore, we added a
second step that takes the dual polygon to be even a number of times. This
second step helps us to obtain better-looking initial estimations of the
polygons for unusual cases. This stage does not affect continuity $C^{1}$ or
$G^{1}$, but improves the look around the extraordinary vertices of high
valent. For example, in the example given in Figure 16, just a few steps of
the Doo-Sabin refinement process produced a star polygon instead of a convex
one. Just making it convex improves the quality significantly, as shown in
Figure 16.

The positions and sizes of the $K$-sided planar polygons are also important to
obtain good results. For estimating polygon centers that are interpolated by
curves, it is possible to use the Doo-Sabin or Catmull-Clark refinement to
estimate positions and sizes (see Figure 17). As shown in Figure 17, it is
possible to change the overall shape of the surfaces by changing the sizes and
positions of planar polygons on the sides of $K$.

(a)  A bulkier and boxier shape.

(b)  A smoother shape.

Figure 17: Effect of changing the positions and sizes of $K$-sided polygons.

Another aspect of this method is that complicated shapes can be described with
a two-manifold structure that consists of only a few vertices, faces, and
edges. To unleash the true power of this method for designers, there is a need
for the development of user interfaces that can allow shape designers to
control final surfaces more intuitively through $K$-sided polygons.

An additional effect that comes free with this approach is to control the
relative orientations of planar polygons with $K$ sides. Figure 18 shows two
such examples that are obtained by rotating each polygon on the same side $K$.
These shapes cannot be created by any other subdivision method since we
explicitly control iso-parameter curves on the surface.

Figure 18: Effect of rotating every $K$-sided polygon.

In this particular implementation, to construct bicubic Bézier control
polyhedra we used a vertex insertion scheme. This is not necessary. We can
even start with a quad-mesh that is obtained by any quadrangulation method
such as Kälberer et al. (2007); Bommes et al. (2009). It is also possible to
define the polygon with $K$ side only for the vertices and faces and construct
bicubic Bézier control polyhedra using the $\sqrt{2}$ subdivision Li et al.
(2004), that is, the dual of the simplest Peters and Reif (1997). Our initial
investigation suggests that the method we have introduced in this paper can
also be used to obtain $C^{2}$-stitched tensor product Bézier patches. An
application of such $C^{2}$-stitched surfaces can be the approximation of
given shapes using Morse-Smale complexes Ni et al. (2004).

(a)  A bulkier and boxier shape.

(b)  A smoother shape.

Figure 19: Effect of changing the positions and sizes of $K$-sided polygons.

Figures 19 and 21 provide two more examples that demonstrate the effect of
changing the positions and sizes of polygons with sides $K$. These examples
show that by providing more control to designers, it can be possible to
control the shapes of holes and handles. This also suggests that we can simply
animate the geometry by changing these polygons with $K$ sides. Figures 22,
23, and 24 further demonstrate the visual quality resulting from our modified
de Casteljau subdivision compared to the Catmull-Clark subdivision.

(a)  Initial polygonal mesh that is created by connecting two cubes with one
edge, which creates a bridge between two cubes.

(b) Corresponding smooth surface.

Figure 20: Smoothing two cubical shapes that are connected by one edge, which
creates a bridge.

(a)  A bulkier and boxier shape.

(b)  A smoother shape.

Figure 21: Effect of changing the positions and sizes of $K$-sided polygons.

(a)  Catmull-Clark suvbdivision.

(b)  Our modified de Casteljau Subdivision.

Figure 22: Comparison with Catmull-Clark subdivision. Note that we get almost
the same shape with smoother hole.

(a)  Catmull-Clark suvbdivision.

(b)  Our modified de Casteljau Subdivision.

Figure 23: Another comparison with Catmull-Clark subdivision. Note that we get
almost the same shape with smoother handle.

(a)  Catmull-Clark suvbdivision.

(b)  Our modified de Casteljau Subdivision.

Figure 24: Another detailed comparison with Catmull-Clark subdivision. Note
that it is clear that the handle is smoother.

## References

  * Bartels et al. [1987] Richard H Bartels, John C Beatty, and Brian A Barsky.  _An introduction to splines for use in computer graphics and geometric modeling_.  Morgan Kaufmann, 1987.
  * De Casteljau [1963] Paul De Casteljau.  Courbes et surfaces à pôles.  _André Citroën, Automobiles SA, Paris_ , 1963.
  * Bezier et al. [1974] Pierre Bezier et al.  Mathematical and practical possibilities of unisurf.  _Computer Aided Geometric Design_ , 1(1), 1974.
  * Akleman and Chen [2006] Ergun Akleman and Jianer Chen.  Insight for practical subdivision modeling with discrete gauss-bonnet theorem.  In _International Conference on Geometric Modeling and Processing_ , pages 287–298. Springer, 2006.
  * Akleman et al. [2017a] Ergun Akleman, Vinod Srinivasan, and Jianer Chen.  Interactive modeling of smooth manifold meshes with arbitrary topology: G1 stitched bi-cubic bézier patches.  _Computers & Graphics_, 66:64–73, 2017a.
  * Newell [1975] M. E. Newell.  _The Utilization of Procedural Models in Digital Image Synthesis_.  Ph.D. Thesis, UTEC SCc-76-218, University of Utah, Department of Computer Science, Salt Lake City., 1975.
  * Peters [2017] Jorg Peters.  On G1 stitched bi-cubic bezier patches with arbitrary topology.  _Computers & Graphics_, 67:11–13, 2017.
  * Akleman et al. [2017b] Ergun Akleman, Vinod Srinivasan, and Jianer Chen.  Response to “on g1 stitched bi-cubic bézier patches with arbitrary topology”.  _Computers & Graphics_, 2017b.
  * Catmull and Clark [1978] Edwin Catmull and James Clark.  Recursively generated b-spline surfaces on arbitrary topological meshes.  _Computer-aided design_ , 10(6):350–355, 1978\.
  * DeRose et al. [1998] Tony DeRose, Michael Kass, and Tien Truong.  Subdivision surfaces in character animation.  In _Proceedings of the 25th annual conference on Computer graphics and interactive techniques_ , pages 85–94. ACM, 1998.
  * Stam [1998] Jos Stam.  Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values.  In _Proceedings of the 25th annual conference on Computer graphics and interactive techniques_ , pages 395–404. ACM, 1998.
  * Li et al. [2004] Guiqing Li, Weiyin Ma, and Hujun Bao.  $\sqrt{2}$ subdivision for quadrilateral meshes.  _The Visual Computer_ , 20(2-3):180–198, 2004\.
  * Peters and Reif [1997] Jörg Peters and Ulrich Reif.  The simplest subdivision scheme for smoothing polyhedra.  _ACM Transactions on Graphics (TOG)_ , 16(4):420–431, 1997.
  * Doo and Sabin [1978] Daniel Doo and Malcolm Sabin.  Behaviour of recursive division surfaces near extraordinary points.  _Computer-Aided Design_ , 10(6):356–360, 1978\.
  * Kochanek and Bartels [1984] Doris HU Kochanek and Richard H Bartels.  Interpolating splines with local tension, continuity, and bias control.  _Acm Siggraph Computer Graphics_ , 18(3):33–41, 1984.
  * Catmull and Rom [1974] Edwin Catmull and Raphael Rom.  A class of local interpolating splines.  In _Computer aided geometric design_ , pages 317–326. Elsevier, 1974\.
  * Yuksel et al. [2011] Cem Yuksel, Scott Schaefer, and John Keyser.  Parameterization and applications of catmull–rom curves.  _Computer-Aided Design_ , 43(7):747–755, 2011\.
  * Kälberer et al. [2007] Felix Kälberer, Matthias Nieser, and Konrad Polthier.  Quadcover-surface parameterization using branched coverings.  _Computer graphics forum_ , 26(3):375–384, 2007\.
  * Bommes et al. [2009] David Bommes, Henrik Zimmer, and Leif Kobbelt.  Mixed-integer quadrangulation.  _ACM Transactions On Graphics (TOG)_ , 28(3):77, 2009.
  * Ni et al. [2004] Xinlai Ni, Michael Garland, and John C Hart.  Fair morse functions for extracting the topological structure of a surface mesh.  _ACM Transactions on Graphics (TOG)_ , 23(3):613–622, 2004.