For practical reasons, one cannot choose a value of
d
and proceed to calculate the corresponding distribution of angles
between the normals at the end of the chords, because of the discrete
nature of the data. Instead, all possible pairs of voxels are considered
and the length of the chord they define is computed as well as the angle
between the corresponding normals. The computed values are used to
create a 2D accumulator array the (
i
,
j
) element of which contains the number of pairs found with chord length
in the range
i
and angle between the normals in the range
j
. If
is the maximum possible length of chord considered, and
N
is the number of bins used for the chord length, then range
i
is defined as
for
. If
M
is the number of bins used to accumulate the values of angle
, range
j
is defined as
for
.
The 2D histogram created in this way is subsequently normalised to be
used as the probability density function of the sampled
pairs that characterise the surface.
Three other issues that arise when we are dealing with discrete data are:
There are two options in 3D: to use either 6, or 26 connectivity. The
surface voxels are identified by using a
neighbourhood (26 connectivity) or a 3D cross neighbourhood (6
connectivity): each voxel is visited in turn: If it has at least one of
its neighbours labelled as background, it is considered as a surface
voxel. 6 connectivity creates thin surfaces made up of voxels that may
share only a common edge. 26 connectivity creates thick robust surfaces
which are densely sampled, because they are made up of voxels that share
common faces. This is the connectivity we shall use in our examples.
The normal to the surface can be calculated in a variety of ways. For example, with the help of a 3D differentiating filter (eg [ 17 ]) applied in the three orthogonal directions, the three components of the surface normal can be identified, as long as the interior of the volume is filled with a label/grey value in sufficient contrast with that used for the exterior voxels. Alternatively, one can consider a small circular neighbourhood around each surface voxel and fit a plane trough all the surface voxels in that neighbourhood using the least square error method. This is the local tangent plane and the surface normal is the normal to that plane. In our experiments the radius of the neighbourhood was chosen to be 4 units. To avoid including in the calculation of the tangent plane voxels which do not belong to the local surface patch, but are included in the neighbourhood due to the high non-convexity of the surface, a standard surface growing method is used locally to identify the right voxels. In yet a third approach, the point where the normal has to be computed is taken to be the vertex of an approximately equilateral triangle made up from voxels around it. In our experiment the side of the triangle was approximately 4 units. All possible such triangles with one vertex at the point of interest are considered and for each triangle its normal is computed. The average of all these normal vectors is assigned to be the normal at the point of interest. We experimented with a spherical surface of radius 60 units. We calculated the normal to the surface using all three methods and constructed the chord-angle histogram described earlier. We then constructed the same histogram by using the analytic formula that describes the surface and compared the results by calculating the mean square error over all the bins. For the 3D filter the error was 17.3%, for the plane fitting 1.6% and for the triangles method 3.5%. In all experiments that follow we adopted the plane fitting method.
When the normal line to the plane has been calculated, we choose which way the normal vector points by examining the label of the first voxel the normal line meets as it leaves the surface and taking advantage of the knowledge of the local surface patch identified by the surface growing algorithm. This is an easy and effective way to decide which way is out of the object and which way is the interior.
Maria Petrou