The Automatic Extraction of the Extremal Points

In practical cases, e1 and e2 can be computed for each point of the 3D image with the equations described in  directly from the differentials of the intensity function of the image f We compute these differentials with linear filtering, using the convolution of the discrete image with the differentials of the Gaussian function !la . The normalization of these filters is not straightforward; we use the responses to simple polynomials, as proposed in . We choose the Gaussian function because it is isotropic, a prerequisite if we are looking for geometric invariants for rigid transformations. Different values of a can be chosen, depending on the level of noise in the 3D images. Changing a is somewhat equivalent to changing the scale at which we look for extremal lines and points.

The hypothesis that the isosurfaces are a good representation of the surface of organs for the case of medical images is a prerequisite: Sometimes, the isosurface can be extracted directly from the 3D image, such as the skin surface in magnetic resonance imaging (MRI) or the bones in X-ray scanner images. For other soft tissues, such as for the brain surface, a presegmentation step is required to isolate the brain from the rest of the data. This can be done with a combination of mathematical morphological operators, filtering, and the search for connected parts or with an automatic "surface edge" extractor, such as the zero-crossing of the image Laplacian. In all cases, the final step of the segmentation is performed using isosurface techniques.

Computation of the Extremal Points in an 8-Voxel Cell

One solution to get the set of extremal points of the 3D image is to compute e1 and e2 for all the voxels ofthe 3D image and then to consider individually each cubic cell, formed with 8 voxels (8-cell), as shown in Fig. 3. There are therefore three values defined for each vertex of the cube: f, e1, and e2. The extremal points in that 8-cell are defined as the intersection of the three implicit surfaces f = I, e1 = 0, and e2 = 0- The method varies according to the type of interpolation or convolution function used to extend continuously the three values at the vertices of the cubic cell to the entire cell- The trilinear interpolation is a good first-order approximation-

The extraction of a polygonal approximation of the crest lines with some warranties about the topology and the orientation of the reconstructed 3D curves is presented with the marching Lines algorithm - Its extension to the extraction of the extremal points was performed in - We briefly recall here the method on a very simple example where the isosurface is a triangle in the cell- This operation can be extended to any configuration of the values of / and e1 while ensuring that the extracted segments form a continuous and closed 3D curve (except when / or e1 is not defined, for instance, at the borders of the image)- The original algorithm also considers orientation problems, which allows us to distinguish between minimum and maximum extremal points-

The first step (Fig- 3, left) is to extract the isosurface within the cell- The isosurface intersects the edges on the cell with the value I- Computing, by linear interpolation along the edges, the points where / = I, we get the three points {Q1; Q2, Q3}- Since we are using a trilinear interpolation within the cell, the intersection of the isosurface with the cell is the triangle

In the second step (Fig- 3, middle), we compute the values of e1 for {Q1, Q2, Q3}, by linear interpolation along the edges of the cubic cell- If they have the same sign, there is no extremal line of e1 in this cell- Otherwise we look for the two points along the triangle edges where the interpolated value of e1 is null: We get the two points {P1; P2}, which form a segment-This is the approximation of the extremal line within the cell-

The last step (Fig- 3, right) is to compute the position of the extremal point- Since P1 and P2 lie on the surface of the cell, we compute the value of e2 at these points with a bilinear interpolation of e2 in the faces- If the two values have the same FIGURE 3 Extraction of the extremal points- An empty circle denotes a positive value, whereas a filled circle indicates a negative one- Adapted from J--P- Thirion- New feature points based on geometric invariants for 3D image registration- International Journal of Computer Vision, 18(2): 121-137, 1996, with permission-

FIGURE 3 Extraction of the extremal points- An empty circle denotes a positive value, whereas a filled circle indicates a negative one- Adapted from J--P- Thirion- New feature points based on geometric invariants for 3D image registration- International Journal of Computer Vision, 18(2): 121-137, 1996, with permission-

sign, there is no extremal point on this cell. Otherwise, as is shown here, we compute its position P by interpolating the zero value along the segment.

Randomized Extraction of Extremal Points

Of course, we could look for extremal points in all the possible cells of the image, excepting regions of null gradient and umbilics. However, it is much more efficient to randomize the search: We start with seed cells, randomly chosen in the 3D image, and discard the ones for which the sign of f — I is the same for all the vertices. Then we compute the values of ex for the 8 vertices of the cell. Once again, a simple test discards the cells that are not crossed by a k1 extremal line (the sign of e1 is the same for the 8 vertices). If there is an extremal line, we extract it from end to end, using the marching lines algorithm (we follow the extremal line "marching" from one cell to the next).

At each point of the polygonal approximation of the crest line, we compute the second extremality e2 by bilinear interpolation. If there is a sign change, we compute the extremal point on the segment of the extremal line that we are currently following.

The randomized implementation of the marching lines allows us to extract the main extremal lines (i.e., the longest ones, which experimentally appeared to be the most reliable ones) of the 3D image, with only very few seeds (with respect to the total number of voxels), randomly distributed in the 3D images. The probability of missing an extremal line is approximately proportional to the inverse of its length. This method reduces drastically the number of computations to perform, compared to the extensive implementation: Typically, one uses 10% of the number of voxels as seeds. Even if the set of generated extremal points is not complete, it is generally sufficient to perform a reliable 3D registration.