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.