Computation Of Normal Vectors Of Discrete 3d Objects: Application To Natural Snow Images From X-Ray Tomography

Estimating the normal vector field on the boundary of discrete 3D objects is essential for rendering and image measurement problems. Most of the existing algorithms do not provide an accurate determination of the normal vector field for shapes that present edges. We propose here a new and simple computational method to obtain accurate results on all types of shapes whatever their degree of local convexity. The presented method is based on the analysis of the gradient vector field of the distance map of the object. Some results on simulated data and snow images from X-ray tomography are presented and discussed.


INTRODUCTION
With the development of modern imaging techniques like microtomography or MRI, scientists are faced with the problem of the interpretation of 3D numerical images.
The calculation of the normal vector field on the boundary of discrete 3D objects is essential for a good visualization.Moreover, the knowledge of the normal vector allows the determination of many interesting parameters to analyze the object.
For example, algorithms for the determination of local curvature (Brzoska et al., 1999a) and surface area (Brzoska et al., 2001) were recently developed at the Centre d'Etudes de la Neige.These two algorithms offer interesting outlooks for modeling snow metamorphism.Nevertheless, they require an accurate computation of the normal vector field on the surface boundary to be relevant.
The normal vector field can commonly be determined by triangulation (Lorensen and Cline, 1987).This method is complex regarding to the accuracy of the obtained vector field.Therefore many imaging scientists developed fully discrete algorithms.These techniques provide good results on smooth surfaces such as sphere, tore and infinite plane surface (Lenoir et al., 1996;Thürmer and Würthrich, 1997).However, they often do not provide an accurate determination of the normal vector field for shapes that present edges (Papier and Françon, 1998).An interesting method (Tellier and Debled-Renesson, 1999) based on the tangential lines computation allows to obtain a realistic rendering on shapes that present edges.However, it does not seem very accurate when applied on small spheres.
We propose here a new and simple computational method to obtain accurate results on all types of shapes whatever their degree of local convexity.A discrete background distance map of the original 3D image is first constructed.For each voxel of the object, an elementary gradient vector of the distance map is then computed.The normal vector of a surface voxel is obtained by summation of each elementary gradient vector in an appropriate neighborhood.This neighborhood is determined for each surface voxel so that it does not contain any singularities of the gradient field.Some results obtained on real and simulated data are shown and commented.

METHOD
Normal vectors are computed at the center of each voxel V P of the surface S of an object O.
-V P is a voxel belonging to the background ⇔ V P ∉O.
-V P ∈S ⇔ V P has a 6-connected neighbor belonging to the background.
Our main idea for the computation of the normal vectors is to use the gradient vector field of the background distance map of the object.

DISTANCE MAP
Let P and Q be the centers of the voxels V P and V Q and d(P, Q) the distance from P to Q.
In a numerical object O, the background distance d b on a point P is defined as follows: (1) For instance, the background distance map can be built using a two-run iterative chamfer algorithm (Borgefors, 1984;Chassery and Montanvert, 1991).In summary, this algorithm increments at each run the discrete distance to the closest already processed voxel of the object.Depending on the size of the used chamfer kernel, some distortion with respect to the Euclidean distance map (Danielsson, 1980;Yamada, 1984;Ragnemalm, 1989) may occur in deeper parts of the map; however, it is negligible near the surface (Rolland, 1991), i.e. on the first layers of voxels.
For our method, we used a second-rank chamfer distance d 5-7-9-11 (d 5 for short), which seems a good compromise between the fastness of the first-rank chamfer distance d 3-4-5 and the accuracy of the Euclidean distance (Verwer, 1991;Borgefors, 2001).

GRADIENT MAP
Let N(P) be the 3D first-rank neighborhood of P and M i a point of N(P).
Let be Mj i ∈ N(P) / j (M i ) = j (P) -1 and M'j i ∈ N(P) / j (M i ) = j (P) + 1, where j = x, y, z are canonical coordinates.For each current voxel V P of O, the components of the gradient vector in P are computed as follows: The gradient map obtained is very sensitive to each detail of the distance map, and to its digitization effects, too.Thus, to compute an accurate normal vector field on the surface S of O, we have to sum these elementary gradient vectors on an appropriate neighborhood.

SUMMATION
A summation of the elementary gradient vectors in a fixed spherical neighborhood can induce significant errors for shapes that present singularities (like the edges or the vertices of cubes).
To overcome this problem, for each voxel V P of the surface S, the construction of the neighborhood is driven by the detection of the singularities of the gradient field.
The problem of the edge detection on a surface is replaced by the search of irregularities inside the gradient map (Fig. 1).
Let us define a critical angle a 0 so that, for an angle a < a 0 , the surface is assumed as having an edge.Given  Let r be the radius of the spherical neighborhood around P where local gradient contributions must be tested.For each point P ∈ S, the summation is calculated as follows: grad (P) is assumed to be the first estimation of the normal vector V(P) While (r ≤ r max ) { V(P) is updated by summing each contribution of the gradient map on the spherical cap of r radius.To be taken into account, each contribution from the current point Q should fulfill: Note that a 0 must be acute enough to obtain a smooth vector field and wide enough to detect discontinuities.In our case, we choose a 0 = 120°.This angle verifies the above conditions and corresponds to the angle of hexagonal ice crystals.

RESULTS
Our algorithm was tested on simulated data and natural snow images processed by X-ray microtomography (Brzoska et al., 1999b).The evaluation of the normal vectors was realized in two ways: -by comparing the obtained normal vectors to the theoretical normal vectors for synthetic objects whose surface is differentiable.
-by visualizing them for others objects.For visualization, u being the illumination vector, the gray level in P was simply computed as follows: Gray_level = -V(P) .u Except in the last test where the effect of r max is analyzed, r max was set to 5 voxels.

SPHERES
Tests were carried out on a set of digitized spheres from 1 to 60 voxels of radius.Spheres were obtained according to the definition proposed by Papier and Françon (1998) and Tellier and Debled-Renesson (1999): A sphere with an integer radius r is defined by (x, y, z) ∈ ZZ 3 so that: ( ) The theoretical normal vector is the vector starting from the center of the sphere through the centers of the voxels.

EFFECTS OF THE NEIGHBORHOOD SIZE
For this test, an object whose normal vector field is well known but relatively complex was chosen.The surface equation of this object: x 6 + y 6 + z 6 = 35 6 allows to define the theoretical normal vector field on each point of its surface as follows: V th (x, y, z) = 6x 5 .i + 6y 5 .j + 6z 5 .k (5) (i, j, k) being the canonical base.
To obtain a titled cube, a rotation of the object and of his theoretical vector field was done.To evidence the effect of the choice of r max , we computed the normal vectors of this object for r max = 0 to 24 voxels.Here are the average errors and standard deviation of the errors obtained:

DISCUSSION
The results obtained for spheres (Fig. 2) are close to commonly obtained results in others works (Lenoir et al., 1996;Papier and Françon, 1998).We can note that for spheres of radius < 10, errors are significant.This can easily be explained by the fact that a sphere of small dimensions is not a sphere but a polyhedron.
Rendering images point out the validity of the method on both convex and concave sharp edges (Fig. 3).The relevance of the computed vector field on edges allows to display facetted snow grains with accuracy (Fig. 4).
The choice of the maximal radius of the neighborhood depends on two constraints (Fig. 5): -the normal vector needs to be computed from a sufficient number of gradient vectors.Otherwise, it is corrupted by digitization effects (Fig. 6b).
-the maximal size of the neighborhood has to be close to the size of the smallest detail we want to process correctly.Otherwise, the vector orientation may be altered by irrelevant contributions (Fig. 6d).
Besides, retained surface voxel of the search neighborhood should remain connected together; this is mostly the case for snow at a resolution of 10 microns.
The choice of a maximal radius close to 4-5 voxels, used for all images, seems a good compromise because of the two above constraints.
We have presented a new and simple computational method to obtain accurate normal vectors on the surface of a numerical object.This method allows to distinguish a surface that contains edges from a smooth one.Thus, relevant vectors can be obtained in both cases.
A way to improve this algorithm is to make the calculation as independent as possible of r max .This could be achieved by summing gradient contributions on the largest "well-balanced" neighborhood.A method for the determination of such a neighborhood is now in progress.
P ∈ S and Q ∈ O, we name b the angle (grad(d b (Q)), grad(d b (P))).If b > (π -a 0 )/2, then the segment [PQ] intersects a singularity of the gradient map (such a singularity corresponds to a local extremum of the distance map).Thus, the local gradient whose angle b is wider than b 0 = (π -a 0 )/2 has to be ignored in the summation.

Fig. 1 .
Fig. 1.Summation of gradient vectors (in 2D).For each voxel, the orientation of the gradient vector is represented by a bar in a circle.If (grad(d b (Q)), grad(d b (P))) > b 0 , then [PQ] intersects a singularity of the gradient map.In this case, the contribution of Q has to be ignored.

Fig. 5 .
Fig. 5. Variation of the angular error with the size of the neighborhood for a surface of equation x 6 + y 6 + z 6 = 35 6 .