Keywords: Point cloud, Morphing, Clustering
The representation of geometric models by large sets of point samples (a.k.a point clouds) constitute one of the canonical data formats for scientific data visualization. We can acquire point clouds from the measurement of some physical process. Point clouds can represent surfaces, volumetric or isosurface data. The availability of the modern 3D scanners brings the possibility to acquire point sets representing the surface of the analyzed solid, that contain millions of sample points.
Point cloud can be for some applications better representation than widely used boundary representation. This holds mainly for very complex models, such as fractal surfaces. It is generally good idea to store these object as point clouds, because the algorithms for conversion to the surface representation (such as polygonal mesh) are very computationally involved and require great amounts of main memory. One reason for the inefficiency of boundary representation is that highly detailed models contain a large number of small primitives, which fill lesser area than a pixel when they are projected and displayed.
Point cloud is the unstructured set of point samples. One point sample is elementary object, specified by its location in 3D space, normal vector, color, transparency and size. Single point sample can be visualized as a small sphere or a point (pixel). As presented in [2], [4] and [5], a point set representing the surface of a model can be rendered as a solid textured surface.
The morphing between two solids is the animation, during which the solid smoothly changes its shape from one shape to another. Our goal is implement and analyze several methods, which may be theoretically applicable to morphing between two point clouds. The methods should be independent of the topology of the models (should consider only the locations, possibly normals of the point samples) and should be able to perform morphing between two differentsized point clouds.
The fundamental principle underlying all introduced methods is the finding of the mapping
The assignment process can be based on various criteria. Probably the most commonly used criterion will be that the total distance, which the individual points travel during the morphing, is minimal. In this case the finding of mapping M is equivalent to the solving of the optimization problem
The trivial algorithm, which finds all possible assignments and chooses the optimal solution based on the optimization criterion, is unusable because of its exponential computational complexity. Where this algorithm can be used for finding the assignment between small point sets (max. 10 points), the models represented by point cloud contain thousands, possibly millions of sample points. We can find the assignment by incorporating the algorithms of artificial intelligence (searching in the state space), genetic algorithms, or space partitioning (clustering). This report introduces several methods based on clustering.
The idea of clustering is to group the points into a smaller number of clusters and find the mapping M between the two sets of clusters. The motivation behind the use of clustering is basically to reduce the problem size, solve the morphing on the higher level (cluster level) and then descend to the lower level. This approach is based on “divide and conquer” programming technique.
Clustering is one of generally used methods for the simplification of the pointsampled geometry. In the process of clustering the unstructured point cloud is divided to several smaller, spatially compact subsets (clusters). The clusters can be in the next step of simplification replaced by single point samples, whose size will reflect the size of replaced cluster. Two methods of clustering are clustering by incremental regiongrowing and hierarchical clustering. In our implementation we will use the latter method. This method is suitable for our needs, because it organizes the point cloud by creating the hierarchy of the subsets of point cloud. All implemented methods are based on the methods for efficient point cloud simplification from [1].
The hierarchical clustering is based on the partitioning of the point cloud. The recursive algorithm divides the point cloud P to several smaller point clouds (clusters). The result of the algorithm is tree, whose nodes represent point sets. The root represents the point cloud P, the leaves represent the terminal clusters (which are further indivisible and whose size is smaller than the specified limit or which contain only one point). The presented algorithms will be shown on the hierarchical clustering using binary space division (BSP tree). Other similar data structures, such as octrees, can be also used [4].
The hierarchical clustering using binary partition performs recursive division of point cloud to two parts by the split plane. The split plane is determined in most cases by the anchor point (usually the center point  centroid of the point set) and the normal vector. The choice of normal of the split plane has great impact on the quality of the morphing. The division algorithm is then applied to two resulting point sets. The cluster is divided only if its size is larger than the specified limit. In the text below, several possible methods for the finding of the normal of splitting plane are shown. As the anchor point of the splitting plane we will choose in most cases the centroid of the point set. For the covariance analysis is the centroid the only plausible choice; for the other methods we could use also other points, such as the center of bounding box/bounding sphere of the point set.
In this straightforward algorithm the normal of split plane is chosen alternately in the direction of x, y and z axis, depending on the depth of divided cluster in the binary tree (see figures 1 and 2).
This method is very easy to implement. The disadvantage of this method is that the point cloud is split by orthogonal cuts, which can be visible during the morphing. For the reduction of this undesirable phenomenon we can apply a random noise function to the normal vector of the split plane.
This method is based on the division of the (axesaligned) bounding box of the point cloud. The parameters of the bounding box are computed from the coordinates of all point samples and the bounding box is split into two equal halves by the plane perpendicular to its longest axis. This method does not necessarily involve the computation of the centroid for each cluster (the anchor point for the splitting plane can be chosen as the center of the bounding box), but we must find the dimensions and location of the bounding box. The bounding box can be specified by two points in 3D space, representing the lower front left and upper back right corners. This method leads to decent results, compared to orthogonal clustering. However, the method is more computationally expensive than the former method, since the finding of the corners of the bounding box takes several times (roughly twice) more time than the finding of the centroid.
The normal vector of the split plane can be found using the covariance analysis of the neighborhood of the centroid of the point cloud. The covariance analysis allows us to estimate various local surface properties, such as normal vector of approximation surface or surface variation. The covariance matrix C of the point cluster P is defined as
If we assume that _{0} < _{1} < _{2}, the plane
The normal vector of the split plane will be defined by the centroid of P and the largest eigenvector of the covariance matrix of P (which is v_{2}). The point cloud is always split along the direction of greatest variation, see the figure 3.
The drawback of this method is increased computational complexity (although the asymptotic complexity remains the same). Besides the centroid we must evaluate the coefficients of covariance matrix C (since the matrix is symmetric, there are six coefficients), find the roots of characteristic polynomial of C and finally compute the eigenvector corresponding to the largest eigenvalue.
The trees T_{s} and T_{d} created by the clustering of starting and final pointcloud model in the next stage serve as input for the algorithm for the generation of assignments M = {(s_{i}, d_{i})}. The algorithm is pretty straightforward: in each step will be assigned one cluster from T_{s} to some cluster from T_{d}. In the beginning the roots of both trees are assigned one to another and the assignment is put to the ordinary queue. In the next step, the assignment (c_{s } , c_{d}) is extracted from the queue. If any of the clusters c_{s}, c_{d } is a leaf, the assignment is put again to the queue, else the children of c_{s} and c_{d} are examined, assigned to each other and these assignments are put to the queue. This step is repeated until the queue contains only the assignments of type (c_{s}, c_{d}), where at least one of c_{s} or c_{d } represents a leaf in the corresponding tree (terminal cluster). For an example of possible assignment see the figure 4.
The way clusters are assigned to each other depends on the chosen metric function. Since we have defined the metric only between two points (which will be in our case an Euclidean distance between two points, (s, d) = s  d), we have to extend the definition, so it could be applicable also to clusters. The metric (a, b) between two clusters a T_{s}, b T_{d} is defined as (p_{a}, p_{b}), where p_{a} and p_{b } are centroids of clusters a and b respectively. Let us assume that both trees are regular with the same degree n (every internal node has exactly n children). During the finding of assignments of children of source node c_{s } to the children of destination node c_{d} we seek such permutation of n elements, which minimizes the following expression:
The algorithm presented above generates the set of assignments M' = {(c_{s}, c_{d})}, where c_{s} and c_{d} are clusters from source and destination point clouds respectively. We need to generate the assignment of individual point samples. This can be done by replacing the assignments between clusters by the sets of assignment between points. The assignment (A, B) M', where A = (a_{1}, a_{2}...a_{k}) and B = (b_{1}, b_{2}...b_{l}), will be replaced by the set of m assignments (a_{i}, b_{j}), where m = max(k, l). Assuming that k < l, this set can be described as mapping of the set B onto the set A. The actual mapping we can choose deliberately or we can apply some criteria based on metric function. The only important feature of this mapping is that it should be uniform in the sense that the number of points from B assigned to individual points a_{i} A should be for all points a_{i} approximately the same. The situation can be simplified if we set for the clustering algorithm the maximum limit for the cluster size to 1, so the assignments between clusters will be always of type 1 : N or N : 1.
By computing the set M of the assignments between individual point samples we have reduced whole rather complicated problem with point cloud morphing to the trivial problem of morphing between the assigned point pairs. Considering the straight paths between assigned points (s_{i}, d_{i}), we can get the coordinates of the points of transition point cloud p_{i} using the linear interpolation
Another possible solution is to use the polar coordinates for normals, so instead of coordinates we have to store the angles (since all normals are unit vectors, we do not have to store the length of vector). The linear interpolation of angles will yield correct results. During the rendering the angles have to be transformed back to the normals. This method will give the same results as the rotation of vectors and is faster and less memory expensive (we have to store two angles instead of three components of normal).
For the purposes of testing we have developed an application that can load two point cloud models, perform the necessary data preprocessing, perform the hierarchical clustering, generate the assignments between point pairs and visualize the actual morphing process. Our application was written in C++ using OpenGL and Qt 2.3.0 and compiled by MSVC 6.0 under Windows. Some data structures and the format of input files were adopted from the existing application QSplat [6] (a visualization system for point cloud models).
The implemented methods for morphing between models represented by a point cloud have been tested on twenty pairs of models. For some of the pairs the computation times of the process of finding of the assignment of points have been measured. The asymptotic complexity of the orthogonal clustering, the clustering based on the covariance analysis and the computation of the assignment of the points have been calculated.
The orthogonal clustering and the axesaligned bounding box division have the O(n ^{.} log _{2}n) complexity. The clustering based on the covariance analysis have also the O(n ^{.} log _{2}n) complexity, they differ only in constant. The computation of the assignment of the points has the O(n) complexity.
The measured durations for the tested pairs of models are given in the table 1. Three pairs of models were tested: the sphere and the Stanford bunny (both models have the same number of points), the African mask and the Christmas star, and finally the lion and the dinosaur. The measured times correspond to the fact that the asymptotic complexities of all tested clustering methods differ only in constant.
We have decreased the computational complexity and the total time spent by finding the assignments by allowing to store the once found assignments of points to file and then load it again.
The progress of morphing of the tested pairs of models using various clustering methods is shown on the figures 610 (morphing of the sphere to the bunny), figures 1115 (morphing of the African mask to the Christmas star) and the figures 1620 (morphing of the lion to the dinosaur). Besides four tested clustering methods (orthogonal clustering, orthogonal clustering with random noise, bounding box division, covariance analysis) there is also shown for the sake of comparison the case when the assignment of points is chosen randomly.

In this section we will discuss the problems that provides the implemented methods. The main problem is the occasional occurrence of cracks in the model during the morphing. The cracks occur because of the space partitioning and the nonconvexity of the shape. Our opinion is that this problem can not be eliminated only on the assumption of the sum of squared distances as the criterion of quality. This criterion is usable only for convex shapes. In a convex shape every surface point can be connected with every other surface points with a line without intersecting the surface of the shape. If a convex shape is divided by a split plane then the results are always two convex shapes. On the other hand a nonconvex shape can be divided by a split plane to any number of convex or nonconvex shapes and this is the main reason for the occurrence of cracks, see figure 5.
Alexa et al. have presented in [3] methods for transformation of a nonconvex Brep shape to a convex shape. The 3D manifold objects are homeomorphic with a sphere, so the unit sphere is the convex shape on which the nonconvex shapes should be transformed. Those methods can be adapted for the point cloud representation.
One class of shapes are the star shapes. For a star shape exist at least one point O inside the shape which can be connected with every boundary point of the shape with a line without intersecting the surface of the shape. The point O is called origin. If the coordinate system is transformed to be the origin O in its center and the coordinates of each boundary point normalized then the star shape has been transformed to the unit sphere. The main problem is to find whether a shape is a star shape or not and find the origin O, especially for point cloud representation where no information about the surface is included.
Second class of shapes are nonconvex nonstar shapes. This class of shapes can be transformed to unit sphere using the barycentric mapping. The barycentric mapping uses a simple idea that every point is placed in centroid of its neighbors. The shape should be downsampled until it is a convex shape. The coordinate system should be transformed to be the centroid of the shape in its center, then the coordinates of the points should be normalized, and the downsampled points should be upsampled on unit sphere respecting the idea of barycentric mapping.
In this case the paths can not be lines. The paths should depend on function which transforms the nonconvex shape to unit sphere.
These problems are rather nontrivial ones and the future work will be focused on these particular problems.
For convex or just a little bit nonconvex models gives the best results the morphing based on the orthogonal clustering. The generally best method for nonconvex models can not be chosen. The testing has proved that the results of all methods depend on the shapes of the source model and the target model.
These methods are not usable for highquality morphing, so the methods have to be improved. The sum of squared distances can not be used as only criterion of quality for nonconvex shapes. Nonconvex shapes have to be transformed to convex shapes.
[1] Pauly, M.; Gross, M.; Kobbelt L. P.: Efficient Simplification of PointSampled Surfaces, IEEE Visualization, 2002
[2] Alexa, M.; Behr, J.; CohenOr, D.; Levin, D.; Silva, C. T.: Computing and Rendering Point Set Surfaces, IEEE Transactions on Computer Graphic and Visualization, 2002
[3] Alexa, M.: Mesh Morphing STAR, Eurographics 2001 State of the Art Reports, Eurographics, 2001
[4] Pfister, H.; van Baar, J.: Surfels: Surface elements as rendering primitives, In Proceedings of SIGGRAPH 2000, pages 335342, SIGGRAPH, 2000
[5] Zwicker, M.; Pfister, H.; van Baar, J.; Gross, M.: Surface Splatting, SIGGRAPH, 2001
[6] Rusinkiewicz, S.; Levoy, M.: QSplat: A Multiresolution Point Rendering System for Large Meshes, In Proceedings of SIGGRAPH 2000, SIGGRAPH, 2000














