166 |
W.A.P. Smith |
Fig. 4.15 Loop subdivision scheme: (a) creation of an edge vertex; (b) calculating new position for original points; (c) triangulation of new edge points
Fig. 4.16 Two iterations of the Loop subdivision scheme applied to an icosahedron base mesh. Figure adapted from [77]
available. The simplest choice is:
|
|
3 |
|
|
|
|
|
|
|
|
|
if di |
= |
3, |
|
|
|
|
|
|
|
|
|
|
|
||||
α = |
16 |
|
|
|
|
|
|
|
|
|
(4.16) |
|||
1 |
|
[ |
5 |
− ( |
3 |
+ |
1 |
2π 2 |
] |
|
|
|||
|
|
|
|
8 |
8 |
4 cos |
di ) |
if di > |
3. |
|||||
|
|
di |
||||||||||||
3.The subdivided surface is given by connecting the new edge vertices and the updated original vertices as shown in Fig. 4.15(c).
After one subdivision, all vertices have degree six except those which were in the original mesh and had a degree other than six. These are the extraordinary vertices. The limit surface is C2 continuous except at the extraordinary vertices. Figure 4.16 shows loop subdivision applied to an icosahedron base mesh.
4.5 Local Differential Properties
Many useful 3D data processing operations require the computation of local differential properties of a surface. These range from the construction of local features such as spin images [26] to the computation of geodesic paths over a manifold [30].
4 Representing, Storing and Visualizing 3D Data |
167 |
However, most surface representations are discrete and contain only an approximate sampling of the underlying surface. The nature of the representation determines the manner in which these properties are computed.
First order properties of the surface (normal vector and tangent plane) are most commonly used for shading (since surface reflectance is a function of orientation). Interpolation shading uses interpolated surface normals to allow a polygonal approximation of a smooth surface to appear smooth when rendered (see Fig. 4.21(g)). Second order properties (principal curvatures and directions) are often used for local characterization of the surface topology, for example the shape index [32]. Occasionally, even third order properties (directional derivatives of the principal curvatures) can be useful.
For surfaces described in functional form (such as implicit surfaces), differential properties can be computed analytically using differential calculus. On the other hand, discrete representations require the assumption that the underlying surface is smooth and differential properties can only be approximated. Such operators should converge asymptotically to the true result as the sampling density increases.
4.5.1 Surface Normals
For uniformly sampled representations such as voxels and depth maps, differential properties can easily be approximated using finite differences which consider adjacent pixels or voxels, corresponding to the local neighborhood about a point. For example, in the simplest case, the surface gradients in the x and y-directions of a surface represented as a discrete depth map, z(x, y), can be approximated using single forward differences:
∂x z(x, y) ≈
where δx and δy then given by:
z(x + 1, y) − z(x, y) |
, |
∂ |
z(x, y) |
≈ |
z(x, y + 1) − z(x, y) , |
|
|||||
δx |
y |
|
δy |
||
|
|
|
|
|
(4.17) |
are the spacings on the pixel array. The surface normal vector is
n(x, y) |
= |
[∂x z(x, y)∂y z(x, y)1]T |
|
. |
(4.18) |
|
[∂x z(x, y)∂y z(x, y)1]T |
|
|||||
|
|
|
Note that computing gradients using only forward differences results in high sensitivity to noise. For this reason, a window about each pixel can be used in conjunction with an appropriate filter to provide a more stable estimate of the gradients [45].
In the case of an arbitrary triangular mesh, such a simple approach is not possible. If the surface is truly piecewise planar, then it is enough to associate a face normal with each triangular facet. For a face composed of vertices v1, v2 and v3, the face normal is given by:
n |
f = |
(v1 |
− v2) × (v1 |
− v3) |
. |
(4.19) |
|
|
|
||||
|
(v1 − v2) × (v1 − v3) |
|
||||
168 |
W.A.P. Smith |
Fig. 4.17 Computation of vertex normal by averaging adjacent face normals
For the normal to be directed to the outside of the object, the vertices must be ordered in an anticlockwise sense when viewed from the outside.
More commonly, the underlying surface is assumed smooth and the surface normal is approximated at each vertex. There are a number of alternative approaches to computing vertex normals [25], which are based on weighted averages of the adjacent facet normals (see Fig. 4.17).
Consider a vertex whose incident edges in anticlockwise order are given by the sequence e1, . . . , en . The first edge is repeated at the end of the sequence so e1 = en. The most well known method which is very efficient to compute uses triangle areas as weights. The area weighted normal to a facet defined by edges ei and ei+1 is given simply by ei × ei+1. Hence, the area weighted vertex normal is:
n−1 |
|
nvArea = ei × ei+1. |
(4.20) |
i=1 |
|
Note that this result requires normalization back to unit length. The efficiency of this approach lies in the fact that the cross product computation factors in the area weight at no extra computational cost. In particular, no trigonometric calculations are required. The downside to this approach is that a triangle with large area but small angle between the edges incident on the vertex contributes disproportionately to the result. This can lead to highly inaccurate normals under certain circumstances and the result depends heavily on the mesh triangulation (see [41] for an example).
An alternative, originally proposed by Thürmer and Wüthrich [70] uses the angle between pairs of edges incident on a vertex as the weight:
|
Angle |
n−1 |
|
ei |
· |
ei 1 |
ei |
× |
ei |
|
1 |
|
|
n |
|
arccos |
|
|
+ |
|
|
+ |
|
. |
(4.21) |
||
|
|
ei ei+1 ei |
× ei+1 |
||||||||||
|
v |
= i 1 |
|
|
|
||||||||
|
|
= |
|
|
|
|
|
|
|
|
|
|
|
Again, the result requires normalization. This approach is relatively simple and accurate but requires trigonometric calculations so is unsuitable for applications involving real-time computation of vertex normals (for example interactive rendering of a deforming surface).
A method proposed by Max [41] offers a compromise. It does not require trigonometric calculations, but is more stable in the special cases where area weighting performs badly. The weight is comprised of the sine of the angle between edge vectors
4 Representing, Storing and Visualizing 3D Data |
169 |
and the edge length reciprocal:
nMax |
|
n−1 |
sin αi |
|
|
ei × ei+1 |
|
|||||
|
|
|
|
, |
||||||||
v |
= i 1 |
ei ei+1 |
|
ei × ei+1 |
|
|||||||
|
|
= |
|
|
|
|
|
|
|
|
|
|
where |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
sin α |
|
ei × ei+1 |
. |
|
|
|||||
Which simplifies to: |
|
|
|
i = |
ei ei+1 |
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
nMax |
|
n−1 |
|
ei × ei+1 |
|
|
|||||
|
|
|
|
. |
|
|||||||
|
v |
|
= i 1 |
ei 2 |
ei+1 2 |
|
|
|||||
|
|
|
|
= |
|
|
|
|
|
|
|
|
(4.22)
(4.23)
(4.24)
Again, the result requires normalization. The derivation follows from an assumption of a locally spherical surface. For meshes representing surfaces which are exactly locally spherical, the result is exact.
4.5.2 Differential Coordinates and the Mesh Laplacian
When a surface is represented using a mesh, a transformation can be made to differential coordinates (δ-coordinates), which provides an alternate route to the computation of differential surface properties as well as a useful representation for other operations (including mesh compression). This is most easily accomplished when the mesh is stored in an adjacency matrix. A vertex vi , {i} KN which is a member of the mesh M N = (KN , S) may be represented in terms of δ-coordinates. This representation describes the difference between the absolute position of the vertex and the center of mass of its adjacent vertices:
x |
y |
z |
|
T |
|
|
1 |
|
|
|
δi = δi |
δi |
δi |
|
|
= vi − |
|
Adj(i) |
|
vj . |
(4.25) |
|
|
|
|
|
|
|
j Adj(i) |
|
||
di = Adj(i) is the degree of the vertex i. This transformation can be represented in matrix form. Given the mesh adjacency matrix A, where
|
|
if {i, j } |
KN |
|
Aij = |
1 |
|
(4.26) |
|
0 |
otherwise, |
|||
and the diagonal degree matrix, D, where Dii = di , the graph Laplacian matrix (sometimes called the topological Laplacian) is defined by L = D − A, i.e.
|
|
if i |
= |
j |
|
di |
|
|
|||
Lij = −1 |
if {i, j } KN |
(4.27) |
|||
|
0 |
otherwise. |
|
||
|
|
|
|
|
|