Calculating Surface Normals and Curvatures
Surface Normals and Curvatures
Starting with the range map
Range map visualization, and many types of automated analysis, require the calculation of local surface orientation (the normal vector) and also — ideally — the accurate detection of surface discontinuities, or edges.
The range map is X — recall that this is a twodimensional array of threedimensional vectors, each representing a local (x,y,z) measurement if such a measurement could be made, or identical to (0,0,0) to indicate an invalid point where no measurement could be made. We will refer to the (x,y,z) measurement at range map index row, col as X_{row,col}.
Traditional approachs make a number of inappropriate assumptions about range data:

Adjacency relationships in the range map correspond to
adjacency relationships in the real world.
This really implies three deeper related assumptions:
 Jump discontinuities do not exist.
 Small local extrema do not exist.
 Sensor sampling was dense enough to capture all surface features.
 The spatial sampling rates along both columns and rows are roughly equal and constant throughout the range map.
The Surface Normal
The surface normal is a unit vector perpendicular to the
local surface.
We can calculate a map of surface normals, a twodimensional array
of surface normal vectors, where
S_{row,col}
is the surface normal at range map index row,col.
The surface normal indicates how the measured surface location varies
while indexing across the range map,
and so it is the first derivative of position:
S = dX/dindex
where index is some measure of location within the map.
The Curvature
The curvature indicates how the surface orientation changes
as one traverses the surface.
Put another way, it indicates how the surface normal change with index
across the range map.
Therefore the curvature is the derivative of the surface normal,
and thus the second derivative of position:
C = dS/dindex = d^{2}X/dindex^{2}
Traditional Methods of Calculating Surface Normal and Curvature
Traditional methods for calculating surface normal and curvature start by convolving the range map X with various 3x3 kernels, then combining the resulting terms. 5x5, 7x7, or even larger kernels could be used to reduce noise; this is effectively the same as first smoothing X with a 3x3, 5x5, or larger smoothing kernel, and then applying the convolutions that calculate partial derivatives.
A variation on this is to observe that each range point belongs to 9 3x3 patches (or to 25 5x5 patches, or 49 7x7 patches, or ...), and to calculate the normal and curvature for each of them. Then, for the point in question, use the measures from the patch that is most planar. If the target is known to be made up of roughly planar regions, this modification may help. But smoothly curving surfaces end up with a faceted appearance if this assumption is inappropriately applied. The result is a "contour map" effect, as seen below in the dorsal view of a fish skull (the skull is seen from above, with the nose of the fish to the left). Note the "contour map" effect in the nearly flat region at the top of the head behind the eye sockets, and the general graininess of the rendered surface.
The Signal Processing View, and What is Wrong With Traditional Techniques
Since normals and curvatures are first and second derivatives, respectively, we must keep in mind that a derivative is a highpass filter. A second derivative is a highpass filter with a faster rolloff.
If we are really discussing filters, it would be appropriate to ask about the cutoff frequency (or wavelength) and the rolloff. For instance, "a highpass filter with 3dB cutoff frequency of 10.8 MHz and a rolloff of 16 dB per decade", as we might specify a conventional filter.
We cannot answer these questions in any meaningful way!
The spatial sampling rates along rows and columns are unequal:
dX/drow != dX/dcol
The spatial sampling rates along either rows or columns varies
throughout the image, so in some sense dX/drow
is not equal to dX/drow when measured at
different parts of the range map.
Furthermore, the traditional methods of calculation were devised in a time when range maps were only as large as absolutely necessary. In other words, the data sampling was as coarse as possible to achieve the goal with a minimum of data collection and processing. With increasing computation and storage capabilities, modern data frequently has dX/drow and dX/dcol close to the spatial sampling limits of sensor geometry and optics.
Applying a 3x3 convolution to finely sampled data yields very localized measures. In signal processing terms, the rolloff frequency is very high. Speaking practically, or esthetically, the result will be of little use as it will be dominated by smallscale highfrequency noise, as seen below.
The Solution for Calculating Surface Normals
First, we must calculate local indexbased "derivatives",
which of course are not really derivatives at all but merely differences.
δX/δrow =
X_{row+1,col}  X_{row,col}
δX/δcol =
X_{row,col+1}  X_{row,col}
Then, calculate local estimates of dX/drow and dX/dcol based on weighted averages of δX/δrow and δX/δcol. These weighted averages are calculated over collections of points that lie within some specified distance of the point in question. The weight varies linearly with distance from that point. The specified distance, or radius of interest, is based on the scale of the features of interest.
In words, to calculate dX/drow for some specific
r,c index in the range map,
let us select some range of interest R,
where R is approximately the size of the target features
that interest us.
Next, find the set of points in the range map within R of the
point at index r,c.
Note that this distance is in realworld coordinates, not in the indices
of the range map, and it consists of finding all range map indices i,j
where:
 X_{r,c} 
X_{i,j}  < R
For each point i,j we must calculate a scaling factor S:
S = 1.0   X_{r,c} 
X_{i,j} / R
In other words, the contribution of points decreases linearly with
distance, reaching zero at distance R
We then sum the scaled contributions of δX/δrow and δX/δcol, dividing each by the sum of the scaling factor S, to estimate dX/drow and dX/dcol on a physically meaningful scale. Note that these derivatives are vectors:
for (each row, col index) {
dX_{row,col}/drow = 0
dX_{row,col}/dcol = 0
sum = 0
for (each i, j index where  X_{r,c}  X_{i,j}  < R) {
S = 1.0   X_{r,c}  X_{i,j} /R
dX_{row,col}/drow = dX_{row,col}/drow + δX_{i,j}/δrow
dX_{row,col}/dcol = dX_{row,col}/dcol + δX_{i,j}/δcol
sum = sum + S
}
dX_{row,col}/drow = (dX_{row,col}/drow)/sum
dX_{row,col}/dcol = (dX_{row,col}/dcol)/sum
}
Given these derivatives
dX/drow
and
dX/dcol,
we can calculate the local surface normal N as their cross product:
N_{i,j} = dX_{i,j}/drow x dX_{i,j}/dcol,
The Solution for Calculating Curvatures
The curvature C can be estimated as follows, where "·" specifies
the dot product operator:
d^{2}X_{r,c}/drow^{2} =
( X_{r+1,c}  X_{r,c})
· N_{i,j} 
( X_{r,c}  X_{r1,c})
· N_{i,j}
d^{2}X_{r,c}/dcol^{2} =
( X_{r,c+1}  X_{r,c})
· N_{i,j} 
( X_{r,c}  X_{r,c1})
· N_{i,j}
This is noisy, so again it should be smoothed over linearly scaled regions
based on actual range from the point of interest,
as was done to calculate normals.
Example Results
The example target is a cuneiform tablet.
Overall tablet dimensions are approx 30x30 mm.
Lines of cuneiform text are spaced approximately 710 mm apart.
Individual signs are approximately 57 mm on a side. Wedges are approximately 1.2 x 5 mm. Circular pits are approximately 3.5 mm in diameter.
The following are rendered images based on normals calculated in different ways. First, using the traditional method, showing the graininess. Then at a variety of scales corresponding to the scales of sign components, signs, and the overall tablet.
Original method — note graininess
R = 0.20 mm, approximate scale of sign components
R = 1.00 mm, approximate scale of signs
R = 3.00 mm, individual signs are lost and overall shape is shown
Further Details Regarding Curvature
The traditional method calculated d^{2}X/drow^{2} and d^{2}X/dcol^{2}, and then combined them. This assumed that interesting features were aligned with either the range map rows or columns. But how should we really do this?
Our method is to calculate four terms — along rows, along columns, and along both diagonals:
 d^{2}X/drow^{2}
 d^{2}X/drowdcol
 d^{2}X/dcoldrow
 d^{2}X/dcol^{2}
Some features will only appear in one of these four measures.
As a single curvature measurement, we could select at each row,col index the one measure with the greatest magnitude. Or, we could simply use the average of the four. While the results are very slightly different, there seems to be no significant difference in data quality. Below are the results for curvature measurements at scales of 0.3 mm, 0.5 mm, 1.0 mm, and 3.0 mm.
R = 0.30 mm
R = 0.50 mm, approximate scale of sign components
R = 1.00 mm, approximate scale of signs
R = 7.00 mm, individual signs are lost and overall shape is shown
Further Results
Below are some images that suggest that the curvature measurements described here can highlight some detail that would unseen in a rendered image. A rendered image must specify some lighting vector, and one combination of lighting and view vectors will always leave some 3D features unseen.
The curvature image (at left) clearly shows that there are four semicircular sign components, while the rendered image (at right) is a bit confusing.
The curvature image (at left) clearly shows that the leftmost column has three wedges, while the middle wedge is difficult to see in the rendered image (at right).
While the curvature image (at left) looks a bit "noisy", it does reveal some sign components that are difficult to see in the rendered image (at right).
Conclusions
The methods described here are more computationally expensive. However, this is of little concern with modern computer speeds, and it will only become less of a concern over time.
It appears that this method would be worth the computational expense for certain applications:
 Accurate visualization of densely sampled range maps
 Cueing a human user to find areas of interest
 Automatically selecting points of interest to be used for automated registration (e.g., ICP)