Xavier Pennec

Nicholas Ayache

Jean-Philippe Thirion

Landmark-based registration using features identified through differential geometry

Registration of 3D medical images consists in computing the "best" transformation between two acquisitions, or equivalently, determines the point to point correspondence between the images. Registration algorithms are usually based either on features extracted from the image (feature-based approaches) or on the optimization of a similarity measure of the images intensities (intensitybased or iconic approaches). Another classification criterion is the type of transformation sought (e.g. rigid or non-rigid).

In this chapter, we concentrate on feature-based approaches for rigid registration, similar approaches for non-rigid registration being reported in another set of publication [START_REF] Subsol | Crest lines for curve based warping[END_REF][START_REF] Subsol | A general scheme for automatically building 3D morphometric anatomical atlases: application to a skull atlas[END_REF]. We show how to reduce the dimension of the registration problem by first extracting a surface from the 3D image, then landmark curves on this surface and possibly landmark points on these curves. This concept proved its efficiency through many applications in medical image analysis as we will see in the sequel. This work has been for a long time a central investigation topic of the Epidaure team [START_REF] Ayache | Medical computer vision, virtual reality and robotics -promising research[END_REF] and we can only reflect here on a small part of the research done in this area.

We present in the first section the notions of crest lines and extremal points and how these differential geometry features can be extracted from 3D images. In Section 2, we focus on the different rigid registration algorithms that we used to register such features. The last section analyzes the possible errors in this registration scheme and demonstrates that a very accurate registration could be achieved.

Features Extraction: Extremal Points and Lines

To extract reliable curves on a surface, most approaches try to generalize the notion of "edges" to smooth surfaces to find the most salient features on the surface: ridges. Prior to the late 1980s and early 1990s, the interest in ridges was mostly theoretical, in areas of mathematics related to catastrophe theory [START_REF] Porteous | Ridges and Umbilics of Surfaces[END_REF][START_REF] Markatis | Some Generic Phenomena in Families of Surfaces in R 3[END_REF][START_REF] Koenderink | Solid Shape[END_REF][START_REF] Porteous | Geometric Differentiation: For the Intelligence of Curves and Surfaces[END_REF][START_REF] Anoshkina | Ridges and Ravines on a Surface and Related Geometry of Skeletons, Caustics and Wavefronts[END_REF]. Crest lines are then defined as the cuspidal edges of a caustic surface, and the link between caustics and curvatures on a surface was established. Practical applications were then discovered by researchers in computer vision, graphics, and medical imaging together with the specification of algorithms to extract ridges. In the following, we are going to focus on the crest and extremal line as introduced in [START_REF] Monga | From voxel to intrinsic surface features[END_REF][START_REF] Monga | Using Partial Derivatives of 3D Images to Extract Typical Surface Features[END_REF] and developed in [START_REF] Thirion | Computing the differential characteristics of isointensity surfaces[END_REF][START_REF] Thirion | The 3D marching lines algorithm[END_REF]. Basically, these curves are (subsets of) the loci of the surface where one of the principal curvatures reaches a locally optimum in the associated principal direction. In these works, the crest lines are extracted using third order derivatives of the image intensities. An alternative approach was presented in [START_REF] Guéziec | Surface representation with deformable splines: Using decoupled variables[END_REF] with the use of a B-spline approximation of the surface.

A different notion of ridges is given in [START_REF] Van Den Elsen | Automatic Registration of CT and MR Brain Images Using Correlation of Geometrical Features[END_REF][START_REF] Eberly | Ridges in Image and Data Analysis[END_REF]: they are defined as the salient structures of the intensity surface defined by I = f (x, y, z). Here, the ridges are surfaces and are more like results of the medial axis transform than the intuitive notion of salient lines on a surface. Cutting et al [START_REF] Cutting | Applications of computer graphics to the evaluation and treatment of major craniofacial malformation[END_REF] also developed a method using a template of ridges to assist in their extraction from image data. This method was extended by Dean et al [START_REF] Dean | Three dimensional MR-based morphometric comparison of schizophrenic and normal cerebral ventricles[END_REF]. A good review of the many definitions of ridges can be found in [START_REF] Belyaev | Principal direction ridges[END_REF].

Definition and Properties Differential Geometry of 3D Surfaces

Let us first recall briefly some results of differential geometry about surface curvatures (a good introduction to these notions can be found in [START_REF] Manfredo | Differential Geometry of Curves and Surfaces[END_REF] or in [START_REF] Koenderink | Solid Shape[END_REF]). In this paper, we call a smooth surface a 3D surface which is continuously differentiable up to the third order. At any point P of such a 3D surface, we can define one curvature per direction t in the tangent plane of the surface. This directional curvature k t is the curvature of the 3D curve defined by the intersection of the plane (P, t, n) with the surface, where n is normal to the surface.

Except for the points where this curvature k t is the same for all the directions t, which are called umbilic points, the total set of curvatures can be described with only two privileged directions, t 1 and t 2 , and two associated curvature values, k 1 = k t1 and k 2 = k t2 , which are called respectively the principal directions and the associated principal curvatures of the surface at point P , as shown in Figure 1. These two principal curvatures are the extrema of the directional curvatures at point P , and (except for umbilic points) one of these two is maximal in absolute value, let us say k 1 : we call this the largest curvature, in order not to be mistaken with the maximal curvature. We simply call the second (principal) curvature the other principal curvature k 2 .

Extremal Lines

The crest lines are intuitively the loci of the surface where the "curvature" is locally maximal. More precisely, we define them as the loci of the surface where the largest curvature, k 1 , is locally maximal (in absolute value) in the associated principal direction t 1 . In [START_REF] Monga | Using Partial Derivatives of 3D Images to Extract Typical Surface Features[END_REF], it is shown that these points can be defined as the zero-crossing of an extremality function e, which is the directional derivative of k 1 in the direction t 1 .

We have proposed another method to compute them in [START_REF] Thirion | The 3D marching lines algorithm[END_REF][START_REF] Thirion | Computing the differential characteristics of isointensity surfaces[END_REF], for the case of iso-intensity surfaces. Our method is based on the use of the implicit functions theorem. Basically, we have shown that the crest lines can be extracted as the intersection of two implicit surfaces f = I and e = 0, where f represents the intensity value of the image, I an iso-intensity threshold, and e = ∇k 1 • t 1 is the extremality function (see Figure 2, left). We have proposed an algorithm, called the Marching Lines, to automatically extract these crest lines. This algorithm can also be used to overcome some orientation problems (mainly due to the fact that the principal directions are only directions and not oriented vectors), by locally orienting the principal directions along the extracted lines.

In fact, for each point of the surface, two different extremality coefficients can be defined, corresponding to the two principal curvatures:

e 1 = ∇k 1 • t 1 and e 2 = ∇k 2 • t 2 (1)
We found experimentally that the maxima (in absolute values) are more stable landmarks than the minima: crests or rifts (maxima) are stable, whereas the loci in a valley where the ground floor is the flattest (minima) are very sensitive to small perturbations in the data.

We call extremal lines all the lines defined as the zerocrossings of either e 1 or e 2 . There is therefore four major different types of extremal lines, depending of whether the corresponding curvature is the largest or the second one and whether it is a local maximum or minimum. Furthermore, the signs of the largest and second curvatures help to distinguish between four additional sub-types of extremal lines, leading to a classification into 16 types. The crest lines are two of them: positive largest curvature maxima (k 1 > 0 and ∇e 1 • t 1 < 0) and negative largest curvature minima (k 1 < 0 and ∇e 1 • t 1 > 0).

Extremal Points

We now define the extremal points as the intersection of the three implicit surfaces: f = I, e 1 = 0 and e 2 = 0. The notions of extremal lines and extremal points are closely related to the notion of corner points, in 2D images, as defined in [START_REF] Kitchen | Gray-level corner detection[END_REF], [START_REF] Noble | Finding corners[END_REF], and [START_REF] Deriche | A computational approach for corner and vertex detection[END_REF]. A study of the evolution in 2D of corner points with respect to the scale can be found in [START_REF] Gauch | Multiresolution Image Shape Description[END_REF]. A similar study on the scale-space behavior of the extremal lines and the extremal points was presented in [START_REF] Fidrich | Multiscale analysis of invariants: Application to volumetric medical image processing[END_REF].

Extremalities e 1 and e 2 are geometric invariants of the implicit surface f = I: they are preserved with rigid transforms (rotations and translations of the object). Therefore, the relative positions of the extremal points are also invariant with respect to a rigid transformation, i.e. for two different acquisitions of the same subject. There are 16 different types of extremal points, depending on the type of extremality: local minimum or maximum of the extremalities e 1 and e 2 and the signs of k 1 and k 2 . This classification can be used to reduce the complexity of the matching algorithm.

However, the intuitive interpretation of extremal points is not straightforward. The extremal lines are 3D curves, for which we are able to compute the curvature, but the extremal points are generally not the points of the extremal lines whose curvature is locally maximal. Even if they are not extremal curvature points, the extremal points are very well defined, and there is no reason for their locations along the extremal lines to be less precise that the lines positions themselves, because the precision of the computation of k 1 and k 2 is approximately the same.

Geometric Characteristics

Let us begin with the points on a surface. We have already seen (Figure 1, left) that any such point could be provided with a trihedron (n, t 1 , t 2 ) formed by the normal to the surface and the two principal directions. As our points are also on extremal lines, we could provide them with the differential characteristics of 3D curves (Figure 1, right), i.e. the Frénet trihedron (t, n c , b), where t is the tangent to the curve, n c its normal and b the binormal.

These two trihedrons are not the same as the extremal lines are generally not lines of curvature. However, as the curve is embedded in the surface, the tangent to the curve t is constrained to be in the tangent plane of the surface spanned by (t 1 , t 2 ). Thus, there are two independent parameters characterizing the relative configuration of the trihedron: we can measure two angles θ = (t, t 1 ) and φ = (n c , n). These characteristics are invariant with respect to rigid transformations.

Two other invariants come from the surface (principal curvatures k 1 and k 2 ). One could also think to add the curvature k, the torsion τ of the curve and the geodesic torsion τ g of the curve with respect to the surface but it appears that k and τ g are completely determined by the surface invariants: k cos φ = k 1 cos 2 θ + k 2 sin 2 θ and τ g = (k 2 -k 1 ) cos θ sin θ. Thus, we are left with the torsion of the curve.

However, the computation of the Frénet trihedron (t, g, b) and the curve torsion τ has to be done on the extremal curve itself after its extraction. If this can be done directly on the polygonal approximation, a much better method is to compute the characteristics on a local B-spline approximation of the curve [START_REF] Guéziec | Smoothing and matching of 3d space curves[END_REF].

Chap. 31

Landmark-Based Registration Using Features Identified Through Differencial Geometry 502

The Automatic Extraction of the Extremal Points

In practical cases, e 1 and e 2 can be computed for each point of the 3D image with the equations described in [START_REF] Thirion | New feature points based on geometric invariants for 3D image registration[END_REF] 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 e -x 2 /2σ 2 . The normalization of these filters is not straightforward; we use the responses to simple polynomials, as proposed in [START_REF] Monga | Using Partial Derivatives of 3D Images to Extract Typical Surface Features[END_REF]. We choose the Gaussian function because it is isotropic, a prerequisite if we are looking for geometric invariants for rigid transformations. Different values of σ can be chosen, depending on the level of noise in the 3D images. Changing σ is somewhat equivalent to changing the scale at which we look for extremal lines and points.

The hypothesis that the iso-surfaces are a good representation of the surface of organs for the case of medical images is a prerequisite: sometimes, the iso-surface can be extracted directly from the 3D image, such as the skin surface in Magnetic Resonance Image (MRI) or the bones in X-ray scanner images. For other soft tissues, such as for the brain surface, a pre-segmentation 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 zerocrossing of the image Laplacian. In all cases, the final step of the segmentation is performed using iso-surface techniques.

Computation of the Extremal Points in a 8-Voxel Cell

One solution to get the set of extremal points of the 3D image is to compute e 1 and e 2 for all the voxels of the 3D image and then to consider individually each cubic cell, formed with 8 voxels (8-cell), as shown in Figure 3. There are therefore three values defined for each vertex of the cube: f , e 1 , and e 2 . The extremal points in that 8-cell are defined as the intersection of the three implicit surfaces f = I, e 1 = 0, and e 2 = 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 tri-linear 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 line algorithm [START_REF] Thirion | The 3D marching lines algorithm[END_REF]. Its extension to the extraction of the extremal points was performed in [START_REF] Thirion | Computing the differential characteristics of isointensity surfaces[END_REF]. We briefly recall here the method on a very simple example where the iso-surface is a triangle in the cell. This operation can be extended to any configuration of the values of f and e 1 while ensuring that the extracted segments form a continuous and closed 3D curve (except when f or e 1 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 (Figure 3, left) is to extract the iso-surface within the cell. The iso-surface intersects the edges on the cell with the value I. Computing, by linear interpolation along the edges, the points where f = I, we get the three points {Q 1 , Q 2 , Q 3 }. Since we are using a tri-linear interpolation within the cell, the intersection of the iso-surface with the cell is the triangle

{Q 1 , Q 2 , Q 3 }.
In the second step (Figure 3, middle), we compute the values of e 1 for {Q 1 , Q 2 , Q 3 }, by linear interpolation along the edges of the cubic cell. If they have the same sign, there is no extremal line of e 1 in this cell. Otherwise we look for the two points along the triangle edges where the interpolated value of e 1 is null: we get the two points {P 1 , P 2 } which form a segment. This is the approximation of the extremal line within the cell.

The last step (Figure 3, right) is to compute the position of the extremal point. Since P 1 and P 2 lie on the surface of the cell, we compute the value of e 2 at these points with a bi-linear interpolation of e 2 in the faces. If the two values have the same 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 e 1 for the 8 vertices of the cell. Once again, a simple test discards the cells which are not crossed by a k 1 extremal line (the sign of e 1 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 e 2 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.

Example of Extracted Extremal Lines and Points

In Figure 4, we can see an example of the lines extracted automatically (with a manual choice of the iso-intensity threshold) in a CT image of a dry skull. Some of the 550 crest lines may be recognized as anatomical landmarks, such as the orbits or the inferior border of the mandible. The lines are colored by the sign of the e 2 extremality. Thus, extremal points are located at the color changes along the lines. There are around 3000 such extremal points.

In an MR image, the surface of the brain is not very well defined by an iso-intensity of the image. A presegmentation step is usually needed 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 Figure 5, we used a segmentation of the surface of the brain and extracted the crest lines on this surface. Lines in red (with a positive largest curvature) roughly correspond to sulci whereas blue lines (with a negative largest curvature) could be interpreted as gyri.

Chap. 31

Landmark-Based Registration Using Features Identified Through Differencial Geometry 504

Rigid Registration

Let us now consider two images of the same modality and of the same patient, but in a different position. We extract extremal lines on both images. The problem is to put into correspondence the two sets of lines (the model and the scene), which is often called the matching step, and to compute the best rigid transformation that superimposes the matched lines.

It is important to note that a global registration algorithm, for instance superimposing the barycenters of all points and the inertia axes, will often fail due to the occlusion. Indeed, the images being taken in different positions, the region of interest are frequently different in the two images, leading to crest lines and extremal points present in one image and not in the other. The images noise will also induce the extraction of spurious lines and points in different parts of the two images.

Curve Registration

Several algorithms adapted from computer vision have been proposed and used over time. In [START_REF] Guéziec | Smoothing and matching of 3d space curves[END_REF], Guéziec matches the crest lines using a combination of geometric hashing [START_REF] Lamdan | Geometric hashing: A general and efficient model-based recognition scheme[END_REF] and Hough transform (see for instance [START_REF] Leavers | Survey: which Hough transform?[END_REF]). The basic idea was to index each point of each model crest line in a hash-table using its invariant characteristics. At recognition time, the invariants of each scene crest line point are used to recover, thanks to the hash table, the possible matches with model points (geometric hashing step). For each match (i.e., couple of model and scene points), a rigid transformation is computed by superimposing the Frénet frames and used to index the match in a new accumulator sampling the rigid transformation space (Hough transform step). Densely populated cells in this second accumulator are detected as rigid body transformations that are candidates to match a large number of crest points. For each such cell, a refined least-squares transformation is computed using the matches indexed in this cell.

Extremal Points Registration using Alignment

With the development of completely automated methods to extract crest lines and the higher resolution of images, the number of crest lines drastically increased, leading to a much higher density of invariants in the hash table. This could lead to an important number of false positives that would overwhelm the correct matches. The maximum complexity would then be reached and the algorithm could even provide a wrong answer. To address this problem, Thirion reduced once again the image information by keeping only a very small number of specific points on the crest lines: the extremal points. Typically, they represent only 16% of the number of crest line points, but we are still left with 2000 to 5000 points in each image. Thirion used in [START_REF] Thirion | New feature points based on geometric invariants for 3D image registration[END_REF] another computer vision based technique: alignment (or prediction-verification) [START_REF] Ayache | Hyper: A new approach for the recognition and positioning of twodimensional objects[END_REF][START_REF] Huttenlocher | Object recognition using alignment[END_REF]. The basic idea is to compute the registration of each triplet of model points with each triplet of scene points, superimpose the two sets of points using this transformation and verify this match using an iterative closest point algorithm (see Section 2.4). However, the search for com-

Handbook of Medical Imaging -Processing and Analysis

IV Registration patible triplets in the model and the scene can be reduced since there are some unary invariants (the principal curvatures k 1 and k 2 ), secondary invariants (e.g. the distance between the two points, or the relative orientation of the surface normals and the principal curvatures) and even ternary invariants (involving the whole triplet of extremal points). Thirion used 18 such invariants for each triplet, pre-computed and stored them in a hash table to retrieve in constant time the compatible triplets. Thus, the complexity of the algorithm is O(n 4 ) since there are n 3 triplets, and a verification of O(n) for each triplet match.

In practice, this complexity is not reached as we can stop as soon as a given number of points is matched after verification (typically 10 %).

Substructure Matching with Frame Features

We came back in [START_REF] Guéziec | Medical image registration using geometric hashing[END_REF] to geometric hashing, but the idea was to use all the geometric information on the features while taking great care of the uncertainty handling for the algorithm to be robust (see [START_REF] Pennec | Toward a generic framework for recognition based on uncertain geometric features[END_REF] for an analysis of recognition algorithms with uncertain geometric features). In addition to the point's position, we can add the normal vector n and the two principal directions t 1 and t 2 of the surface to constitute a local coordinate system, or a frame. In this context, each medical image is modeled by a set of frames and the matching problem is to find the correspondences between two subsets of frames that are in the same configuration in the two "images", up to a global rigid transformation.

Invariant Representation: Preprocessing

Step To obtain an invariant representation with respect to the global position and orientation of the considered structure, we can express the configuration of all frames relative to one frame (called the basis). For efficiency, this representation is stored in a hash table and, for correctness, we include the uncertainty of each invariant. As only part of the frames are in the same configuration in the two images, the one chosen as the basis may not be present in the other image. The preprocessing step is thus repeated with each frame as the basis.

Recognition

Step Choosing a frame of the second structure (the scene) as the basis, we compute the invariant representation and retrieve, thanks to the hash table, what are the compatible model frame couples. If the basis belongs to a common substructure, then a significant number of frames are in the same configuration with respect to it. We then match the model and scene bases (Fig. 6).

This process is repeated for every extremal point as the basis to find its possible matches in the model and we only keep the matches that are above a given threshold (typically 10 % of the number of extremal points).

Clustering Compatible Matches and Verification

For each individual match, we maintain during the recognition step an estimation of the associated transformation by fusing the transformations between confirming frames. To group matches belonging to the same (rigid) substructure, we run a very simple clustering algorithm on the associated transformation. Indeed, up to measurement errors, frames should undergo a similar transformation within a single substructure. Each cluster is then refined by an iterative closest neighbor technique where we enforce symmetry of the matches and verify their validity with a χ 2 test.

Matching Crest Lines In order to reduce once again the complexity of the algorithm, we exploited in this method the structure of the extremal points: they belong to crest lines. The principle is to consider each model crest line as a different object. Index all model lines in the same hash-table, we retrieve at recognition time the model lines corresponding to each scene crest line.

However, different crest line matches can correspond to different transformations. Thus, we run once again our clustering algorithm on the transformations to find out the compatible line matches and we obtain a single transformation from the model to the scene image.

ICP on Frames

When images are close enough, one can use still another algorithm: the Iterative Closest Point [START_REF] Besl | A method for registration of 3-d shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF]. The basic principle is the following. For each scene point, we look for the closest point in the model with the current transformation, compute a new rigid transformation with these matches, and iterate the process until convergence.

Of course, since we have more geometric information than just the point position, we use a generalization: the Iterative Closest Feature [START_REF] Pennec | Toward a generic framework for recognition based on uncertain geometric features[END_REF]. The idea is to use a higher dimensional space for the closest point search. In our case, the space is made of the extremal point position, the trihedron (n, t 1 , t 2 ), and the unary invariants k 1 and k 2 . The important point is to set an appropriate metric on this space in order to combine efficiently the different units of measurement. In our algorithm, this is done using the inverse of the covariance matrix of the features. This matrix can be re-estimated after convergence and the whole process iterated. However, we did not observe a critical influence of the covariance matrix values, as soon as it is approximately respecting the variation range of the different components.

Chap. 31

Landmark-Based Registration Using Features Identified Through Differencial Geometry 506

Figure 6: Preprocessing: the 6D invariant vector associated with every couple of model frames is computed with its error zone and used as an index for the couple in the hash table. Recognition: for each scene frame couple, we compute the 6D invariant vector and retrieve through the hash table every compatible model frame couple. For each such couple, we tally a vote for the matching of the reference frames (here the match (F m i , F s j ) scores 2).

Examples of Rigid Registrations

Registration of CT Images of the Skull Figure 7 presents an example of the registration of two CT images of the a dry skull in a Plexiglas box in two different positions. We used the geometric hashing algorithm on frames (Section 2.3). As the transformation between the two images is close enough to the identity, the ICP algorithm also gives very similar results. About 75 crest lines are matched with more than 4 extremal points among the 550 lines in each image leading to a total of 550 matched extremal points (only on the 75 matched lines). Using the techniques described in Section 3.2, we have computed that the typical object accuracy (the expected standard RMS error on image super-imposition due to the transformation in the area of the matched features) is 0.04 mm, whereas the typical corner accuracy is 0.1 mm. This is to be compared with the voxel size: 1 x 1 x 1.5 mm.

Registration of MR Images of the Head Figure 8 is an example of the registration of two MR T1 images of the same patient. In this case, 240 crest lines are matched among approximately 2100 in each image, for a total of 860 matched extremal points among 3600 in each image (about 25 %). We used the zero-crossing of the Laplacian to define the interest surfaces. Thus, there are crest lines all over the head. However, if some of the matched lines are located on the surface of the skin (we can recognize the nose and the eyes), most of them are located on the surface of the brain. The typical object accuracy of the registration is 0.06 mm for a typical corner accuracy of 0.125 mm. Once again, the accuracy is far below the voxel size (0.97 x 0.97 x 1.5 mm).

Robustness and Uncertainty Analysis

Once we have registered the images, i.e. found matches and a rigid transformation, the question is: how confident can we be with this result? There are two main types of errors in feature-based registration algorithms. Firstly, the matches could be completely wrong and we simply recognized by chance n features in approximately the same configuration. This is called a gross error in statistics and a false positive in recognition. But even if we got the matches right, the features we are using to compute the registration are corrupted by noise and induce a small error, or uncertainty, on the transformation parameters. In this section, we analyze in turn these two types of error.

Robustness Analysis

Since our features are noisy, we had to allow for a certain error when matching them. In the registration algorithm of Sections 2.3 and 2.4, this is computed from the covariance matrices. The existence of such an error zone allows us to match features that by chance fall in this area. When this probability is sufficiently high, individual false matches can combine themselves (conspiracy) to produce an important matching score. However, we should note that such a false positive is a correct recognition from a feature point of view (a globally consistent matching) but is incorrect from the image or data point of view. This problem is inherent to the ascendent organization of information: some important information can be dismissed by simplifying the image as a set of features.

Chap. 31

Landmark-Based Registration Using Features Identified Through Differencial Geometry 508

In [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF], we developed a generic method to evaluate the probability of matching τ features just by chance. The basic framework is the following. Let us assume for the moment that the transformation between the two images is fixed. First, we compute the selectivity η, which is the probability that a random feature (uniformly distributed in the image) falls in a given error zone. Then, we compute the probability that at least τ of the m scene features fall in one of the n model error zones. In our analysis, computations are simplified by assuming that all features are randomly distributed. Now, we will accept a match if there exists one transformation such that at least τ features are put into correspondence. Thus, to obtain the mean number of false positives, we just have to integrate over all possible transformations. Let d be the "diameter" of the images, we get the following estimation:

P (2πd) 3 3   1 -e -(nmη) τ j=0 (nmη) j j!   .
In the example of Section 3.3, we compute that the selectivity is η pt = 2.10 -6 if we just use the position of the extremal points and η f r = 1.5.10 -8 if we model them using frames. The diameter of the image is d 400mm and we extracted around 2,500 extremal points in each image. We plot in Figure 9 the number of false positives P with these values. The actual matches found involve about 500 features and the probability of being a false positive is thus practically zero. However, we must be careful that the "object" we registered is not always the one we wanted, even if this is definitely not a false positive: there can be several different rigid motions in a single image, for instance the skull and the brain in MR images.

From Feature to Transformation Uncertainty

Here, we assume that the matches are right. However, measurement errors on the features induce an estimation error on the transformation. We developed in [START_REF] Pennec | A framework for uncertainty and validation of 3D registration methods based on points and frames[END_REF][START_REF] Pennec | Feature-based registration of medical images: Estimation and validation of the pose accuracy[END_REF] a method where we register the features, estimate the noise on the frames and propagate this noise to estimate the uncertainty of the rigid transformation.

Feature Uncertainty For extremal points (modeled as frames), we proposed a "compositive" model of noise.

Basically, this means that the noise is assumed to be identical on all extremal points in the local frame (i.e., with respect to the surface normal and the principal directions). This has to be compared with the standard additive noise model on points where we assume an identical noise with respect to the image axes. In the case of the MR images of the next Section, this leads to an interesting observation: we draw in Figure 10 a graphical interpretation of the covariance matrix estimated on the extremal points after registration.

We obtain an approximately diagonal covariance matrix with standard deviations σ t1 σ t2 2 deg, σ n = As far as the trihedron part is concerned, this means that the surface normal is more stable than the principal directions, which is expected since the normal is a first order differential characteristic of the image and the principal directions are second order ones.

For the position, the coordinate along the normal is once again more stable than in the tangent plane for the same reasons. The 3D standard deviation of the position is σ = 1.04, which is in agreement with the additive noise model on points. However, for the additive noise model, the estimated covariance is isotropic. Thus, using an adapted model of noise on frames allows us to extract much more information about the feature noise. This constitutes an a-posteriori validation of our "compositive" model of noise on extremal points.

Transformation Uncertainty Now the problem is to propagate the feature uncertainty to the rigid transformation. Let f represent a rigid transformation and χ the observed data. The optimal transformation f minimizes a given criterion C(f, χ) (for instance the least-squares or the Mahalanobis distance). Let Φ(f, χ) = ∂C(f, χ)/∂f . The characterization of an optimum is Φ( f , χ) = 0. Now, if the data are moving around their observed values, we can relate the new optimal parameters using a Taylor expansion. Let Ĥ = ∂Φ/∂f and ĴΦ = ∂Φ/∂χ be the values of the second order derivatives of the criterion at the actual values ( χ, f ). We have

Φ( χ + δχ, f + δf ) ĴΦ .δχ + Ĥ.δf = 0 ⇔ Σ f f = E (δf.δf T ) = Ĥ(-1) ĴΦ Σ χ χ ĴT Φ Ĥ(-1)
Thus, we can express (an approximation of) the covariance of the resulting transformation using the covariance on features and the criterion derivatives. However, a covariance matrix on a rigid transformation is quite hard to understand since it mixes angular values for the rotation and metric values for the translation. To characterize the transformation accuracy with a single number, we can compute the uncertainty (expected RMS error) induced on a set of representative points by the registration uncertainty alone (without the uncertainty due to the feature extraction). In our case, two sets are particularly well suited: the position of the matched extremal point represents the localization of the object of interest, whereas the corners of the image symbolize the worst case. In the example below, we find for instance a typical boundary precision around σ corn = 0.11 mm and a typical object precision far below the voxel size: σ obj = 0.05 mm for echo-1 registrations. The values are even a little smaller for echo-2 registrations: σ corn = 0.10 and σ obj = 0.045 mm.

Validation Index Last but not least, we need to validate this whole chain of estimations to verify if our uncertainty prediction is accurate. We observed that, under the Gaussian hypothesis, the Mahalanobis distance between the estimated and the exact transformation (or between two independent estimations of the same transformation) should be χ 2 6 distributed (if the covariance matrix on the estimation is exact). To verify this, the idea is to repeat a registration experiment N times and to compute the empirical mean value

I = μ2 = 1 N µ 2
i and the variance σ 2 I of this Mahalanobis distance. The values for an exact χ 2 6 are respectively 6 and 12. We can also verify using the Kolmogorov-Smirnov test (K-S test) that the empirical distribution corresponds to the exact distribution. The validation index I [START_REF] Pennec | A framework for uncertainty and validation of 3D registration methods based on points and frames[END_REF] can be interpreted as an indication of how the estimation method under-estimates (I > 6) or over-estimates (I < 6) the covariance matrix of the estimated transformation. It can also be seen as a sort of relative error on the error estimation.

We run several sets of tests with synthetic data and verify that our uncertainty estimations very perfectly validated for more than 15 extremal point matches. Now, the question we want to answer is: is it still valid for real data?

Validation with Real Data

The experiment is performed using multiple 2D contiguous Magnetic Resonance images (MRI) which constitute a 3D representation of the head. The images are part of an extensive study of the evolution of the Multiple Sclerosis (MS) disease performed at the Brigham and Woman's Hospital (Harvard Medical School, Boston) by Dr. Guttmann and Prof. Kikinis. Each patient underwent a complete head MR examination several times during one year (typically 24 different 3D acquisitions). The aim is to register precisely all the images acquired at multiple time points in order to segment the lesions and evaluate very accurately their evolution.

Each acquisition provides a first echo image and a second echo image (typically 256 x 256 x 54 voxels of size .9375 x .9375 x 3mm). The two images represent the same slice of T2 weighted signal imaged at different echo times. Thus, they are expected to be in the same coordinate system. This protocol was designed to optimize the contrast in the two channels for easier tissue segmentation. Considering two acquisitions A and B, the registrations of echo-1 images (A 1 to B 1 ) and echo-2 images (A 2 to B 2 ) give two relatively independent estimates of the genuine transformation from A to B. The comparison of these two transformations using the Mahalanobis distance gives a real validation index which can be tested for the accuracy of the uncertainty estimation.

In this experiment, the images being close enough, we used the iterative closest feature algorithm. Typically, we match 1000 extremal points out of the about 3000 extracted with a residual mean square error (RMS) of about 1mm.

Direct Validation Shows Biases

With n different acquisitions, we can run n * (n-1)/2 registrations per echo. In a first experiment, we compared directly the registrations between the corresponding echo-1 and echo-2 im-ages. The resulting validation index clearly indicates that the transformations do not agree (μ = I > 50 instead of 6). However, our registration method cannot detect systematic biases.

To discover the biases, we ran a series of experiments where we repeated the same registration while varying the algorithm parameters. This confirms that the observed uncertainty is similar in size and shape to the predicted one. Moreover, other experiments show that the inter-echo-1 and the inter-echo-2 registrations are compatible but the two groups significantly differ (Figure 13). Thus we concluded that there was a systematic bias between echo-1 and echo-2 registrations. Additional experiments showed that the bias was different for each registration. Estimation of the Biases To estimate the biases, we first observed that the transformation from image A 1 to image B 2 can be written

f A1B2 = f B •f AB1 = f AB2 •f A .
If measurements where perfect, the bias f A could be expressed for any other image Z:

f A = f (-1) AZ2 • f Z • f AZ1 .
Since measurements are noisy, we obtain an estimator of the bias f A by taking the Fréchet mean value [START_REF] Pennec | Computing the mean of geometric features -application to the mean rotation[END_REF]:

fA = arg min f    Z =A dist 2 f , f (-1) AZ2 • f Z • f AZ1    .
In this formula, each acquisition bias depends upon the others. Thus, we begin with null biases (identity transformations) and iteratively estimate each bias until convergence.

We effectively obtain a different bias for each acquisition that significantly differs from the identity. However, from a more global point of view, all the biases could be modeled as an "additional" noise on the transformation with an identity mean and standard deviations of σ r = 0.06 deg on the rotation (not significantly different from 0) and σ x = 0.09, σ y = 0.11 and σ z = 0.13 mm on the translation (significantly different from 0). Very similar values were observed for other patients.

Validation with Bias Although the biases appear very small, they are sufficient to explain the previous errors in the registration accuracy prediction. Indeed, taking the biases into account, the real validation index between acquisition A and B becomes:

I AB = µ 2 (f B • f AB1 , f AB2 • f A ) .
Since the biases are estimated from the registration values, using their uncertainties in this formula would bias the validation index toward low values. Thus we consider them as deterministic. The mean value and standard deviation of this new index across all registrations are now very close to their theoretical value (see table 1).

Origin of the Bias

Most of the extremal points we match are situated on the surface of the brain and the ventricles. These surfaces appear differently in echo-1 and echo-2 images due to the difference in contrast. Other artifacts such as chemical shift or susceptibility effects (see for instance [START_REF] Jezzard | Physical basis of spatial distortions in magnetic resonance images[END_REF]) may also account for the observed bias as they influence the detection of extremal points. Indeed, the two echoes are acquired with different receiver RF bandwidth to improve the signal/noise ratio [START_REF] Kikinis | Quantitative follow-up of patients with multiple sclerosis using MRI: Technical aspects[END_REF]. Therefore, the chemical shift and the susceptibility effect are different in the two echoes.

We plan to correlate the biases with diverse quantities in the images in order to understand their origin. Ultimately, we would like to predict the biases in each acquisition before registration. This would allow the definite validation of the registration accuracy prediction.

Conclusion

We presented in this chapter a method to extract reliable differential geometry features (crest lines and extremal points) from 3D images and several rigid registration algorithms to put into correspondence these features in two different images and to compute the rigid transformation between them. We also presented an analysis of the robustness with the computation of the probability (or mean number) of false positives and an analysis of the accuracy of the transformation.

This method proves to be very powerful for monomodal rigid registration of the same patient imaged at different times, as we show that an accuracy of less than a tenth of voxel can be achieved. In the last experiment, we showed that this uncertainty estimation technique is precise enough to put into evidence systematic biases on the order of 0.1 voxel between features in echo-1 and echo-2 images. Once corrected, multiple experiments on several patients show that our uncertainty prediction is validated on real data. This registration technique is currently used in many Chap. 31

Landmark-Based Registration Using Features Identified Through Differencial Geometry 512 medical applications, such as the registration of a preoperative MR used for planning and MR with a stereotactic frame for neuro-surgery navigation (European Community project Roboscope), or the registration of a series of acquisitions over time of images of Multiple Sclerosis patients to study the disease's evolution (European Community project Biomorph). Several tracks have been followed to generalize this work to non-rigid registration. Feldmar [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF] used the principal curvatures to register surfaces with rigid, affine, and locally affine transformations. Subsol designed an algorithm for non-rigid matching of crest lines. In [START_REF] Subsol | Crest lines for curve based warping[END_REF], he used this method to warp 3D MR images of different patients' brains in the same coordinate system, and even to warp an atlas onto a given 3D MR image of a brain in order to segment it. In [START_REF] Subsol | A general scheme for automatically building 3D morphometric anatomical atlases: application to a skull atlas[END_REF], he used the same method to construct automatically a morphometric atlas of the skull crest lines from several acquisitions of different patients' CT images, with applications in cranio-facial surgery.

Figure 1 :

 1 Figure 1: Differential geometry of 3D curves and surfaces. Left: principal directions and curvatures of a surface. Right: Frénet trihedron of a 3D curve and first differential invariants: curvature and torsion.

Figure 2 :

 2 Figure 2: Left: Crest lines as the intersection of two implicit surfaces. Right: Definition of the extremal points as the intersection of three implicit surfaces.

Figure 3 :

 3 Figure 3: Extraction of the Extremal Points. An empty circle denotes a positive value, whereas a filled circle indicates a negative one.

Figure 4 :

 4 Figure 4: Left: An axial slice of a 3D CT image of a dry skull in Plexiglas. Middle and right: the crest lines extracted on this image. The iso-intensity was manually chosen to delimit the skull. Original 3D Image courtesy of GE-CGR, Buc, France.

Figure 5 :

 5 Figure 5: Left: A sagittal slice of a 3D MR image. Middle and right: 3D views of the extremal lines extracted superimposed on the surface of the brain. Original MR images and segmentation courtesy of Prof. R. Kikinis, Brigham and Women's Hospital, Boston.

Figure 7 :Figure 8 :

 78 Figure 7: Example of registered crest lines between two CT skull images of the same phantom acquired in two different positions. Extremal points are represented by a color change from yellow to blue on the lines. Left: Front view with all crest lines from the two skulls after registration. Middle: Left view of the matched crest lines. Right: Closeup on the occipital foramen on the right. In this last image, the width of a line is a tenth of a voxel, which shows the very precise registration of these extremal points. One can also see that the trihedron part of the matched frames is very well conserved.

Figure 9 :

 9 Figure 9: Qualitative estimation of the number of false positives involving at least τ matches in MR images of 2500 features. Comparison between frames and points: we need roughly 5 times more point matches than frame matches to obtain the same probability (10 frames and 56 point matches for a probability of 10 -10 ).

Figure 10 :Figure 11 :

 1011 Figure 10: Graphical interpretation of the "compositive" noise model estimated on extremal points. The uncertainty of the origin (point X) is 4 times larger in the tangent plane than along the surface normal. The uncertainty of the normal is isotropic, whereas the principal directions t 1 and t 2 are 3 times more uncertain in the tangent plane.

Chap. 31 Figure 12 :

 3112 Figure 12: Left: Example of MS images. The same slice of one acquisition in echo-1 (left) and echo-2 (right). Right: evolution of an image row going through a lesion across 24 time points over a year. Left: without registration; Right: after registration and intensity correction. Original 3D Images courtesy of Dr. Charles Guttmann and Prof. Ron Kikinis from the Brigham and Woman's Hospital (Harvard Medical School, Boston).

Figure 13 :

 13 Figure 13: This diagram represents three acquisitions A, B, and C with the three echo-1 images (A1, B1, C1) and the three echo-2 images (A2, B2, C2). The echo-1 and echo-2 registrations are significantly different (µ 2 (fAB 1 , fAB 2 ), µ 2 (fAC 1 , fAC 2 ), µ 2 (fBC 1 , fBC 2 ) ¿ 50) but the intra-echo-1 and intra-echo-2 registrations are compatible (µ 2 (fBC 1 • fAB 1 , fAC 1 ) 6 and µ 2 (fBC 2 • fAB 2 , fAC 2 ) 6). This led us to assume a global bias for each acquisition between echo-1 and echo-2 images, represented here by the transformations fA, fB, and fC .

  Handbook of Medical Imaging -Processing and Analysis.I. Bankman Editor. Academic Press. 2000.

	Chapter 31
	Landmark-Based Registration
	Using Features Identified
	Through Differential Geometry
	Xavier Pennec
	Nicholas Ayache
	Jean-Philippe Thirion
	INRIA -Projet Epidaure

Table 1 :

 1 Theoretical and observed values of the real validation index with bias for different patients. The number of registrations (which is also the number of values used to compute the statistics on the validation index) is directly linked to the number of images used. Results indicate a very good validation of the registration accuracy prediction: the mean validation index is within 10% of its theoretical value and the K-S test exhibits impressively high values.

		Handbook of Medical Imaging -Processing and Analysis	IV Registration
	Theoretical values	Ī 6	√	σ I 12 = 3.46	K-S test 0.01 -1	num. im. n ≤ 24	num. reg. n * (n -1)/2
	patient 1	6.29		4.58	0.14	15	105
	patient 2	5.42		3.49	0.12	18	153
	patient 3	6.50		3.68	0.25	14	91
	patient 4	6.21		3.67	0.78	21	210

This draft paper is provided to ensure timely dissemination of scholarly and technical work. Copyright and all rights therein are retained by authors or by other copyright holders. All person copying this information are expected to adhere to the terms and constraints invoked by each author's copyright. This work may not be reposted without the explicit permission of the copyright holder.

Acknowledgments

The authors wish to thank the whole epidaure team and particularly their colleagues (or ex-colleagues) who worked on differential geometry features and registration: S. Benayoun, M. Fidrich, A. Gourdon, A. Guéziec, O. Monga, and G. Subsol. This research would have not been possible without the CT data from GE-CGR (Buc, France), the MR data of Dr. C. Guttmann and Prof. R. Kikinis (Brigham & Women's Hospital and Harvard Medical School, Boston, USA).


Marius George Linguraru

Michael Brady

An Anisotropic Diffusion Approach for Early Detection of Breast Cancer

Keywords: digital mammography, breast cancer, microcalcifications, h int image representation, shot noise, anisotropic diffusion

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction 1.1 The Impact of Breast Cancer

The incidence of cancer in the Western world is enormous. The impact that this fierce affliction has for a large percent of the population is a cultural phenomenon. Fear of cancer is an unfortunate reality that our auditive system has got used to! It has massive connotation for many domains, since the hospital lost its copyright for the term "cancer"' and now shares it with general sciences and even humanities, such as philosophy, literature and arts. Nevertheless, science is struggling to find new ways of improving worrying statistics, which show an increase in the impact that cancerous diseases are having especially over Europe and Americas.

Recent figures show that breast cancer is suffered by a high percent of the overall cancer incidence in women, where it accounts for approximately 25% of all female cancers. With about 600.000 new cases every year, breast cancer is second after lung cancer, the most feared form of cancerous death in women of all ages. The disease has a far higher incidence in Europe (especially Western Europe) and North America. In the far Orient and black Africa the mortality rate due to breast cancer is much lower, although there has been a substantial increase in the number of new cases. Such statistics represent intensive food for thought for researchers trying to trace genetic and environmental causes that lead to developing the disease.

Breast cancer's incidence is low in women under 30 years old and then increases with age. Between 40 and 50, women face a doubling of the rate of incidence which continues to increase over the age of 50, but more slowly. Younger women are encouraged to check the status of their breasts by simple palpation. Unfortunately, most women cannot reliably palpate a tumour smaller than 1 cm, therefore more examinations are required. In the UK, a screening programme is necessary for women over 50, since mammography is ineffective before menopause. There are several criteria that need to be fulfilled before starting a screening programme [START_REF] Evertsz | Soft Copy Environment for Screening Mammography -SCREEN[END_REF], such as:

• the disease to be screened must be very common and a treatment must be available for it, since there are very high costs involved and it would be of no use to look for a nontreatable affliction;

• the detection method must be robust, reliable and lead to good results for the overall screening process;

• must have high specificity;

• must be accepted by patients, since the method would not be cost-effective without a large number of patients to be examined. The screening programme mainly represents a mean of detecting cancer signs from an earlier stage and subsequently improving the prognosis for patients. If a mammogram presents any features that seem suspicious to the radiologist, the patient will be asked to attend an assessment clinic where more investigations will be performed by means of medical imaging and consulting. Figures show that 8% of women are recalled for further investigations, most of them not presenting any malignancy. The screening programme should, according to optimistic statements, almost double the chances of survival in women developing breast cancer. There is intensive debate as to whether the breast screening clinical assessments should be performed every two years instead of the three year period presently considered acceptable, as it has been noticed that the assessed women sometimes develop cancer over a period shorter than three years ("interval cancer").

When speaking about breast cancer, one generally refers to a set of breast afflictions, children of the same manyheaded monster, which have one common wicked status: the presence of malignancy. Moreover, we are dealing with different stages of the same disease and the treatment given to affected women differs considerably from one stage to another. So does the prognosis, which represents a statistical measure of the chances of having a positive outcome of the treatment the patient is undergoing. The most relevant factors a doctor considers when deliberating on the prognosis are the size and type of the tumour, the presence of metastases, the stage of the disease, the status of the axillary lymph nodes and patient's age and medical condition.

A Brief Anatomy of the Breast

For a better understanding of the subject, we first overview the breast anatomy. The female breast has two separate functional parts:

• the milk-producing tissue (the epithelial tissue) which appears very dense in mammograms due to the high percent of calcium it contains;

• the other types of breast tissue, including water-dense non-fat tissue (which also appears light, but not as bright as calcium), breast fat-tissue (which is transparent in X-ray imaging and corresponds to the milk-producing tissue converted to fat after menopause), fascia (the connective tissue), and muscles. Figure 1 shows a comparison between the appearance of the mammography of a breast of a young woman and an image of a breast of a post-menopausal woman.

The routes preferred by cancer in its spreading process are the blood and lymph vessels (which make the axillary nodes an important feature in signalling metastatic diseases), but the direct invasion of the surrounding tissue may have the same effect. Therefore, one would naturally speak about spreading (invasive) and non-spreading (in-situ) forms of cancer.

The Need for Image Segmentation

Image processing is a challenging but difficult task. Working with mammograms is especially complex due to the complicated appearance of the structures of interest in this particular type of image representation. Although a mammogram is a good picture of the breast, this is hardly sufficient when searching for small or complex anatomical parts, such as microcalcifications, masses or curvilinear structures in the process of early detection of breast cancer.

Statistics show that approximately 25% of cancers are missed and about 80% of the biopsies are performed for benign cases.

The complicated anatomy of the breast is an inevitable source of the highly textured structure of the mammograms. It provides a most difficult to analyse input for the radiologist, who is expected to distinguish very subtle abnormalities out of this mass of structural ambiguity. The variability between any two different cases adds to the difficult task that the human decision maker faces.

Image Quality

The quality of the image depends on several physical factors.

The most important are the different breast deformation during the compression process, the time of exposure, the breast thickness and various other imaging factors. Since the X-ray dosage cannot be too high, for patient safety reasons, there is a compromise between dosage and the signal-to-noise ratio (SNR) of the mammogram. This is even more degrading to the quality of the image, which must reflect the superimposed structures of the breast (a 3D structure) on a 2D projection.

Noise is added to these features. This presentation of the mammogram could be easily reduced to that of a large textured noisy image, which still represents the best tool for early cancer detection to be used around the world. Section 2 will present some major improvements in image acquisition developed by the Medical Vision Laboratory at the University of Oxford.

Digital Mammography

Mammography is predominantly film-based and represents the outcome of three subsequent processes, image acquisition, image storage and viewing. The overall process is dependent on the same medium, the film. Although many technical problems have retarded the digitisation of mammography, and still prevent it from being used widely, directly digital mammography is now available [START_REF] Evertsz | Soft Copy Environment for Screening Mammography -SCREEN[END_REF]. The digital systems make a clear distinction between image acquisition, image storage and image visualisation. There are several clear advantages that digital imaging has to offer:

• digital image acquisition is expected to reduce the X-ray dose in the imaging process, leading to less risk for the patient; it also improves the contrast resolution; • digital image storage reduces the cost of the operation since it does not require film and chemicals, provides fast and reliable retrieval from the archive and allows the collection of large image databases;

• digital image visualisation enables the use of Computer Aided Diagnosis (CAD) for automatic detection, data documentation and the universe of image processing; it greatly improves the quality of the displayed image over that of the rigid film. Early detection can be accomplished only with high quality detection and processing tools in the different levels of the screening programmes. Directly digital mammography is the new trend in X-ray mammography and will bring serious improvements in the development of screening programmes. The soft-copy reading environment [START_REF] Evertsz | Soft Copy Environment for Screening Mammography -SCREEN[END_REF], a computer workstation which displays digital mammograms for the radiologist to read, is the tool that will make digital mammography available in hospitals for much better results in detecting the structures of interest within the breast.

Future Trends

There are several implemented techniques, which bring improvements in the field of medical imaging. Most of them have proved not to be satisfactory for the purpose that they are used, such as automatic detection and working with mammograms. Moreover, although the resulting images they provide look quite impressive, it has not been definitely proved that radiologists work better on these images. Therefore, there is sufficient room for improvements and further developments in image processing. Some present trends in the field include:

• the development of soft-copy reading workstations [START_REF] Evertsz | Soft Copy Environment for Screening Mammography -SCREEN[END_REF], as the tool for the future use of digital mammograms;

• friendlier user-interfaces (touch-screen, automatic report generation, robust display) which would only require a minimal intervention from the human factor involved, the radiologist;

• the development of training-systems with immediate feedback and the use of larger databases;

• the development of real-time applications for making the best use of the image processing methods in clinical applications;

• software integration of the existing algorithms is an inherent condition in building strong performant medical systems;

• finding more reliable, more robust and faster image processing algorithm that will have a great impact on the future of medical analysis.

H int Model

The h int representation [START_REF] Highnam | Mammographic Image Analysis[END_REF] is a physical-based approach to mammographic analysis, an image normalisation method based on a complete understanding of the imaging process. Since the quality of mammograms is so highly dependent on the imaging conditions, the h int model is an alternative quantitative representation of the breast tissue. Figure 2 shows a depiction of the h int surface of a breast.

The intensity or attenuation value of a pixel in a mammogram is determined by the quantity of X-ray beam absorbed by the tissue present between the X-ray source and the respective pixel. The h int value of a pixel represents the thickness of the breast tissue of interest underlying between the X-ray source and the actual pixel. By interesting tissue one must consider the non-fatty tissue present in the breast, such as glandular, cancerous and fibrous tissue, which have high attenuation. Hence, the h int representation is not dependent on the imaging procedure the same way the intensity value is. Other types of tissue present within the breast are the fatty tissue with low attenuation and the calcifications with very high attenuation, since they contain a concentrated level of calcium. They lead to the definition of h fat and h calc , the thickness of the fatty tissue and calcification. Calcifications are very small anatomical features of the breast, so we can compute the total thickness of the compressed breast with:

H = h int + h fat (1)

Figure 2. A depiction of the h int surface of a breast

Generating H int

Although formula (1) has a very simple appearance, the computation of h int is a rather complex process. Here are the steps to be followed to generate h int [START_REF] Highnam | Mammographic Image Analysis[END_REF]:

• convert the pixel value P(x,y) into film density D(x,y) and deblur; this is achieved by considering the linear relationship between P(x,y) and D(x,y) followed by the use of the modulation transfer function; • convert the film density D(x,y) into energy imparted to the intensifying screen E pse imp (x,y); the film-screen calibration data is relevant for this purpose; the energy imparted image appears as an inverted version of the density image, where dark parts correspond to high energies;

• compensate E pse imp (x,y) for the intensifying screen glare by using the point-spread function of the intensifying screen;

• compensate E pse imp (x,y) for the anode-heel effect and diverging X-ray beam, by knowing the variation between the incident photon flux between different spatial locations on the film;

• estimate the scattered radiation E s imp (x,y), since this component of the imparted energy contains no information about the breast tissue but influences the neighbourhood of the pixel;

• estimate the extra-focal radiation E e imp (x,y) component, which is relevant at the curved breast edge where photons arriving with low angles can reach the intensifying screen;

• compute the primary radiation E p imp (x,y);

E p imp (x,y) = E pse imp (x,y) -E e imp (x,y) -E s imp (x,y) (2)
• convert E p imp (x,y) into h int using the conversion equations developed by Highnam and Brady [START_REF] Highnam | Mammographic Image Analysis[END_REF].

Strengths of the H int Model

The h int representation is a robust and reliable method, resulting in a floating point form, which corresponds to the thickness of interesting tissue in the breast. By removing most of the unwanted effects of the imaging process, such as glare, scatter radiation, anode-heel effect and extra-focal radiation, the output of the method presents a much more adequate representation of the real anatomical structure of the breast. Hence, the h int representation of the breast is mainly a 3D surface built from the h int values of the image pixels. By removing the image parameters, the h int images stand as normalised images of the breast. image.

An h int representation can be easily visualized as an image, since the h int values are in float format, where brighter parts correspond to thicker parts of the breast or calcifications as in Figure 3. Moreover, the depicted surface of an h int representation of the breast can show important anatomical features, such as masses as hills in a less dense background, while the background is mainly flat.

Abnormal Structures in Breast Anatomy

classification of tmours are open subjects for future research. Section 3.2 reviews the state of the art in the detection and classification of masses in X-ray mammography.

Detection of Microcalcification

Microcalcifications represent some of the earliest signs of breast cancer and they account for half of the non-palpable lesions that appear in mammograms. Calcifications are either warnings of malignancy or just benign signs of aging. They are encountered in approximately 25% of mammograms and appear as bright spots or clusters of such spots, due to the high density of calcium. To differentiate malignant from benign microcalcifications, radiologists believe that they take into consideration several criteria, such as shape, size and arrangement.

According to their size, calcifications can be classified into macrocalcification or coarse calcifications (when their size usually exceeds 1 mm) and microcalcifications. While it is believed that large single regular-shaped calcifications are benign, small clustered whorled calcifications are more likely to signify malignancy, as noticed by Caseldine et al} [START_REF] Caseldine | Breast Disease for Radiographers[END_REF] and Le Gal et al. [START_REF] Gal | Mammographic Features of 455 Invasive Carcinomas[END_REF]. Although the detection and classification of calcifications are two fields that have improved significantly in recent years, there is still no robust differentiation between benign and malignant calcifications and a large enough number of calcifications are still not detected by radiologists or by detection algorithms. Figure 4 The aim of automatic detection is to find a sufficiently reliable algorithm to be used in clinical practice. Such algorithms are meant to assist the radiologist in making decisions and to improve the overall sensitivity (the proportion of true positives, TP) and specificity (1-the proportion of false positives, FP) of the detection process. Some of the factors that drastically influence the TP and FP figures are:

• the variability of the anatomy of the breast; every mammogram has different features related to different tissue types and correspondingly variable brightness in the mammographic appearance;

• the imaging conditions, such as shot noise, patient movement, low contrast in mammograms, due to low X-ray dosage and glare;

• faint microcalcifications lost in a dense background, the superposition of certain breast structures. In order to overcome the above-mentioned compromising factors, most of the conventional detecting algorithms comprise three main stages:

• a pre-processing step based on filtering the image; the filter is meant to detect and remove the noise and enhance the structures of interest, the microcalcifications; Chan et al} [START_REF] Chan | Image Features Analysis and Computer-Aided Diagnosis in Digital Radiography: 1.Automated Detection of Microcalcifications in Mammograms[END_REF] and Nishikawa et al} [START_REF] Nishikawa | Computer-Aided Detection of Clustered Microcalcifications: an Improved Method for Grouping Detected Signals[END_REF] • a clustering step using a fixed size kernel to eliminate noise points and isolated calcifications and identify clusters. Although most of the conventional detection methods are based on the above algorithm, the literature includes a couple of novel detection methods. Yam et al. [START_REF] Yam | Denoising h int Surfaces: a Physics Based Approach[END_REF], [START_REF] Yam | Detecting calcifications in mammograms using the h int representation[END_REF] present a physics-based approach, which uses both grey level and h int images. The method extracts regions corresponding to the size of microcalcifications that satisfy a criterion related to the change of surface with the height. The ratio between the actual volume of the extracted blob and the volume of the interesting tissue (h int ) is the second feature that is introduced.

Another innovative detection method is Karssemeijer's statistical approach [START_REF] Karssemeijer | Adaptive Noise Equalisation and Recognition of Microcalcification Clusters in Mammograms[END_REF], [START_REF] Karssemeijer | Stochastic Model for Automated Detection of Calcifications in Digital Mammograms[END_REF]. An adaptive noise equalisation algorithm was developed to deal with the variation of noise characteristics in an image. It makes use of Bayesian techniques combined with Markov random fields.

Detection of Masses

Masses represent a special category of breast abnormalities and their detection is a subject that is intensely debated and researched in digital mammography. Generally, masses appear as bright regions in a mammogram with rather well contoured boundaries. Nevertheless, it is particularly difficult to identify some mass-like structure, since the interesting dense tissue within the breast may have similar density values and overlap with masses. In the above-mentioned case, it is a real challenge to identify the shape and the edges of the mass. An important source of FPs is the overlapping of two dense regions, due to the projective nature of a mammogram, which would look brighter that its surroundings and be confused with a mass. For this reason, two views are taken from the breast, a cranio-caudal (CC) and a medio-lateral-oblique (MLO) view. Figure 5 shows two examples of samples of mammograms containing masses. a.

b.

Figure 5. (a) A mammogram sample showing a benign mass; (b) a malignant mass.

When related to masses, automated image analysis relies on three main applications:

• locating abnormal regions in a single mammogram, when, as a result of several features calculated, possibly pathological regions of the breast are extracted and pointed out to the radiologist; • matching bilateral pairs, making use of the same-view left and right breast mammograms of the same woman at the same time, acknowledging the approximate symmetry of the two breasts; • matching temporal pairs, using the same-view of the same breast mammograms of the same woman at two different times, searching the major changes that appear between the two mammograms; [START_REF] Marias | Multi-Scale Landmark Selection for Improved Registration of temporal Mammograms[END_REF] Amongst the most relevant achievements in detecting breast lesions, te Brake and Karssemeijer [START_REF] Brake | Single and Multiscale Detection of Masses in Digital Mammograms[END_REF], [START_REF] Karssemeijer | Detection of Stellate Distortions in Mammograms[END_REF] consider the irregular shape that so-called speculated masses may have and the often noted stellate appearance of them, as radiating patterns of linear spicules. This method uses gradient operators for determining the relevant line orientation. Two features are considered: the total number of pixels with an intensity gradient pointing to a centre and the likelihood of having found a speculated structure, which becomes higher with the increasing number of directions in which pixels are oriented towards a region.

Chan et al. [START_REF] Chan | Improvements of Radiologists' Characterization of Mammographic Masses by Using Computer-aided Diagnosis: An ROC Study[END_REF], [START_REF] Sahiner | Computerized Characterization of Masses on Mammograms: The Rubber BAnd Straightening TRansform and Texture Analysis[END_REF] use a rubber-band-straightening transform (RBST) to transform the band of pixels surrounding a mass into a Cartesian plane (the RBST image), which has the advantage of making the mass margins almost parallel, with spicules to perpendicular the length of the rectangle. A large number of features are computed and the most relevant of them for used for Fisher's linear discriminant that will conclude with the likelihood of malignancy.

Highnam and Brady [START_REF] Gal | Mammographic Features of 455 Invasive Carcinomas[END_REF] comment on the advantages of using h int representations in the way that a mass would correspond to a hill-like structure surrounded by a smoother region corresponding to fibroglandular tissue. The ultimate aim of the detection of masses is their interpretation; therefore, the shape and arrangement of the salient region are relevant features in this process.

A Review of Anisotropic Diffusion

It has been often proved that filtering methods can substantially improve the quality of the image of interest by means of extinguishing artefacts and unwanted information from the original image, but also simplifying the appearance of otherwise complicated anatomical structures. The notion of image filtering has already been encountered in the previous section as an image pre-processing step used by some of the overviewed algorithms. Conventional filtering methods include a background smoothing stage (convolution with a low pass filter) followed by enhancement of the structures of interest (high pass filtering) and the subtraction of the two newly obtain images. This method could not deal with the large variability of the anatomical features that must be considered.

There are numerous algorithms based on evolving partial differential equations (PDE) for noise removal and image enhancement, but few of them have been previously tested on medical imaging. Furthermore, the extent of applying PDE filters to medical images seems to be generally related to ultrasound and MRI.

Anisotropic diffusion has its origins on the classical nonlinear diffusion filter developed by Perona and Malik in 1987 [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], which is based on a PDE in divergence form. It is the cornerstone for new developments in multi-scale image analysis aiming to simplify the image appearance while enhancing structures of interest, such as edges or coherent structures.

The Diffusion Process

One of the most commonly used methods for smoothing an image f : R 2 → R is by convolving it with a Gaussian with standard deviation s:

K s (x) = (1/2p s 2 ) • exp (-|x| 2 /2 s 2 ) (3)
The image f is transformed into a family of gradually smoother versions for a higher number of iterations (t>0). An increasing scale will simplify the appearance of the original image. There are several limitations of this method, as observed in [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF]:

• although convolution with a Gaussian reduces noise, it also blurs important anatomical structures of the images, such as edges; • linear diffusion dislocates edges when changing from finer to coarser scales.

The diffusion process was proposed as an alternative to smoothing images by a Gaussian kernel, which does not preserve edges. Since it derives from a process of equilibrating concentration differences, it can be expressed through a continuity equation of Fick's law [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF]:

j = -D •∇u (4) δ t = div(D •∇u) (5)
D is called the diffusion tensor, a positive definite symmetric matrix which represents the relation between the concentration gradient ∇u and the flux which aims to compensate this gradient [START_REF] Catté | Image Selective Smoothing and Edge Detection by Nonlinear Diffusion[END_REF]. In image processing, the concept of concentration is replaced by that of grey level. The diffusion tensor may be replaced by a positive scalar-valued diffusivity g. If j and ∇u are parallel, the diffusion is called isotropic. In the anisotropic case, j and ∇u are not parallel. Equation ( 5) is called the diffusion equation. If the diffusion tensor is space-dependent, then the diffusion is called inhomogeneous, while a constant diffusion tensor is related to a homogeneous diffusion.

In order to overcome the correspondence problem (the coarse-to-fine tracking difficulties), the inhomogeneous linear diffusion filtering introduces |∇ f| as edge detector, where high values of the detector indicate the presence of edges in the image. The diffusivity function g was set to:

g(|∇ f| 2 ) = 1/sqrt(1+|∇ f| 2 /λ 2 ) ( 6
)
and the diffusion equation reduces to:

δ t u = div(g(|∇f| 2 ) ∇u) (7)
Introducing a feedback in the diffusion process, by adapting the diffusivity to the gradient of u(x,t) --the actual image, rather than the original f(x) --the diffusion equation becomes nonlinear and therefore the diffusion filtering will be nonlinear and isotropic:

δ t u = div(g(|∇u| 2 ) ∇u) ( 8
)
The Perona-Malik Model [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] represents the first nonlinear diffusion filter. It provides stable edges over a long number of iterations based on a rapidly decreasing diffusivity, but will only enhance those edges for which the gradient is larger than the contrast parameter ?.

g(|∇ f| 2 ) = 1/(1+|∇u| 2 /λ 2 ) (? > 0) (9)
Catté et al. [START_REF] Catté | Image Selective Smoothing and Edge Detection by Nonlinear Diffusion[END_REF] introduced the Gaussian convolution of u: u s =K s *u and the result of it was: [START_REF] Highnam | Mammographic Image Analysis[END_REF] This new form of the diffusion equation solved the spatial regularizing problem of the inhomogeneous filtering, meaning that the solution of the nonlinear filtering method of images aims to a steady state. Moreover, a new parameter is introduced in the process, the scale parameter s. The process is now controlled by the three parameters, t(time), ?(contrast) and s(scale), which tremendously reduce the impact of the choice of diffusivity over the whole process and make the use of it more flexible and robust. Although the contrast parameter will work similar to the Perona-Malik model, the scale parameter will make the filter less sensitive to smallsize structures, such as noise, by increasing the kernel of the Gaussian s.

δ t u = div(g(|∇u s | 2 ) ∇u)

Nonlinear Anisotropic Diffusion

The main improvement introduced by the nonlinear anisotropic filters is the smoothing along the isophote and not across it when the value of the gradient is large [START_REF] Kornprobst | Nonlinear Operators in Image Restoration[END_REF]. While for low gradients smoothing is performed in the usual way, at edges diffusion is inhibited. Weickert [START_REF] Weickert | Theoretical Foundations of Anisotropic Diffusion in Image Processing[END_REF], [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF] introduces a system of eigenvectors v1, v2 of the diffusion tensor D. v1 and v2 are orthonormal and

v1 ¦ ∇u s , ( 11
)
v2 ⊥∇u s ( 12
)
The corresponding eigenvalues are:

λ 1 = g( |∇u s | 2 ), (13)
λ 2 = 1. (14)
In general ∇u s is not parallel to one of the eigenvectors of D for σ>0 and Weickert's model behaves really anisotropically. If σ tends to 0, the process tends to behave like the classical Perona-Malik model.

The edge-enhancing diffusion model proposed by Weickert is the one used in our initial experiments, where the diffusion across edges is performed according to the following eigenvalue:

 1 |∇u s | = 0, λ 1 =   1 -exp(-3.31488/(|∇u s |/λ) 8 ) |∇u s | > 0. ( 15
)
An example showing different stages of applying anisotropic diffusion to a sample of mammogram displaying a large microcalcification and some noise is shown in Figure 6.

Discussion

Anisotropic diffusion overcomes some major limitations of linear and nonlinear isotropic filters:

• enhances noisy edges and flow-like structures (would be useful in the detection of curvilinear structures);

• inhibits diffusion at edges; • more flexible due to the larger number of parameters, but not large enough to alter the robustness and accuracy of this method;

• converges to a steady state for t → ∞ respecting a maximum-minimum principle; • its solution is unique and continuously dependent on the initial image. The use of anisotropic diffusion, as observed in [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF], spreads from computer aided quality control to postprocessing fluctuating data, target tracking in infrared images and blind image restoration, to only enumerate some applications. But most applications are concerned with filtering medical images, mainly ultrasound and magnetic resonance images.

It is however not sufficient to use ``blind'' filtering methods when the features we need to preserve in an image are so precise and specific. The use of it a priori or even a posteriori knowledge in the diffusion process must be a future trend in the use of PDE-based filters, a field with continuous evolution which has to bring many advantages in the overall development of medical vision.

Filter Model

One of the major characteristics that we use in approaching microcalcification detection is the genuine difference that should be visible in the shape of microcalcification versus noise in mammographic imaging. While microcalcifications are anatomical structures with slightly blurred edges, which appear in mammograms due to the effect of X-ray beams passing through the breast anatomical structure, noise tends to have extremely sharp edges, as in Figure 7. The appearance of h int images would be extremely noisy mainly due to the removal of the glare effect, extra-focal and scattered radiation (which accounts for up to 40% of the total radiation exiting the breast). Although calcifications are small structures and they may be computationally omitted, the high attenuation of calcium makes them stand as a definite exception in the h int representation. Since microcalcifications are high in h int images, only the most prominent spots of noise may lead towards FPs, the smaller ones being easily removed by the diffusion process.

We chose an anisotropic diffusion-based filter, which attempts to blur the input mammographic image while preserving significant intensity changes. The process relies on the use of a set of different parameters, e.g. time ( t), contrast (k), size (s), and it is critical to find the right choice of parameters that will lead to good repeatable results. The nonlinear anisotropic diffusion technique proves to be highly flexible de to the variability of its parameters which help in covering a rather extensive set o possibilities with respect to the output one can get by filtering medical images, as Table 1 shows. The diffusion tensor for the anisotropic filtering is based on Weickert's approach having the corresponding eigenvalues, λ 1 and λ 2 in ( 16), where u s represents the Gaussian convolution of the original image u.

Diffusion and Surface Plots

We show some results from applying non-linear anisotropic diffusion filtering on samples of h int representations of real mammograms containing microcalcifications in Figure 8. They demonstrate the accurate effect of de-noising h int images performed by our method while preserving only calcifications and significant noise points. We also considered images with very high likelihood to present false positives (FPs). Such an example is presented in Figure 9.

First Results

The detection method, for both calcifications and noise was based on the association one can make between the original h int mammograms containing the structures of interest and the surface we built from the filtered images after just a few iteration steps. The surface we present would show either hillshaped structures for microcalcifications or sharp-edged formations for noise in the locations corresponding to the structures of interest. The algorithm was tested initially on a set of 13 samples of average 32-bit h int mammograms at 50 mm resolution containing 10 pre-labelled coarse calcifications and several artefacts. Samples were preferred, rather than whole mammograms, in order to reduce processing time. All the calcifications were correctly detected. No FPs were detected during our initial experiments. The free-response receiver operating characteristic (FROC) curve is shown in Figure 10a.

The algorithm was further tested on a set of 20 samples of 32-bit h int mammograms at a resolution of 50 mm containing 27 pre-labelled isolated microcalcifications and various bits of noise. The TPs fraction was 92.6% for a number of 0.

A Non-Parametric Approach

We designed an adaptive Gaussian derivative filter to apply to the de-noised, glare-removed h int images, which results in a gradient map (Figure 11.b) that highlights suspicious regions as regards to microcalcification detection. The same filter outputs the value of k for the subsequent diffusion. Having the gradient map and the value of k, we apply the anisotropic diffusion filter to the h int image over a constant number of iterations t, with a pre-defined value of s and the corresponding value of k. The new image will generally be blurred (Figure 11.c), excepting some suspicious regions which will have their edges preserved. With the diffused images ready, we need to ensure that the artefacts emphasised in the gradient maps will be eliminated. A new filtering step is introduced at this stage to actually depict the microcalcifications in the processed image. It is based on the same computation as the anisotropic diffusion, but also incorporates some adaptive thresholding suited to the image characteristics and the properties of microcalcifications. This last filter will conclude with a black-and-white-map-of detection (BWMD -Figure 11.d).

In summary, our method uses three filters in sequence: • an adaptive Gaussian filter;

• an anisotropic diffusion filter;

• an adaptive thresholding filter.

The method eliminates the initial parametric dependence and improves the overall robustness of the detection. We The overall detection rate for microcalcifications became 91.3% TPs for a number of 0.2 FPs per image. Although the detection rate is slightly lower than in some of our previous experiments, this is the price to be paid for improving robustness.

Further Work

• Breast Edge Correction: the high gradient value around the boundary of the breast is a major source of errors; the solution would be bringing the background to an h int value similar to that of the average breast tissue. • Shot Noise Removal: shot noise is the major source of errors in the detection of microcalcifications; the h int generation can detect shot noise and it can be therefore extracted by interpolation from the images we use, reducing the number of FPs.

• Wavelet Implementation: anisotropic diffusion is not sensitive to extremely small structures, therefore, by using wavelets we should improve the performance of the algorithm in dealing with very small salts of calcium; a combination of anisotropic diffusion and wavelets could be the key to our problem. While our current method could detect the general shape of a cluster and most of its calcifications; further developments should search within the boundaries of the cluster for particularly small and faint microcalcifications.

Figure 1 .

 1 Figure 1. (a) A pair of mammograms taken from a young woman showing mainly dense tissue in the breast; (b) a pair of mammograms taken from a post-menopausal woman where the milk-producing tissue has converted to fat.

Figure 3 .

 3 Figure 3. (a) An original mammogram and (b) its corresponding h int .image.

  shows two examples of samples of mammograms containing calcifications. a.b.

Figure 4 .

 4 Figure 4. (a) A mammogram sample showing two isolated coarse calcifications; (b) a microcalcification cluster.

Figure 6 .

 6 Figure 6. (a) The original sample of a mammogram with a large calcification on the upper-right side; (b) the diffused image with k=5, s=0.6 and t=20; (c) the diffused image with k=5, s=0.5 and t=20, where only small structures are kept and their edges enhanced; (d) further diffusion the image of (c) k=5,s=0.5 and t=500, the large bit of noise if removed.

Figure 7 .

 7 Figure 7. (a) The plot of a filtered intensity image containing a microcalcification at position 150; (b) the plot of a filtered image containing noise at position 147.

Figure 8 .

 8 Figure 8. (a) A pre-processed h int image containing a microcalcification on the left side and a large spot of noise on the lower right side; (b) the diffused version of the (a) image with k=15, s=0.6 and t=5; (c) the surface of the original h int image in (a)with highly noisy appearance; (d) the surface of the diffused h int image, the microcalcification appears as a hill with smoother edges than those of noise.

Figure 9 .

 9 Figure 9. (a) The original pre-processed h int image containing only noise structures; (b) the diffused version of the (a) image with k=15, s=0.6 and t=3; (c) the surface of the original h int image in (a); (d) the surface of the diffused h int image, where all structures have very sharp edges and are labelled as noise.

1

 1 FPs per image. We further applied an implementation of Yam et al.'s algorithm on the filtered images of the same set of microcalcifications. We obtained a 100% fraction of TPs with a number of 0.3 FPs per image. The FROC curve of the detection using the combination of the anisotropic diffusion filter and the algorithm implemented by Yam et al. (1999) is shown in Figure 10b. a.b.

Figure 10 .

 10 Figure 10 . (a) The FROC curve for the set of coarse calcifications; (b) the FROC curve of the combined detection method for the set containing different types of microcalcifications.

  show three examples of detection: Figure11-the detection of a Ductal Carcinoma In Situ; Figure12-the detection of subtle microcalcifications; Figure13-avoiding FPs in detection by eliminating curvilinear structures (CLS)[START_REF] Evans | Detecting and Removing Curvilinear Structures from Mammograms[END_REF],[START_REF] Kovesi | Image Features from Phase Congruency[END_REF].

Figure 11 .

 11 Figure 11. (a) An h int sample containing a Ductal Carcinoma In Situ; (b) the gradient map of the (a) image with k=3; (c) the diffused image of (a); (d) the BWMD where the white regions correspond to calcium.

Figure 12 .

 12 Figure 12. Shows a difficult case of detection: (a) an h int sample containing two subtle calcifications; (b) the corresponding gradient map; (c) the diffused image of (a) after using the computed value of k; (d) the BWMD with the depiction of the two microcalcifications.

Figure 13 .

 13 Figure 13. (a) An h int sample containing one subtle microcalcification and a prominent CLS; (b) the BWMD of (a) showing an FP; (c) the same image as (a) after CLS removal and contrast enhancement; (d) the BWMD of (b) with the correct detection of the microcalcification.

  use a difference image, the result of subtracting a signal-enhanced image by a signalsuppressed image; • a segmentation step based on adaptive thresholding or local contrast; Nishikawa et al. use both global thresholding and morphological erosion;

Table 1

 1

		Blur	Anatomical features	Edges
	kä	ä	ae	ae ae
	σä	ä	ae	ae
	tä	ä	ae ae	Well preserved

This section aims to overview briefly some of the most relevant papers that relate to the subject of this research. Section 3.1 presents the most important achievements to date in detecting microcalcifications. The detection and