### Chapter6: Surfaces From Point Clouds

The objects scanned in the previous examples are solid, with a well-defined boundary surface separating the inside from the outside. Since computers have a finite amount of memory and operations need to be completed in a finite number of steps, algorithms can only be designed to manipulate surfaces described by a finite number of parameters. Perhaps the simplest surface representation with a finite number of parameters is produced by a finite sampling scheme, where a process systematically chooses a set of points lying on the surface.

The triangulation-based 3D scanners described in previous chapters
produce such a finite sampling scheme. The so-called * point
cloud*, a dense collection of surface samples, has become a
popular representation in computer graphics. However, since point
clouds do not constitute surfaces, they cannot be used to determine
which 3D points are inside or outside of the solid object. For many
applications, being able to make such a determination is critical. For
example, without closed bounded surfaces, volumes cannot be
measured. Therefore, it is important to construct so-called *
watertight* surfaces from point clouds. In this chapter we consider
these issues.

### Representation and Visualization of Point Clouds

In addition to the 3D point locations, the 3D scanning methods described in previous chapters are often able to estimate a color per point, as well as a surface normal vector. Some methods are able to measure both color and surface normal, and some are able to estimate other parameters which can be used to describe more complex material properties used to generate complex renderings. In all these cases the data structure used to represent a point cloud in memory is a simple array. A minimum of three values per point are needed to represent the point locations. Colors may require one to three more values per point, and normals vectors three additional values per point. Other properties may require more values, but in general it is the same number of parameters per point that need to be stored. If \(M\) is the number of parameters per point and \(N\) is the number of points, then point cloud can be represented in memory using an array of length \(NM\).

#### File Formats

Storing and retrieving arrays from files is relatively simple, and storing the raw data either in ASCII format or in binary format is a valid solution to the problem. However, these solutions may be incompatible with many software packages. We want to mention two standards which have support for storing point clouds with some auxiliary attributes.

**Storing Point Clouds as VRML Files**

The Virtual Reality Modeling Language (VRML) is an ISO standard
published in 1997. A VRML file describes a scene graph comprising a
variety of nodes. Among geometry nodes, `PointSet`

and `IndexedFaceSet`

are used to store point clouds. The
`PointSet`

node was designed to
store point clouds, but in addition to the 3D coordinates of each
point, only colors can be stored. No other attributes can be stored in
this node. In particular, normal vectors cannot be recorded. This is a
significant limitation, since normal vectors are important both for
rendering point clouds and for reconstructing watertight surfaces from
point clouds.

The `IndexedFaceSet`

node was
designed to store polygon meshes with colors, normals, and/or texture
coordinates. In addition to vertex coordinates, colors and normal
vectors can be stored bound to vertices. Even though
the `IndexedFaceSet`

node was not
designed to represent point clouds, the standard allows for this node
to have vertex coordinates and properties such as colors and/or
normals per vertex, but no faces. The standard does not specify how
such a node should be rendered in a VRML browser, but since they
constitute valid VRML files, they can be used to store point clouds.

**The SFL File Format**

The SFL file format was introduced with
Pointshop3D `[ZPKG02]`

to provide
a versatile file format to import and export point clouds with color,
normal vectors, and a radius per vertex describing the local sampling
density. A SFL file is encoded in binary and features an extensible
set of surfel attributes, data compression, upward and downward
compatibility, and transparent conversion of surfel attributes,
coordinate systems, and color spaces. Pointshop3D is a software
system for interactive editing of point-based surfaces, developed at
the Computer Graphics Lab at ETH Zurich.

#### Visualization

A well-established technique to render dense point clouds is point splatting. Each point is regarded as an oriented disk in 3D, with the orientation determined by the surface normal evaluated at each point, and the radius of the disk usually stored as an additional parameter per vertex. As a result, each point is rendered as an ellipse. The color is determined by the color stored with the point, the direction of the normal vector, and the illumination model. The radii are chosen so that the ellipses overlap, resulting in the perception of a continuous surface being rendered.

### Merging Point Clouds

The triangulation-based 3D scanning methods described in previous chapters are able to produce dense point clouds. However, due to visibility constraints these point clouds may have large gaps without samples. In order for a surface point to be reconstructed, it has to be illuminated by a projector, and visible by a camera. In addition, the projected patterns needs to illuminate the surface transversely for the camera to be able to capture a sufficient amount of reflected light. In particular, only points on the front-facing side of the object can be reconstructed (i.e., on the same side as the projector and camera). Some methods to overcome these limitations are discussed in the Conclusion. However, to produce a complete representation, multiple scans taken from various points of view must be integrated to produce a point cloud with sufficient sampling density over the whole visible surface of the object being scanned.

#### Computing Rigid Body Matching Transformations

The main challenge to merging multiple scans is that each scan is produced with respect to a different coordinate system. As a result, the rigid body transformation needed to register one scan with another must be estimated. In some cases the object is moved with respect to the scanner under computer control. In those cases the transformations needed to register the scans are known within a certain level of accuracy. This is the case when the object is placed on a computer-controlled turntable or linear translation stage. However, when the object is repositioned by hand, the matching transformations are not known and need to be estimated from point correspondences.

We now consider the problem of computing the rigid body transformation \(q=Rp+T\) to align two shapes from two sets of \(N\) points, \(\{p_1,\ldots,p_N\}\) and \(\{q_1,\ldots,q_N\}\). That is, we are looking for a rotation matrix \(R\) and a translation vector \(T\) so that

\[ q_1=Rp_1+T\quad\ldots\quad q_N=Rp_N+T\;. \]The two sets of points can be chosen interactively or automatically. In either case, being able to compute the matching transformation in closed form is a fundamental operation.

This registration problem is, in general, not solvable due to measurement errors. A common approach in such a case is to seek a least-squares solution. In this case, we desire a closed-form solution for minimizing the mean squared error

\begin{equation} \phi(R,T)=\frac{1}{N}\sum_{i=1}^N \|Rp_i+T-q_i\|^2\;, \label{eq:matching-error-RT} \end{equation}over all rotation matrices \(R\) and translation vectors \(T\). This yields a quadratic function of 12 components in \(R\) and \(T\); however, since \(R\) is restricted to be a valid rotation matrix, there exist additional constraints on \(R\). Since the variable \(T\) is unconstrained, a closed-form solution for \(T\), as a function of \(R\), can be found by solving the linear system of equations resulting from differentiating the previous expression with respect to \(T\).

\[ \frac{1}{2}\frac{\partial\phi}{\partial T} = \frac{1}{N}\sum_{i=1}^N (Rp_i+T-q_i) = 0 \Rightarrow T= \overline{q}-R\overline{p} \]In this expression \(\overline{p}\) and \(\overline{q}\) are the geometric centroids of the two sets of matching points, given by

\[ \overline{p}=\left(\frac{1}{N}\sum_{i=1}^Np_i\right)\qquad \overline{q}=\left(\frac{1}{N}\sum_{i=1}^Nq_i\right)\;. \]Substituting for \(T\) in Equation \ref{eq:matching-error-RT}, we obtain the following equivalent error function which depends only on \(R\).

\begin{equation} \psi(R)=\frac{1}{N}\sum_{i=1}^N \|R(p_i-\overline{p})-(q_i-\overline{q})\|^2 \label{eq:matching-error-R} \end{equation}If we expand this expression we obtain

\[ \psi(R)= \frac{1}{N}\sum_{i=1}^N \|p_i-\overline{p}\|^2 -\frac{2}{N}\sum_{i=1}^N (q_i-\overline{q})^tR(p_i-\overline{p}) +\frac{1}{N}\sum_{i=1}^N \|q_i-\overline{q}\|^2\;, \]since \(\|Rv\|^2=\|v\|^2\) for any vector \(v\). As the first and last terms do not depend on \(R\), maximizing this expression is equivalent to maximizing

\[ \eta(R)=\frac{1}{N}\sum_{i=1}^N (q_i-\overline{q})^tR(p_i-\overline{p}) = \hbox{trace}(RM)\;, \]where \(M\) is the \(3\times 3\) matrix

\[ M = \frac{1}{N}\sum_{i=1}^N (p_i-\overline{p})(q_i-\overline{q})^t\;. \]Recall that, for any pair of matrices \(A\) and \(B\) of the same dimensions, \(\hbox{trace}(A^tB)=\hbox{trace}(BA^t)\). We now consider the singular value decomposition (SVD) \(M=U\Delta V^t\), where \(U\) and \(V\) are orthogonal \(3\times 3\) matrices, and \(\Delta\) is a diagonal \(3\times 3\) matrix with elements \(\delta_1\geq\delta_2\geq\delta_3\geq 0\). Substituting, we find

\[ \hbox{trace}(RM)= \hbox{trace}(RU\Delta V^t)= \hbox{trace}((V^tRU)\Delta)= \hbox{trace}(W\Delta)\;, \]where \(W=V^tRU\) is orthogonal. If we expand this expression, we obtain

\[ \hbox{trace}(W\Delta)=w_{11}\delta_1+w_{22}\delta_2+w_{33}\delta_3 \leq\delta_1+\delta_2+\delta_3\;, \]where \(W=(w_{ij})\). The last inequality is true because the components of an orthogonal matrix cannot be larger than one. Note that the last inequality is an equality only if \(w_{11}=w_{22}=w_{33}=1\), which is only the case when \(W=I\) (the identity matrix). It follows that if \(V^tU\) is a rotation matrix, then \(R=V^tU\) is the minimizer of our original problem. The matrix \(V^tU\) is an orthogonal matrix, but it may not have a negative determinant. In that case, an upper bound for \(\hbox{trace}(W\Delta)\), with \(W\) restricted to have a negative determinant, is achieved for \(W=J\), where

\[ J=\left( \begin{matrix} 1 & 0 & 0 \cr 0 & 1 & 0 \cr 0 & 0 & -1 \end{matrix} \right)\;. \]In this case it follows that the solution to our problem is \(R=V^tJU\).

#### The Iterative Closest Point (ICP) Algorithm

The Iterative Closest Point (ICP) is an algorithm employed to match two surface representations, such as points clouds or polygon meshes. This matching algorithm is used to reconstruct 3D surfaces by registering and merging multiple scans. The algorithm is straightforward and can be implemented in real-time. ICP iteratively estimates the transformation (i.e., translation and rotation) between two geometric data sets. The algorithm takes as input two data sets, an initial estimate for the transformation, and an additional criterion for stopping the iterations. The output is an improved estimate of the matching transformation. The algorithm comprises the following steps:

- Select points from the first shape.
- Associate points, by nearest neighbor, with those in the second shape.
- Estimate the closed-form matching transformation using the method derived in the previous section.
- Transform the points using the estimated parameters.
- Repeat previous steps until the stopping criterion is met.

The algorithm can be generalized to solve the problem of registering multiple scans. Each scan has an associated rigid body transformation which will register it with respect to the rest of the scans, regarded as a single rigid object. An additional external loop must be added to the previous steps to pick one transformation to be optimized with each pass, while the others are kept constant---either going through each of the scans in sequence, or randomizing the choice.

### Surface Reconstruction from Point Clouds

Watertight surfaces partition space into two disconnected regions so that every line segment joining a point in one region to a point in the other must cross the dividing surface. In this section we discuss methods to reconstruct watertight surfaces from point clouds.

#### Continuous Surfaces

In mathematics surfaces are represented in parametric or implicit
form. A parametric surface \(S=\{x(u):u\in U\}\) is defined by a
function \(x:U\rightarrow\mathbb{R}^3\) on an open subset \(U\) of the plane. An
implicit surface is defined as a level set
\(S=\{p\in\mathbb{R}^3:f(p)=\lambda\}\) of a continuous function
\(f:V\rightarrow\mathbb{R}\), where \(V\) is an open subset in 3D. These
functions are most often smooth or piecewise smooth. Implicit
surfaces are called *watertight* because they partition space into
the two disconnected sets of points, one where \(f(p)>\lambda\) and a
second where \(f(p)<\lambda\). Since the function \(f\) is continuous,
every line segment joining a point in one region to a point in the
other must cross the dividing surface. When the boundary surface of a
solid object is described by an implicit equation, one of these two
sets describes the inside of the object, and the other one the
outside. Since the implicit function can be evaluated at any point in
3D space, it is also referred to as a scalar field. On the other hand,
parametric surfaces may or may not be watertight. In general, it is
difficult to determine whether a parametric surface is watertight or
not. In addition, implicit surfaces are preferred in many
applications, such as reverse engineering and interactive shape
design, because they bound a solid object which can be manufactured;
for example, using rapid prototyping technologies or
numerically-controlled machine tools, such representations can define
objects of arbitrary topology. As a result, we focus our remaining
discussion on implicit surfaces.

#### Discrete Surfaces

A discrete surface is defined by a finite number of parameters. We
only consider here polygon meshes, and in particular those polygon
meshes representable as `IndexedFaceSet`

nodes in VRML files. Polygon meshes are
composed of *geometry* and topological
*connectivity*. The geometry includes vertex coordinates,
normal vectors, and colors (and possibly texture coordinates). The
connectivity is represented in various ways. A popular representation
used in many isosurface algorithms is the *polygon soup*, where
polygon faces are represented as loops of vertex coordinate
vectors. If two or more faces share a vertex, the vertex coordinates
are repeated as many times as needed. Another popular representation
used in isosurface algorithms is the \texttt{IndexedFaceSet} (IFS),
describing polygon meshes with simply-connected faces. In this
representation the geometry is stored as arrays of floating point
numbers. In these notes we are primarily concerned with the array
\texttt{coord} of vertex coordinates, and to a lesser degree with the
array \texttt{normal} of face normals. The connectivity is described
by the total number \(V\) of vertices, and \(F\) faces, which are
stored in the \texttt{coordIndex} array as a sequence of loops of
vertex indices, demarcated by values of \(-1\).

#### Isosurfaces

An isosurface is a polygonal mesh surface representation produced by
an isosurface algorithm. An isosurface algorithm constructs a
polygonal mesh approximation of a smooth implicit surface
\(S=\{x:f(x)= 0\}\) within a bounded three-dimensional volume, from
samples of a defining function \(f(x)\) evaluated on the vertices of a
volumetric grid. Marching Cubes
```
[LC87]
```

and related algorithms
operate on function values provided at the vertices of hexahedral
grids. Another family of isosurface algorithms operate on functions
evaluated at the vertices of tetrahedral grids
```
[DK91]
```

. Usually, no additional information about the function is
provided, and various interpolation schemes are used to evaluate the
function within grid cells, if necessary. The most natural
interpolation scheme for tetrahedral meshes is linear interpolation,
which we also adopt here.

#### Isosurface Construction Algorithms

An isosurface algorithm producing a polygon soup output must solve three key problems: (1) determining the quantity and location of isosurface vertices within each cell, (2) determining how these vertices are connected forming isosurface faces, and (3) determining globally consistent face orientations. For isosurface algorithms producing IFS output, there is a fourth problem to solve: identifying isosurface vertices lying on vertices and edges of the volumetric grid. For many visualization applications, the polygon soup representation is sufficient and acceptable, despite the storage overhead. Isosurface vertices lying on vertices and edges of the volumetric grid are independently generated multiple times. The main advantage of this approach is that it is highly parallelizable. But, since most of these boundary vertices are represented at least twice, it is not a compact representation.

Researchers have proposed various solutions and design decisions (e.g., cell types, adaptive grids, topological complexity, interpolant order) to address these four problems. The well-known Marching Cubes (MC) algorithm uses a fixed hexahedral grid (i.e., cube cells) with linear interpolation to find zero-crossings along the edges of the grid. These are the vertices of the isosurface mesh. Second, polygonal faces are added connecting these vertices using a table. The crucial observation made with MC is that the possible connectivity of triangles in a cell can be computed independently of the function samples and stored in a table. Out-of-core extensions, where sequential layers of the volume are processed one at a time, are straightforward.

**Figure 6.1**

*In isosurface algorithms, the sign of the function at the grid vertices determines the topology and connectivity of the output polygonal mesh within each tetrahedron. Mesh vertices are located on grid edges where the function changes sign.}*

Similar tetrahedral-based
algorithms
`[DK91][GH95][TPG99]`

,
dubbed Marching Tetrahedra (MT), have also been developed (again using
linear interpolation). Although the cell is simpler, MT requires
maintaining a tetrahedral sampling grid. Out-of-core extensions
require presorted traversal schemes, such as in
`[CS97]`

.
For an unsorted
tetrahedral grid, hash tables are used to save and retrieve vertices
lying on edges of the volumetric grid. As an example of an isosurface
algorithm, we discuss MT in more detail.

0 | 0000 | [-1] |

1 | 0001 | [2,1,0,-1] |

2 | 0010 | [0,3,4,-1] |

3 | 0011 | [1,3,4,2,-1] |

4 | 0100 | [1,5,3,-1] |

5 | 0101 | [0,2,5,3,-1] |

6 | 0110 | [0,3,5,4,-1] |

7 | 0111 | [1,5,2,-1] |

8 | 1000 | [2,5,1,-1] |

9 | 1001 | [4,5,3,0,-1] |

10 | 1010 | [3,5,2,0,-1] |

11 | 1011 | [3,5,1,-1] |

12 | 1100 | [2,4,3,1,-1] |

13 | 1101 | [4,3,0,-1] |

14 | 1110 | [0,1,2,-1] |

15 | 1111 | [-1] |

0 | (0,1) |

1 | (0,2) |

2 | (0,3) |

3 | (1,2) |

4 | (1,3) |

5 | (2,3) |

**Table 6.1**

*Look-up tables for tetrahedral mesh isosurface evaluation. Note that consistent face orientation is encoded within the table. Indices stored in the first table reference tetrahedron edges, as indicated by the second table of vertex pairs (and further illustrated in Figure 6.1 ). In this case, only edge indices \(\{1,2,3,4\}\) have associated isosurface vertex coordinates, which are shared with neighboring cells.*

**Marching Tetrahedra**

MT operates on the following input data: (1) a tetrahedral grid and (2) one piecewise linear function \(f(x)\), defined by its values at the grid vertices. Within the tetrahedron with vertices \(x_0,x_1,x_2,x_3\in\mathbb{R}^3\), the function is linear and can be described in terms of the barycentric coordinates \(b=(b_0,b_1,b_2,b_3)^t\) of an arbitrary internal point \(x = b_0\,x_0+b_1\,x_1+b_2\,x_2+b_3\,x_3\) with respect to the four vertices: \(f(x)=b_0\,f(x_0)+b_1\,f(x_1)+b_2\,f(x_2)+b_3\,f(x_3)\), where \(b_0,b_1,b_2,b_3\geq 0\) and \(b_0+b_1+b_2+b_3=1\). As illustrated in Figure 6.1 , the sign of the function at the four grid vertices determines the connectivity (e.g., triangle, quadrilateral, or empty) of the output polygonal mesh within each tetrahedron. There are actually \(16=2^4\) cases, which modulo symmetries and sign changes reduce to only three. Each grid edge, whose end vertex values change sign, corresponds to an isosurface mesh vertex. The exact location of the vertex along the edge is determined by linear interpolation from the actual function values, but note that the \(16\) cases can be precomputed and stored in a table indexed by a 4-bit integer \(i=(i_3i_2i_1i_0)\), where \(i_j=1\) if \(f(x_j)>0\) and \(i_j=0\), if \(f(x_j)<0\). The full table is shown in Table 6.1. The cases \(f(x_j)=0\) are singular and require special treatment. For example, the index is \(i=(0100)=4\) for Figure 6.1 (a), and \(i=(1100)=12\) for Figure 6.1 (b). Orientation for the isosurface faces, consistent with the orientation of the containing tetrahedron, can be obtained from connectivity alone (and are encoded in the look-up table as shown in Table 6.1). For IFS output it is also necessary to stitch vertices as described above.

Algorithms to polygonize implicit surfaces
`[Blo88]`

,
where the implicit functions are provided in analytic form, are
closely related to isosurface algorithms. For example, Bloomenthal and
Ferguson
`[BF95]`

extract
non-manifold isosurfaces produced from trimming implicits and
parameterics using a tetrahedral isosurface algorithm.
`[WvO96]`

polygonize
boolean combinations of skeletal implicits (Boolean Compound Soft
Objects), applying an iterative solver and face subdivision for
placing vertices along feature edges and points. Suffern and
Balsys~`[SB03]`

present an
algorithm to polygonize surfaces defined by two implicit functions
provided in analytic form; this same algorithm can compute
bi-iso-lines of pairs of implicits for rendering.

**Isosurface Algorithms on Hexahedral Grids**

An isosurface algorithm constructs a polygon mesh approximation of a level set of a scalar function defined in a finite 3D volume. The function \(f(p)\) is usually specified by its values \(f_{\alpha}=f(p_\alpha)\) on a regular grid of three dimensional points

\[ G=\{p_\alpha:\alpha=(\alpha_0,\alpha_1,\alpha_2)\in [ n_0]\!\times\! [ n_1 ]\!\times\! [ n_2]\} \;, \]
where \([ n_j]=\{0,\ldots,n_j-1\}\), and by a method to
interpolate in between these values. The surface is usually
represented as a polygon mesh, and is specified by its * isovalue*
\(f_0\). Furthermore, the interpolation scheme is assumed to be linear
along the edges of the grid, so that the isosurface cuts each edge in
no more than one point. If \(p_\alpha\) and \(p_\beta\) are grid points
connected by an edge, and \(f_\alpha>f_0>f_\beta\), the location of the
point \(p_{\alpha\beta}\) where the isosurface intersects the edge is

**Marching Cubes**

One of the most popular isosurface extraction algorithms is *
Marching Cubes* ```
[LC87]
```

. In this algorithm the points
defined by the intersection of the isosurface with the edges of the
grid are the vertices of the polygon mesh. These vertices are
connected forming polygon faces according to the following
procedure. Each set of eight neighboring grid points define a small
cube called a * cell*

Since the function value associated with each of the eight corners of a cell may be either above or below the isovalue (isovalues equal to grid function values are called singular and should be avoided), there are \(2^8=256\) possible configurations. A polygonization of the vertices within each cell, for each one of these configurations, is stored in a static look-up table. When symmetries are taken into account, the size of the table can be reduced significantly.

#### Algorithms to Fit Implicit Surfaces to Point Clouds

Let \(U\) be a relatively open and simply-connected subset of
\(\mathbb{R}^3\), and \(f:U\rightarrow \mathbb{R}\) a smooth function.
The gradient \(\,\nabla f\,\) is a vector field defined on \(U\).
Given an * oriented point cloud*, i.e., a finite set \({\cal
D}\) of point-vector pairs \((p,n)\), where \(p\) is an interior point
of \(U\), and \(n\) is a unit length 3D vector, the problem is to find
a smooth function \(f\) so that \(f(p)\approx 0\) and \(\nabla
f(p)\approx n\) for every oriented point \((p,n)\) in the data set
\({\cal D}\). We call the zero iso-level set of such a function \(\{ p
: f(p)=0 \}\) a * surface fit*, or * surface
reconstruction*, for the data set \({\cal D}\).

We are particularly interested in fitting isosurfaces to oriented point points. For the sake of simplicity, we assume that the domain is the unit cube \(U=[0,1]^3\), the typical domain of an isosurface defined on an hexahedral mesh, and the isolevel is zero, i.e., the isosurface to be fitted to the data points is \(\{p:f(p)=0\}\), but of course, the argument applies in more general cases.

**Figure 6.2**

*Early results of Vector Field Isosurface reconstruction from oriented point clouds introduced in*

`[ST05]`

.
Figure 6.2 shows results of surface reconstruction from an
oriented point cloud using the simple variational formulation
presented in `[ST05]`

, where
oriented data points are regarded as samples of the gradient vector
field of an implicit function, which is estimated by minimizing this
energy function

where \(f(x)\) is the implicit function being estimated, \(\nabla f(x)\)
is the gradient of \(f\), \(Hf(x)\) is the Hessian of \(f(x)\),
\((p_1,n_1),\ldots,(p_m,n_m)\) are point-normal data pairs, \(V\) is a
bounding volume, and \(\lambda_1\) and \(\lambda_2\) are regularization
parameters. Minimizing this energy requires the solution of a simple
large and sparse least squares problem. The result is usually unique
modulo an additive constant. Given that, for rendering or
post-processing, isosurfaces are extracted from scalar functions
defined over regular grids (e.g., via Marching Cubes), it is worth
exploring representations of implicit functions defined as a regular
scalar fields. Finite difference discretization is used in
`[ST05]`

, with the
volume integral resulting in a sum of gradient differences over the
edges of the regular grid, yet Equation \ref{eq:lsf1} can be
discretized in many other ways.