Gilles Bertrand

Hugues Talbot

Petr Dokladal

Isabelle Bloch

Yoshinobu Sato

Maximum Intensity Projection DT Distance Transform GVF Gradient Vector Flow SNR Signal-to-Noise Ratio SE Structuring Element SV Spatially-Variant SVMM Spatially-variant Mathematical Morphology GT Ground Truth CFA Color Filter Array TP True Positive TPR True Positive Rate FP False Positive FPR False Positive Rate ROC Receiver Operating Characteristics

Keywords:

To those who were there and who are always here.

A C K N O W L E D G M E N T S

During my Master thesis studies specialising in Image Analysis at the University of Uppsala, Sweden, I got interested in the French culture. A simple google search on the image analysis laboratories in France with a sequence of lucky coincidences have led me to three and a half years of research on the inspiring medical imaging applications and many other changes in my personal life.

In addition to getting closer to knowing the French culture and language during this time, I have been lucky to met many wonderful people, some staying in my life as love and dear friends, others staying as prints in my heart and memory.

I would like to thank my thesis referees Isabelle Bloch and Yoshinobu Sato for accepting to sit in the jury of this thesis. I am honored by this and grateful for their time and effort. I especially would like to thank Professor Sato for travelling such a long distance.

I am would like to thank Gilles Bertrand, the director of the laboratory A3SI for his welcoming and generous help. Moreover, I am grateful for the environment present in the laboratory that unites all the members and makes the work in here so pleasant.

The author is also very grateful to Fernand Meyer for all his efforts in the collaboration between the University Paris-Est and Ecole de Mines de Paris, Center of Mathematical Morphology and his kind help, especially during the final six months of the thesis.

I am endlessly grateful to my thesis supervisor Hugues Talbot for the many things he has done and said to me since we came to know each other. In particular and in the first place, it is mainly thanks to his time and effort in various bureaucratic battles that this thesis took place. Apart from the supervising, code-writing, and idea-crisis help, I am grateful for his friendship and kind support.

I am also very grateful to Petr Dokladal, my second supervisor, for all the help and effort he has put into this joint work. Moreover, I am grateful for the personal lessons I have learnt, maybe without showing my appreciation enough.

The author would like to thank Nicolas Passat, Michel Couprie and John Chaussard for their contribution with materials, methods and generous insights.

I would like as well to thank the medical personnel from Colmar Hospital, especially Mariano Musacchio and Michel Lagneau for their interest in our ideas, kind provision in medical expertise and data.

I would especially like to thank Vincent Bismuth for the provided expertise in guidewire images, materials and image data that were used in this thesis.

The author is also very grateful to all the people of the two respectful laboratoies for their help, support, collaboration and kindness without whom this work would not have been so pleasurable.

Notably, I am thankful to Yukiko Kenmochi for her friendly listening, understanding and shared times. v I am grateful to the "comrades" students for their friendship and help during all these years, notably: Yohan Thibault, John Chaussard, Roland Levillain, Ania Jezierska, Jean Cousty, Benjamin Raynal, Camille Couprie, Harold Phelippeau and others with whom we shared the fun time.

I would like to thank my friends who were there and who are still there, notably Virginie and Johann.

I have been lucky to be a daughter of my parents that have supported me on this long, steep education road full of sacrifices. I am endlessly grateful for their efforts and sacrifices that they have done and continuously do for me.

Last but not least, I am grateful to my companion, Laurent, for being the most believing, optimistic and encouraging. Even in harder times, it has been so much better with you.

vi A B S T R A C T

The motivation of this work is filtering of elongated curvilinear objects in digital images. Their narrowness presents difficulties for their detection. In addition, they are prone to disconnections due to noise, image acquisition artefacts and occlusions by other objects. This work is focused on thin objects detection and linkage. For these purposes, a hybrid second-order derivative-based and morphological linear filtering method is proposed within the framework of scale-space theory. The theory of spatially-variant morphological filters is discussed and efficient algorithms are presented.

From the application point of view, our work is motivated by the diagnosis, treatment planning and follow-up of vascular diseases. The first application is aimed at the assessment of arteriovenous malformations (AVM) of cerebral vasculature. The small size and the complexity of the vascular structures, coupled to noise, image acquisition artefacts, and blood signal heterogeneity make the analysis of such data a challenging task. This work is focused on cerebral angiographic image enhancement, segmentation and vascular network analysis with the final purpose to further assist the study of cerebral AVM.

The second medical application concerns the processing of low dose X-ray images used in interventional radiology therapies observing insertion of guide-wires in the vascular system of patients. Such procedures are used in aneurysm treatment, tumour embolization and other clinical procedures. Due to low signal-to-noise ratio of such data, guide-wire detection is needed for their visualization and reconstruction. Here, we compare the performance of several line detection algorithms. The purpose of this work is to select a few of the most promising line detection methods for this medical application.

keywords: angiography, directional filtering, multi-scale filtering, vessel enhancement, detection of thin objects, adaptive morphology. vii R É S U M É Le but de ce travail est de filtrer les objets fins et curvilinéaires dans les images numériques. Leur détection est en soit difficile du fait de leur finesse spatiale. De plus, le bruit, les artefacts de l'acquisition et les occlusions induites par d'autres objets introduisent des déconnexions. De ce fait, la reconnexion des objets fins est également nécessaire. Dans ce but, une méthode hybride à base de dérivés secondes et de filtrage linéaire morphologique est proposée dans le cadre de la théorie espace-échelle. La théorie des filtres morphologiques spatialement variants et des algorithmes sont présentés.

Du point de vue applicatif, notre travail est motivé par le diagnostic, la planification du traitement et le suivi des maladies vasculaires. La première application étudie les malformations artério-veineuses (MAV) dans le cerveau. L'analyse de telles données est rendue difficile par la petite taille, la complexité des vaisseaux couplés à diverses sources de bruit et à leur topologie, sans compter les artefacts d'acquisition et l'hétérogénéité du signal sanguin. Ainsi, nous nous sommes intéressés à l'amélioration et la segmentation des images angiographiques cérébrales dans le but d'aider à l'étude des MAVs cérébrales.

La seconde application concerne le traitement des images X-ray à faible dose utilisées en radiologie interventionnelle dans le cas d'insertion de guides dans les vaisseaux de patients. De telles procédures sont utilisées dans les traitements des anévrysmes, obstructions des tumeurs et d'autres procédures. Dû au faible ratio du signal-bruit, la détection des guides est indispensable pour leurs visualisations et leurs reconstructions. Dans ce travail, nous comparons la performance des algorithmes de filtrage d'objets linéaires. Le but est de sélectionner les méthodes de détection les plus prometteuses dans le cadre de cette application médicale.

A thin, noisy brain blood vessel. There is no low-gradient zone in this object, disconnections are due to noise. . . . . . . . . . . . . . . Figure 5 Volume rendering of MRA volume, from left to right: original image, ground truth, entropy-based measure from [QBD + 09],

Hessian-based measure. (Illustration is from [TdTF + 07]). . . . . . Figure 6 Geomteric vessel models: (a) Hessian eigenvalues' shape-space according to [START_REF] Lin | Stenosis detection using a new shape space for second order 3d-variations[END_REF]. (b) Superellipsoid model by [TdTF + 07].(Illustrations are from according articles.) . . . . . . . . . . . . . . . . . . . . . . . Figure 7 Schema of segmentation with level-sets. (a) 2D case (left -evolving curve, right -curve evolution with level sets). (b) Codimensiontwo curve evolving on tubular isolevel set Γ of C. Right: the tangent to C at p (the normal plane), the external vector → d and its projection onto the normal plane.(Illustration is from [LFG + 99]). . Figure 8 MRA segmentation results using: CURVES (left), capillary active contours (right).(Illustration is from [START_REF] Yan | MRA image segmentation with capillary active contour[END_REF]). . . . . . . . . . . . . . Figure 9 Tracking scheme from [START_REF] Boldak | An improved model-based vessel tracking algorithm with application to computed tomography angiography[END_REF]. (Illustration is from the corresponding publication. . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 10 A sequence of superellipsoidal fits traversing vessel. (Illustration is from [TdTF + 07]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 11 Brain vessels segmentation based on grey-level hit-or-miss transform. Illustration is from [START_REF] Naegel | Grey-level hit-or-miss transforms -Part II: Application to angiographic image processing[END_REF]

Illustration of the concept of firefront. Fire propagates from the border of the object toward its interior, in an isotropic manner.

The point C is a quench point and is the center of a maximal ball (source: [START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 16 Illustration of a maximal ball in the Euclidean case. The ball B max is maximal because no other ball can contain it and simultaneously be included in the set. In contrast, the ball B 2 is not maximal (source: [START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 17 Illustration of medial axes (source: [START_REF] Blum | A Transformation for Extracting New Descriptors of Shape[END_REF]). . . . . . . . . . . . . . Figure 18 Discrete medial axis: (a) the silhouette of a bear; (b) its discrete medial axis, with the 4-connected neighborhood as unit ball; and (c) the ultimate erosion with the same unit ball . . . . . . . . . . . Figure 19 The thinning process on an example 2D shape. Boundary points are marked "B" at the beginning of each iteration and then removed if they are simple (source: [START_REF] Nd Cornea | Curve-skeleton properties, applications, and algorithms[END_REF]). . . . . . . . . . . . . . Figure 20 Example of a distance transform (right) and its original image (left). Figure 21 A graph based on the skeleton. . . . . . . . . . . . . . . . . . . . . . Figure 22 Extraction of junctions and segments from binary image data. The original vascular structure (left) is thinned to its centerline (center), voxels with one or more than 2 neighbors are marked in red as junction points. Segments are identified by following the centerline voxels from one junction to the next one (right). . . . . . Figure 23 An illustration vascular tree using Strahler Ordering (left), where

Slice of a MRA data of the head (left) and a close-up of the AVM. . . Figure 38 An object and its symmetric (left) and asymmetric 1D isthmus skeleton (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 39 An object and its asymmetric 1D isthmus skeleton (left) and its graph based on the skeleton (right) . . . . . . . . . . . . . . . . . . . Figure 40 Rendering of the synthetic vessel-like object with the added Gaussian noise of κ = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 41 Slices of phantom images with four levels of Gaussian noise (from left to right): 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . Figure 42 Comparison of average ROC accuracies for Frangi's or Sato's vesselness and their respective hysterysis thresholding for original phantom image and with four levels of Gaussian noise: 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 43 From left to right: Morpho Hessian with Frangi's normalized and non-normalized vesselness or Sato's measure on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4]. . . . . . . . . . . . . . Figure 44 From left to right: Frangi's and Sato's vesselness on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4]. . . . . . . . . . . . . . Figure 45 Comparison of average ROC accuracies for Frangi's or Sato's vesselness according to results at four discrete scales and one multi-scale (σ: 1, 2, 3, 4, and [1 -4]). . . . . . . . . . . . . . . . . . . Figure 46 Comparison of average ROC accuracies for segmentation results with region growing on SV-enhanced and Frangi's and Sato's vesselness markers, anisotropic diffusion, and original image.

Comparison is performed on the original phantom image and the phantom corrupted with four levels of Gaussian noise: 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 47 Region growing segmentation results on (a)original phantom data with Frangi's and (b) Sato's vesselness (right) seeds (orange) mapped with the ground truth (white). . . . . . . . . . . . . . . . Figure 48 Comparison of average ROC accuracies of region growing segmentation results based on four methods according to scales σ = {1, 2, 3, 4, Frangi's vesselness (yellow) and Sato's with the same thresholds. . Figure 56 Whole head volume segmentation (on the left). Its skeleton with the AVM (yellow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 57 (a) MIP and the selected AVM. (b) AVM with its skeleton. (c) AVM with more vessel connections and (d) their skeleltons. . . . . . . . Figure 58 (a) MIP and the selected AVM. (b) AVM with its skeleton. (c) AVM with more vessel connections and (d) their skeleltons. . . . . . . .

Figure 59

Example of a guidewire image (cropped and contrast-enhanced). . Figure 60 Ridge detectors with second-order (left) and fourth-order (right) derivatives. . .

I N T R O D U C T I O N
The line has in itself neither matter nor substance and may rather be called an imaginary idea than a real object;

and this being its nature it occupies no space.

-Leonardo Da Vinci

The primary subject of this thesis is thin objects analysis applied to medical images. In this chapter we introduce the basic problem of thin objects segmentation and motivate its importance for many applications. We take into account modern medical imaging modalities in a clinical setting. At the same time, as such modalities generate larger and larger image data which gets harder to be examined, it becomes increasingly important to develop efficient procedures to be able to analyse such data in a timely fashion. This motivates the development of more efficient and robust image analysis methods, tailored to the problems encountered in medical images. The motivation of this thesis is directed towards the particular problem of detection and enhancement of thin objects such as vessels or medical instruments inserted into them as catheter and guidewire. However, the generality of the problem can lead to potential usefulness in other areas of image analysis.

thin objects analysis

Greek philosopher and mathematician Euclid was the first one known to geometrically define a curve (or a line) described as "that which has a length without width". This definition illustrates the two most important features of a line: it is extended in some direction and its width is small. It is generally difficult to detect such objects due to these characteristics and even despite progress in resolution improvement, precision and overall quality digital imaging.

Therefore, image processing of linear objects has two main challenges: 1) object enhancement and noise removal for visual quality improvement, 2) segmentation and classification for (semi-)automatic measurements.

However, to perform the first operation, object enhancement, one needs to find these objects in the image. Detecting these objects can be more or less reduced to the image processing operation called segmentation. Segmentation means, in general, partitioning an image between semantically consistent areas into objects of interest on the one hand and the rest, the background, on the other. There are numerous ways to achieve this, however, no general purpose method exists, a fact that requires development of specific application-oriented tools. Therefore, in every application, the primary focus is the definition of objects of interest and the search for their particular features.

medical image applications

The spectacular rise of medical imaging during the 20th century, mainly induced by physics breakthroughs related to nuclear magnetic resonance and X-rays sources and sensors, linked with the general availability of computers, has led to the development of numerous advanced imaging modalities devoted to visualise vascular structures. The analysis of such angiographic images is of great interest for various clinical applications.

Initially designed to generate 2D data, these imaging modalities progressively led to the acquisition of 3D images, enabling clinicians to visualise volumes directly. Even four-dimensional data (3D images changing over time, i.e. movies) can sometimes be acquired. However, such an increase in size and dimensionality provides major technical as well as cognitive challenges. How to store, transmit, visualize and find relevant information are questions related to such data. In addition, the possibly low signal-to-noise ratio and the potential presence of artefacts (such as generated by tomographic reconstruction) make the analysis of such images a rather challenging task. In order to assist radiologists and clinicians, it is therefore necessary to design software tools enabling to extract as well as possible the relevant information from these kind of data.

The highly-desired methods are those that can automatically detect diseases, lesions, tumours and other bodily structures, and visualize them out of this large collection of images. However, the bigger challenge here is to develop methods that can be trusted. This is especially important in medical applications: we expect low false detections rate and even fewer false negatives. Therefore, developing algorithms for medical image analysis requires thorough validation studies to make the results usable in practice. This process involves communication between at least two different worlds: the patient-centered medical world, and the computer-centered technical world, which itself consists of the applied mathematicians whose job is to discover new algorithmic ways to analyse medical data, and the computer scientists who make these algorithms work in practice. The symbioses of these worlds are a rare find and they require significant efforts from all sides to join on a common objective.

The aim of this thesis is to develop image analysis techniques for the assessment of medical imaging applications. In particular, the main application is aimed at segmentation of cerebral vasculature from 3D angiographic image data. The process starts with the detection of the vascular structure and its enhancement. Then, the segmentation and skeletonization procedures generate a computer model of the vessel network which can be used in study of vascular diseases and their visualization. The motivation for this work is to aid clinicians in image studies by providing enhanced, more precise data and its visualization for medical decisions. This involves multi-disciplinary communication.

So, our side goal of the work is to establish links and identify important and relevant medical problems.

The second medical application concerns the treatment of low dose X-ray images used in interventional radiology therapies observing insertion of catheter and guidewires in the vascular system of patients. Guidewire detection assists in such applications as visualization enhancement, guide-wire reconstruction and respiratory motion tracking. To detect guide-wire from such "2D plus time" images, an essential initial step is filtering prior higher level techniques. The quality of the guidewire enhancement algorithm boosts the performance of the final application. In this work, we compare the performance of several line enhancement algorithms for the application of guidewire detection. The purpose of this work is to select a small number of the most promising line detecting methods for further processing and use in medical applications.

Therefore, in this thesis, as a prior, we are interested in detection of thin objects from noisy images for further image enhancement, segmentation and analysis of such objects, e.g. connected in networks. Then, we are interested in applying these techniques in the medical applications corresponding to 2D and 3D problems.

Cerebral arteriovenous malformations and their assessment

Interactions of blood happens at different scales, where large-scale flow is connected to cellular and sub-cellular biology. On the largest scale, blood flows from the heart through arteries to all other parts of the body. Then, arteries branch and become smaller as they reach other tissues until their size decreases to capillaries, which are at cellular scale. The capillary bed is where the exchange process of oxygen and nutrients on the one hand, and waste on the other, actually happens. Then blood is collected through the veins, vessels that bring blood back to the heart. Injured blood vessels in the heart, can lead to chest pain and heart attacks; of the arteries of the neck and brain, to strokes.

An arteriovenous malformation (AVM) is a complex tangle of abnormal arteries and veins that are connected directly without a capillary bed. This "knot" is called nidus (its schematic illustration can be seen in Figure 2).

There is typically high blood flow through the nidus of the AVM, but it is unknown if this flow is a cause or an effect of the abnormal vessels, or both. One hypothesis is that high-pressure blood uses the path of least resistance. Another is that the AVM itself uses blood vessels. In any case, the blood goes through the AVM and not through available capillary beds. This redirection is called a shunt. With time, due to the shunting, the AVM dilates. This dilation weakens veins making them susceptible to haemorrhage and feeder arteries becoming susceptible to aneurysms. A haemorrhage in the brain is a type of stroke where a blood vessel ruptures and bleeds into the brain. Each time blood leaks into tissues, these are damaged. This results in loss of temporary or permanent normal function. The amount of damage depends on how much blood was leaked.

AVM can also occur in other parts of the body: spleen, lung, kidney, spinal cord, liver, intercostal space, iris, and spermatic cord. However, the relationships between a nidus and all aforementioned vessels are more complex. The Arteriovenous Malformation (AVM) can be very simple, featuring a single nidus of one feeder and at least one draining vein. Conversely, they can have multiple compartments: several separated feeders and draining veins. The last notion is referred to as hidden compartments -those that are invisible even on static angiography. Hidden compartments can be detected with serial super-selective digital subtraction angiography or serial high-resolution MR angiography [START_REF] Moftakhar | Cerebral arteriovenous malformations. Part 2: physiology[END_REF]. Therefore, it can be very difficult to effectively picture the 3D intra-cerebral vascular tree.

Moreover, unfortunately, conventional imaging techniques do not supply sufficient information. Digital subtraction angiograms (DSA) provide some useful information, especially when considering time information, but due to the projection overlap, it is difficult to determine the location, connections and the vessel parent-child relationship. Magnetic resonance three-dimensional images are noisy and still do not provide, at present, sufficient connectivity information directly.

Angiography images

The first angiographic experiment dates back from 1896, only a few months after Roentgen discovered X-rays. Vessel visualization examinations entered clinical practice with the invention of Digital Subtraction Angiography (DSA) in the 1970s, which has become and still remains the gold standard for vascular studies.

The main imaging techniques for vessel imaging are: 1) X-ray/digital subtraction angiography (DSA); 2) magnetic resonance angiography (MRA) ; 3) computer tomography angiography (CTA); 4) ultrasound angiography, and 5) angiography as a way of fusing the different medical imaging modalities.

From the clinical point of view, digital subtraction angiography (DSA) is considered the most reliable and accurate method for vascular imaging. However, in DSA, a catheter is used to inject a radio-opaque contrast agent directly into the arteries. The catheterization is an interventional procedure presenting the risk of severe complications. Moreover, the contrast medium, typically containing iodine, can cause serious side effects, such as allergy-like reactions or renal failure in some patients.

As an alternative, MRA and CT offer procedures that are safe for the patient, with reasonably good spatial resolution, even if not currently matching that of DSA. On one hand, the MR and CT techniques lack the ability to locate tiny vessels and the morphological estimation of stenoses and aneurysms. On the other hand, DSA lacks 3D information, which is available with MR and CT techniques.

Here, we make use of MRA data and in the following we present MR principles of work and its angiography techniques. magnetic resonance angiography There exist three major techniques to perform angiography with MRI:

1. Time-of-Flight (TOF). This method uses the fact that, as blood circulates even during the acquisition, its magnetization properties differ from the surrounding, introduction non-moving tissue. For optimal conditions, the image plane has to be orientated perpendicular to the main flow.

Phase Contrast (PC).

This technique is based on the linear relationship between the phase of moving spins and their velocity. The MR signal is directly related to blood flow velocity.

3. Contrast-Enhanced (CE). Like in DSA, contrast agents have been developed to be used in MRA. They are injected intravenously and reduce the relaxation time of the surrounding blood, making vessels appear brighter on images. Often based on Gadolinium, they stay in the vascular bed for only a short period of time (minutes). It is important to note that, unlike DSA, MRA shows the effect of the contrast agent on proton relaxation around contrast molecules, and does not show the contrast agent itself. By timing the injection so that the contrast medium reaches the vessels of interest during image acquisition, it is possible to achieve exclusive arterial enhancement. The total volume of injected contrast medium is called the bolus, and the synchronization of injection and scan time is called bolus chase.

Flow-based techniques, like TOF-MRA and PC-MRA are of limited practical use. They require long acquisition times, and the quality of images depend on the orientation of vessels and blood flow patterns. In addition, the size of the Field-of-View is limited, and many parameters have to be set in order to obtain optimal results. CE-MRA is therefore the preferred technique in most vascular areas.

The quality of the created vascular model depends on several factors, including field strength, pulse sequence supported by specialized hardware, segmentation and/or post processing, level of expertise and attention to details of the operator. The field strength is important and generally, the higher the field strength the higher spatial resolution and/or signal to noise ratio. More vessels are discernible on 3T scans than on 1.5T ones. Similarly, on 7T more vessels are visible than on 3T scans; moreover, some single vascular spots on 3T are in fact two neighbouring vessels on 7T. However, the field strength itself is not a sufficient factor to generate good models.

Compared to the size of the nidal vessels (typical diameter 0.1 mm) the in-plane spatial resolution is of the order of 0.5 mm and thus not all vessels can be characterized.

In the case of MRI, the noise is typically modelled by a Ricean distribution, which, for larger means, or equivalently high Signal-to-Noise Ratios (SNR), can be approximated by a Gaussian.

The increased sensitivity, higher resolution, and additional contrast mechanism (e.g., by employing multiple coils) available at higher field strengths offer advantages when acquiring MRA volumetric data. 3D TOF MRA is the dominant MR technique for imaging the intracranial arterial vasculature.

The drawback of the MRA method for arterial shape reconstruction is its inability to capture a volume where blood is stagnant, since the method is based on detecting a moving fluid. MRI is a very elaborate technique. Many sources of defects in an image exist, for instance due to the scanner hardware, the imaging sequence, and even the subject inside the scanner. These defects are called artefacts. Some of them are relevant to our study.

The noise present in the MRI image is primarily due to thermal noises in the patient [START_REF] Nowak | Wavelet-based Rician noise removal for magnetic resonance imaging[END_REF]. Noise present in the raw, complex MR signal acquired in the Fourier domain presents a Gaussian distribution, but the transformation to a magnitude image changes the Gaussian distribution of the data to a Rician distribution [SDDVdL + 99]. Noise mainly depends on: 1) voxel size 2) acquisition time 3) main magnetic field intensity.

In MR data, a voxel actually contains many protons and the signal is averaged over the whole voxel volume. When a voxel contains more than one tissue type, the signal contribution of this voxel is a mix of the signals from the different tissue types the voxel contains. This is called Partial Volume Effect (PVE). It is a common property of every imaging system. However, as it causes blurred boundaries and signal variations within a volume, it is an important effect that has to be taken into account by the image analysis process.

Coronary heart disease

Vascular heart diseases are the first causes of death in Europe. One of the main reasons is stenosis, calcifications causing vessel erosion or dilation. The consequences are usually serious especially in case of damaged coronaries, arteries that provide oxygen to heart.

X-ray imaging makes it possible to diagnose and assist a treatment in such failures. During the procedure, the clinician inserts one or more flexible narrow tubes called catheter through small skin incisions in groin. These tubes are pushed through blood introduction vessels into the vessels of the heart, neck or brain. While viewing an x-ray monitor, called a fluoroscope, the doctor steers the catheter through the blood vessels. A fluoroscope machine, also called a C-arm, is an arc shaped piece of equipment that generates x-rays from one side and photographs them on the other side (Fig. 3).

Such procedures are used in aneurysm treatment, tumour embolization and other clinical procedures.

Interventional images

For diagnostic purposes, the X-Ray imaging system uses relatively high radiation resulting in well contrasted images. However, interventional exams are carried out under limited radiation (5 to 10 times less than in diagnostic images), especially due to the fact that the acquisition rate is typically between 7.5 and 30 images per second.

The critical importance is on accurate positioning of the guide wire and the catheter with respect to the vasculature, taking into account patient motion and the low signal to noise ratio of the images. The limited quality of the image is defined by the low dose constraint in fluoroscopy in order to minimize the radiation exposure of the patient and radiologist.

Because of the physics of X-Rays, the clinical images are transparent. There can be seen the tissues of heart superimposed to the diaphragm, the lungs, the spine, the ribs, and the interventional devices.

Apart from these transparency artefacts, there are some sources of noise in such images. The first one is due to the physics of X-rays and the sensor characteristics, leasing to a Poisson or Poisson+Gaussian noise model and which is spatially correlated by the detector of these rays (flat panel) and is temporarily uncorrelated (from one image to another).

Other sources of noise are electronic, deviation constant noise (white). Another type of noise is called structure noise which corresponds to constant structures present in the image but has no relation to the patient. Those are, e.g. defects of the detector that are not well corrected.

outline

The rest of this thesis is organised as follows. In Chapter 2 we will present the review of thin objects filtering and segmentation methods and their combinations. The methods are organised by their mathematical models and extraction schemes. We present the methods principal properties and characteristics. The review of the literature on vascular network representation and analysis is presented in Chapter 3, where the topological notions and utility are discussed, continuing into the user-friendly aspect of their utility -visualization.

In the following Chapter 4, we will present our hybrid morpho-Hessian filter for thinobject detection and filtering method based on the second-order directional Gaussian 1.3 outline derivatives in the Hessian form. Further, we are making use of several probability measures for the membership of pixels to thin objects.

For the enhancement and reconnection of tubular objects we are using the spatially variant morphology, for which the theory and corresponding algorithms are thoroughly formulated with the purpose of filtering an image differently at various positions. We then propose an efficient algorithm for spatially-variant mathematical morphology operations.

In Chapter 5, we combine these two filtering methods and perform their performance evaluation with different parameters first on a 3D simulated data and then applied to clinical angiographic data. The Hessian-based method is evaluated according to two popular vesselness measures. According to the best filtering results, we perform further vascular objects segmentation and skeleton extraction. The topological representation of the vascular network is aimed to help in visualization and study of the cerebral vasculature and AVM.

In the following Chapter 6, we compare the performance of several line enhancement algorithms for the application of guidewire detection for very low signal-to-noise ratio X-ray images. The purpose of this work is to select a small number of the most promising 2D line detection methods for medical applications.

In this thesis, we propose a new hybrid filtering method based both on local intensities and on local neighbourhood tubular model utilizing the local directions for thin objects enhancement and reconnection with the spatially-variant morphological operations. This filter being a novel combination of methods performs like a reconnecting diffusion filter. We introduce a segmentation and vascular network analysis strategy with an aim of aiding the clinicians in their angiography data studies. We also compare the specialised thin objects detection methods from noisy images applied for guide-wire segmentation.

The human mind delights in finding pattern -so much so that we often mistake coincidence or forced analogy for profound meaning.

-Stephen Jay Gould Line-like objects are widely searched for in image analysis problems. In biomedical images those are vessels, neurones, bones, muscles and skin fibres, cracks of tongue, motion and structures of cells. Along with a multitude of other applications, such as satellite images, print lines, road and material cracks, cosmic space and aerial images features extraction. Our work is motivated by diagnosis, treatment planning and follow-up of vascular diseases.

The interest can lie in a whole object representation (e. g.fibre surface reconstruction) or in a certain characteristic of it (e.g. topology, size, curvature). Moreover, extraction of linear objects and segments is an essential step for the segmentation of whole line networks and/or parts of larger objects.

It is generally difficult to filter thin objects, even despite the progress in resolution improvement, precision and overall quality of digital imaging. In order to improve linear object, a detection of such features is appropriate. Even more so, due to the inherent properties of these objects: i.e. thin, long and curvilinear, they are prone to noise. Real-world and image reconstruction artefacts often imply disconnection along the object.

Indeed, image processing of linear objects has two main dual challenges: (i) object enhancement and noise removal for visual quality improvement, (ii) segmentation and classification for (semi-)automated measurements. Both tasks are difficult due to the specificity of the data. Questions on how to enhance thin vessels while simultaneously not increase the level of noise too much, how to remove noise and not to lose useful information are still relevant.

Most segmentation methods are designed for large and homogeneous regions. Here, we define thin objects in image as semantically consistent objects that exhibit at least one dimension much smaller than the others. We focus particularly on elongated thin objects, locally curve-or line-like.

In this chapter, in the first section, a short review of filtering methods of such objects is presented. In the following sections, more curvilinear objects filtering techniques are reviewed as a part of or in connection with segmentation methods.

thin objects filtering

The filtering techniques of curvilinear objects can be divided into two categories : (i) those that eliminate noise while preserving objects, (ii) and others that enhance curvilinear objects while avoiding noise amplification. In both cases, the prior detection of such objects is implied. All of the curvilinear detection methods can be classified in two types: bottom-up (local) and top-down (global). Commonly "bottom-up" methods include preprocessing and segmentation in order to extract objects edges ( [START_REF] Canny | Finding edges and lines in images[END_REF]) and centrelines (ridges as in Haralick [START_REF] Haralick | Ridges and valleys on digital images[END_REF]). The detection methods based on deformable models or similar can be classified as "top-down" ones.

The "bottom-up" class of methods can in its turn be classified in two categories: (i) those that utilize some filters to extract boundary or ridge of the vessels, followed by further refinement (e. g.[Kol95, FNVV98, SNA + 97]); (ii) and those that employ tracking strategy by some given or detected seeds in the vessels.

An other way to classify line-extraction approaches is based on how the image is considered. The first approach detects lines on a pixel-by-pixel basis by only taking in account the intensity of the image, e. g.: 5. Ridge detection [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF] 6. Some wavelets [START_REF] Unser | Higher-order Riesz transforms and steerable wavelet frames[END_REF] These methods are based on oriented filter banks and consist in finding the orientation corresponding to the maximum response of a filter bank. Each filter is obtained from a rotation of a basis filter. The accuracy and selectivity of these methods usually depend on the number of filters, on the size of the computing support and on the kernel basis filters [START_REF] Perona | Deformable kernels for early vision[END_REF][START_REF] Chen | Orientation space filtering for multiple orientation line segmentation[END_REF].

The second approach for local line-detection considers the image as a function f(x, y) and lines are detected as ridges and are opposites of valleys. Examples of those include:

1. Derivative-based methods: Indeed, classical filters in the literature assume locally isotropic objects, i.e. such that a "large enough" window will still fit into them (this is the case for the median, averaging, linear convolutions, morphological filters with standard structuring elements) [START_REF] Gonzalez | Digital Image Processing[END_REF][START_REF] Heijmans | Morphological filters for dummies[END_REF]. Those that do not make this assumption still filter only in areas of low gradient in a region of interest [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Smith | Susan a new approach to low level image processing[END_REF][START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. Elongated objects may in fact present no part with a suitable low gradient, due to both noise and object edges, as is readily apparent in Figure 4.

For elongated objects, it is often assumed that one dimension is long, everywhere tangent to the object. Filtering can be applied along this direction without necessarily destroying the object entirely. Directional (matched) filters, Hough, and Radon transforms can be used to detect straight lines and other features of known shape, but curvilinear objects are much more difficult to trace.

Orkisz et al. [OBM + 97] used a kind of filter bank called "sticks", which can be seen as a set of directional structuring elements. Similar approaches were also proposed by [CJOJ99, KS96, CH95, DPD95a, DPD + 95b]. The last two combined the outputs of directional operators without an explicit extraction of the vessel local orientation. The main disadvantage of the methods in this category is that they can hardly detect vessels thin objects filtering and segmentation methods in a wide range due to the fixed scale analysis. Although they can be extended to multiple scales by using sticks of variable length, the computation time would increase significantly.

The proposed approach by [TKL + 09] incorporates the use of line-like directional features present in an image, extracted by a directional filter bank, to obtain more precise Hessian analysis in noisy environment and thus can correctly reveal small and thin vessels.

Such filters have their analogues in the mathematical morphology field, they are described in section 2.2.9.

Isotropic diffusion filtering removes noise as well as thin objects and hence is not suited for vascular images. Several authors proposed anisotropic approaches where local orientation of structures is first estimated and filtering is locally done along this orientation in order to preserve small objects and not to blur boundaries of larger structures. This approach implies a two-stage processing: (i) orientation estimation and (ii) filtering.

Non-linear anisotropic filtering has also been applied to vessel enhancement [KMA97, Kri02, OBM + 97]. In [START_REF] Krissian | Directional anisotropic diffusion applied to segmentation of vessels in 3d images[END_REF], smoothing was carried out by anisotropic diffusion in the direction of the least principal curvature. One of the methods proposed in [VDS + 92] applied a morphological operation of crest detection, using directional (linear) tools, and selected the direction giving the strongest response.

In the case of the linear scale-space framework, which should be suitable to this problem, edge and ridge detection methods were proposed [START_REF] Bigün | IEEE First International Conference on Computer Vision[END_REF][START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF][START_REF] Danielsson | A new shape space for second order 3d-variations[END_REF] utilising the Hessian or the structure tensor [START_REF] Deguchi | Detection and enhancement of line structures in an image by anisotropic diffusion[END_REF]. Anisotropic diffusion is often used for this task [Per98, MVN06] using the tensorial information for diffusing only within the object, while becoming isotropic outside the oriented objects.

Directional second derivatives can also be brought together in a Hessian matrix, so as to exploit the matrix's eigenvectors and eigenvalues . These methods are discussed in detail in a further section. Similar geometrical considerations concerning the principal curvatures led to the use of the Weingarten matrix eigenvalues. However, a common problem of the methods based on derivatives is their sensitivity to noise.

There were several attempts to characterize the anisotropic properties of the vascular images by non-linear combinations of outputs of directional filters. In [START_REF] Chen | An algorithm for MR angiography image enhancement[END_REF], a set of directional mean-filters was used. The difference between the strongest and the weakest response was used as anisotropy measure. Indeed, the mean intensity along a vessel should be larger than in perpendicular directions. Poli [START_REF] Poli | An algorithm for real-time vessel enhancement and detection[END_REF] proposed a computationally efficient algorithm based on a set of linear filters, obtained as linear combinations of properly shifted Gaussian kernels, sensitive to vessels of different orientation and radius.

detection and segmentation methods

One of the main image analysis tasks is segmentation or classification, which are both taken to mean distinguishing regions of the image belonging to different objects. In practice, this must be performed with only a few facts known about the objects present in the image. To these aims, the importance of edge detection has been demonstrated multiple times. Some of the theories insist on finding edges at an early stage on a brightness representation of the image, while segmentation and other stages follow afterwards. Others find discontinuities and a potential focus of attention as a bi-product of the perceptual organization process which is based on a ridge detection [START_REF] Subirana-Vilanova | Multi-scale vector-ridge-detection for perceptual organization without edges[END_REF].

Indeed, there exist comparatively few segmentation strategies proposed that manage without edge information. Those fall into two classes. In the first class, one of the classical methods is based on region growing [HR78, HP74, HS85, Cle91, AB94]. These methods are based on placement of several "seeds" defining different classes or objects in the images and then "growing" them until a complete region is found. The growing decisions are based on neighborhood intensities. This scheme has two drawbacks: it is local and shape constraints can be difficult to incorporate.

The second class of segmentation strategies that work without edges are based on computations that find discontinuities while preserving certain region properties such as smoothness [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF][START_REF] Terzopoulos | Regularization of inverse visual problems involving discontinuities[END_REF][START_REF] Poggio | Parallel integration of vision modules[END_REF]. These methods are scale dependent and in fact in some cases depend on reliable edge detection. The first scale studies on the discontinuity level [WT83, [START_REF] Koenderink | The structure of images[END_REF][START_REF] Perona | Detecting and localizing edges composed of steps, peaks and roofs[END_REF] do not explicitly represent regions.

Moreover, it is a known fact that edges alone are not enough in order to distinguish objects that are difficult to find.

Considering that two main line distinguishing features are width and orientation, it is suitable to analyse these two characteristics to start off the segmentation of lines. The early methods applied non-linear pixel operators by looking at neighbouring pixel intensities in order to extract roads from aerial images [START_REF] Fischler | Detection of roads and linear structures in low-resolution aerial imagery using a multisource knowledge integration technique[END_REF][START_REF] Rosenfeld | Edge and curve detection for visual scene analysis[END_REF]. However, these methods did not take into account the fact that objects can be represented at different scales.

The segmentation of the vascular structures is a particularly challenging task. This is due to its particular image acquisition modalities, three dimensional nature, and presence of other tissues and quantity of image noise and artefacts in proportion to vessel signal. Moreover, physical and anatomical properties of the tissue are highly variable in size, appearance, geometry and topology, and even more so in pathological cases like aneurysms, stenoses, calcifications and arteriovenous malformations.

More detailed reviews on vascular segmentation can be found in [LABFL09, QBD + 09, KQ04, SLRL02a, SLRL02b].

Becoming more and more popular are combinations of several techniques, e.g. as a pre-, post-filtering steps, in order to obtain the desired segmentation result. From an image processing point of view, segmentation consists of partitioning an image into an object, i.e., a structure of interest, and a background, i.e., the remainder of the image volume.

The segmentation of vascular structures from 3D images is an even more difficult task. Here, the notion of segmentation is considered in a wide sense. In the context of angiographic imaging, we consider that vessel segmentation embeds (i) methods that detect either whole vessels (i.e., their lumen and/or walls) or their medial axes, and/or (ii) methods that perform low-level processing or high-level knowledge extraction (e.g., vein/artery discrimination [LUSO01, TS02, BSVN03] or vessel labelling [CFM01, HEM + 99]).

We also consider some methods which could be classified as filtering ones, since their purpose is to perform vessel enhancement, which consists mainly of denoising, but also of vessel reconnection (for example in case of stenosis, or of signal loss [DLPB99, PRB + 05a]).

As discussed above, the difficulty to perform vessel segmentation is due to the sparseness of data, and the possible presence of irrelevant signal (other tissues, artefacts or noise). The choice of a segmentation method is often linked to the type of considered images, the vessel(s) to study and the clinical purpose. The next section discusses different methodological segmentation strategies relevant to the problem at hand.

Surveys of vessel segmentation methods

Several surveys devoted to 3D vascular segmentation have been proposed during the ten last years or so. The oldest ones [START_REF] Suri | A review on MR vascular image processing: Skeleton versus nonskeleton approaches: Part II[END_REF][START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] are now globally incomplete, the reader may refer to the most recent ones [START_REF] Orkisz | Models, algorithms and applications in vascular image segmentation[END_REF][START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF] for a more exhaustive overview of the literature. The survey proposed in [START_REF] Suri | A review on MR vascular image processing: Skeleton versus nonskeleton approaches: Part II[END_REF] focuses on vessel segmentation from MRA images, and divides them into skeleton methods (with an interest in medial axes) and non-skeleton ones (that aim to segment whole vascular structures). These two families are then refined, based on the involved image processing strategies. Part I of this survey [SLRL02a] also describes MRA acquisition techniques.

An other classification is proposed in [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] which deals more generally with vessel segmentation from any kind of data independently of their dimension, or acquisition technique. The review of image analysis for angiographic image of [START_REF] Orkisz | Models, algorithms and applications in vascular image segmentation[END_REF] gives a comprehensive, but non-exhaustive introduction on methods per clinical vascular application. While one of the most recent survey [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF] mainly refers to 3D vessel segmentation from MRA and CTA, it divides its description into (i) the a priori information which can be used for segmentation, (ii) the basic tools using this information for detecting vessels, and (iii) the methodological frameworks involving these tools, as well as a discussion on pre-and post-processing considerations.

In the next section, we introduce the segmentation methods divided into eight main families corresponding to the main image processing strategies on which they rely: region-growing, deformable models, model-based filtering, path finding, vessel tracking, statistical approaches, differential analysis, and mathematical morphology.

Region growing methods

Region-growing has been one of the first strategies considered for image segmentation [START_REF] Zucker | Region growing: Childhood and adolescence[END_REF], and in particular medical/angiographic ones, dedicated to seeded region-growing segmentation. Basically, region growing relies on two elements: one (or several) seed(s) [START_REF] Adams | Seeded region growing[END_REF] assumed to belong to the structure of interest to be segmented, and a propagation criterion, enabling to segment the object from the seed, by iterative addition of adjacent voxels.

In the case of vessel segmentation, seeds are generally defined interactively inside vessels. The possible definition of several seeds can straightforwardly lead to an application of region-growing to vessel separation, and in particular, to vein/artery discrimination. In such a case, a set of seeds is defined for arteries and veins, respectively. A competitive region-growing is then performed, based on ad hoc propagation criteria (e.g., a measure of grey-scale connectedness in [START_REF] Tizon | Segmentation with gray-scale connectedness can separate arteries and veins in MRA[END_REF]). Note that, by duality, regiongrowing also provides solutions to segment vessels by skeletonisation. In such a case, the growing process starts from a seed being a subset of the background (which can then be automatically defined), and generally includes topological constraints in the propagation criterion [DLPB99, PRB + 05b].

The seeds can be detected automatically, especially in the case where they constitute the root of a vascular tree [START_REF] Naegel | Using grey-scale hit-or-miss transform for segmenting the portal network of the liver[END_REF]. The propagation criterion is commonly based on intensity criteria, related to the high-intensity vascular signal.

However, more sophisticated properties can also be embedded in this segmentation strategy. In particular, it has been proposed to consider a priori knowledge related to the shape and size of the vessels to be segmented [START_REF] Naegel | Grey-level hit-or-miss transforms -Part II: Application to angiographic image processing[END_REF], or to their topology [PRB + 05a]. The correctness of the orientation of the vessels during the segmentation process has also been considered by proposing "wave propagation" strategies [ZJE + 95], which aim to constrain the segmentation front to remain normal to the vessel axis. It may be noticed that this kind of approaches has been further used for vessel tracking methods (discussed hereafter in the section). The concept of wave propagation has also further led to the development of methods related to both deformable models (level-sets) and path-finding approaches, namely, fast-marching methods [START_REF] Manniesing | Local speed functions in level set based vessel segmentation[END_REF].

Region-growing methods rely on a simple algorithmic framework, which make their development and use quite easy and induces a low (generally linear) computational cost. In addition, they guarantee termination which is not systematically available for other non-monotonic strategies. However, the connectivity hypothesis intrinsically associated to this strategy constitutes a weakness, since the method may fail in segmenting vessels in case of vascular signal loss (due to partial volume effect, or flowing artefacts, for instance). A contrario, the use of a criterion being too permissive may lead to leakage phenomena, and a final over-segmentation of vessels, requiring to anticipate this effect [MSvdG + 07].

In this context, region-growing methods have often been preferentially devoted to the segmentation of large and/or well-contrasted vessels (for which intensity and connectivity hypotheses are generally reliable). Otherwise, dual object and background competitive region-growing has been largely used since proposed in [START_REF] Yi | A locally adaptive region growing algorithm for vascular segmentation[END_REF].

A spatially adaptive propagation can bring advantages, such as topological coherence, e.g., in order to detect bifurcations [LRSB03, MSvdG + 07]. Branches can be segmented and reconnected easier by a local approximation of the vessel centerline, which can be derived from the front centers [START_REF] Lorenz | Simultaneous segmentation and tree reconstruction of the coronary arteries in MSCT images[END_REF], by skeletonization of the local binary mask [START_REF] Manniesing | Local speed functions in level set based vessel segmentation[END_REF], or by path tracking [START_REF] Kirbas | Vessel extraction in medical images by 3d wave propagation and traceback[END_REF]. The last method is closely related to minimal path methods [START_REF] Appleton | Efficient path openings and closings[END_REF].

Region-growing methods largely depend on the robustness of their features. To handle potential under-and over-segmentation issues, adaptive parameters such as inclusion thresholds are desirable [START_REF] Lorenz | Simultaneous segmentation and tree reconstruction of the coronary arteries in MSCT images[END_REF][START_REF] Manniesing | Local speed functions in level set based vessel segmentation[END_REF][START_REF] Yi | A locally adaptive region growing algorithm for vascular segmentation[END_REF]. Such adaptive schemes rely on iterative, adaptive nature of the propagation process. In [LRSB03, ESS + 04] parameters are adapted for each branch independently, where at the end of each vessel the algorithm searches for possible vessels around to reconnect with intensities adapted to the neighborhood. Even though adaptive criteria are normally studied heuristically, they are simple and good examples of collaborations between the extraction schemes, the features and the appearance models.

Derivatives-based methods

Vessels are generally bright structures among a dark background. If an image is viewed as the discrete analogue of a function from R 3 to R, vessels then appear as the maxima of this function. Consequently, it may be possible to detect them by analysing the differential properties of the image.

In order to deal with the discrete/continuous issue involved by this strategy, the considered (discrete) image is convoluted with a series of Gaussian derivatives with different standard deviations and in different directions, the obtained responses being combined into a matrix.

In the case of first derivatives analysis, this matrix, which is the covariance matrix of gradient vectors [START_REF] Bigün | IEEE First International Conference on Computer Vision[END_REF][START_REF] Agam | Vessel tree reconstruction in thoracic CT scans with application to nodule detection[END_REF], is called the structure tensor.

First-order derivatives

Except for vessel segmentation, the first derivatives have also been involved in diffusion filtering, which consists in propagation of information in the orientations suggested by these derivatives [START_REF] Mahadevan | Robust model-based vasculature detection in noisy biomedical images[END_REF].

The structure tensor matrix can be analysed with the purpose of analysis of the local distribution of gradient vectors [START_REF] Weickert | Coherence-enhancing diffusion filtering[END_REF]. Some works show the ability to extract more than one principal direction [START_REF] Agam | Vessel tree reconstruction in thoracic CT scans with application to nodule detection[END_REF][START_REF] Agam | Expectation maximization approach to vessel enhancement in thoracic ct scans[END_REF]. Therefore they can distinguish other shapes than tubular, i.e., masses and surfaces.

With the first derivatives it is possible to obtain directional information which can be used for a propagation of information and/or smoothing, diffusion filtering. An isotropic smoothing has been performed with first derivatives as in [START_REF] Mahadevan | Robust model-based vasculature detection in noisy biomedical images[END_REF][START_REF] Sofka | Retinal vessel centerline extraction using multiscale matched filters, confidence and edge measures[END_REF] within homogeneous regions in order to filter-out noise. This procedure can also blur-out borders and small objects of objects if those are not detected well.

Anisotropic filtering can help preventing this by propagating the information along the vessels, which again requires correct detection of the principal object direction [START_REF] Deguchi | Detection and enhancement of line structures in an image by anisotropic diffusion[END_REF]. Some works [BB08, VMA08] used the gradient vector field [START_REF] Verdú-Monedero | Spatially-variant directional mathematical morphology operators based on a diffused average squared gradient field[END_REF] in a diffusion scheme [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF]. This scheme regularizes and anisotropically diffuses gradients of the object boundaries, which is absent at the medial axes of the object. This has been solved by integrating second-order derivatives discussed further. Vectorial gradient flux was also used for anisotropic diffusion for thin object filtering and enhancement [START_REF] Krissian | Flux-based anisotropic diffusion applied to enhancement of 3d angiograms[END_REF]. Gradient oriented flux has been recently used in more complex segmentation schemes as in deformable methods described in Section 2.2.5.

As discussed further, the first-order derivatives can be used in cross-sectional vessel analysis [KMA + 00, WNV00, FPAB03, AGT08] as oriented features.

Second-order derivatives

An other popular approach in the detection of vascular patterns is the use of secondorder derivative information to characterize the local image geometry. Canny [START_REF] Canny | Finding edges and lines in images[END_REF] proposed to take the second-order derivatives of the Gaussian function as a filter and detect lines at points where the convolution of the function with the lines gives a maximum response.

The Hessian matrix is the most common tool to capture such information. Its extension to multiple scales was proposed in [START_REF] Koller | Multiscale detection of curvilinear structures in 2d and 3d image data[END_REF]. In [START_REF] Antoine Maintz | Evaluation of ridge seeking operators for multimodality medical image matching[END_REF] there has been performed an evaluation of Gaussian derivative functions and models. Steger [START_REF] Steger | Extracting curvilinear structures: A differential geometric approach[END_REF][START_REF] Steger | An unbiased detector of curvilinear structures[END_REF] and Hladuvka [START_REF] Hlad | Exploiting the hessian matrix for content-based retrieval of volumedata features[END_REF] have proposed first functions for curvilinear structures detection based on Hessian matrix eigenvalues.

An image can be convoluted with second-order Gaussian derivatives with their reponses forming Hessian matrix. The main idea behind eigenvalue analysis of the Hessian matrix is to extract one or more principal directions of the local structure of the image. This gives the direction of the minimal curvature, the principal one in the tubular structure and a high curvature in the vessel cross-section plane, which makes the filter more efficient to the contrary with multiple orientations line filters.

Compared with the image gradient, whose response is independent of the shape and local structures of boundaries, the Hessian matrix can capture the shape characteristics of objects, such as tubes, planes, blob surfaces or noise. In particular, different combinations of these eigenvalues were proposed to enhance points likely to belong to vessels. Indeed, tubular structures should give rise to the associated eigenvector that is tangential to the vessel axis and to the associated eigenvectors that lie within the plane locally orthogonal to the vessel. Different combinations of these eigenvalues were proposed to enhance points likely to belong to vessels. If appropriately designed and applied at multiple scales, such combinations, often called vesselness function, should give the strongest response at one particular scale corresponding to the plate-, blob-like and/or tubular objects [? FNVV98, KMA + 00] (see Figure 6(a) for an example of an eigen-value elliptic space).

These methods can be performed in multi-scale frameworks in order to detect objects of different sizes. It has to be noticed that the choice and number of the considered scales is particularly important in such methods and are further discussed in the section.

Based on second-order derivative filters and their directional information, there is a fair amount of works performing an anisotropic diffusion for thin objects enhancement [Wei99, MVN06, NRTW07]. The amount and orientation of diffusion depend on the local vessel likeliness. It is also shown that using vessel-enhancing diffusion as a preprocessing step improves level-set based segmentation of the cerebral vasculature, in particular segmentation of the smaller vessels of the vasculature. Recently, Descoteaux [START_REF] Descoteaux | Bone enhancement filtering: Application to sinus bone segmentation and simulation of pituitary surgery[END_REF] reported a multi-scale bone enhancement measure by introducing a sheetness measure from multi-scale eigenanalysis in the geometrical flow framework.

However, second-order based methods have some weaknesses, e.g. scale-space, bifurcations detection and efficiency, that are discussed further. They have difficulties especially when assuming one dominant direction, which does not necessarily exist in branching points or other extreme cases like aneurysms and stenoses.

Even then, derivative-based methods provide efficient solutions for detecting vessels, especially in a multi-scale framework, and have then often been considered for the design of segmentation methods based on model filtering (presented in further sections) or for the guidance of deformable models, for instance.

scale-space

The second-order derivatives techniques commonly are implemented in multi-scale frameworks to detect vessels of different sizes. Hessian matrices are classically computed in a Gaussian linear scale-space, through convolutions with Gaussian derivatives of different standard deviations [START_REF] Lindeberg | Scale-space theory in computer vision[END_REF]. Scale selection is most often performed by selecting the maximum response over multiple scales.

Due to the fact that size of vessels ranges a lot in different parts of the image. The choice and number of scales is especially important to such groups of methods that perform local filter matching: model-based and derivative-based methods. Too small scales emphasize small details, while producing spurious responses to noise, tangent objects and other features. Too large scales deform original objects at detection and lose details.

On the other hand, it is computationally expensive to process along all possible scales. For the derivative-based techniques, it has been proposed to use dynamic scale selection [Lin93, Kol95, KGSD95, LCB + 97, SRV01, Xu04, SED + 04]. They permit not only to detect the local width, but also to track objects whose width changes along the structure.

The dynamic scale methods reduce filtering time by using smaller scale measurements when appropriate. In [SED + 04] an enhanced efficiency was achieved via dynamically adjusting the filter scale, although segment terminal points were identified by a human operator.

Moreover, scale is correlated to quantity of noise: at smaller ones, there is more of it. This knowledge has been incorporated in [START_REF] Wyatt | Minimum reliable scale selection in 3d[END_REF].

In time-of-flight MRA data intensity ridges inherently exist along the central axis of small vessels, while large vessels may appear brighter near their edges because of laminar flow. In contrast, in CT data, the intensity profile of a large vessel may be flat. It has been said in [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF] that use of a larger scale straightens the curves in a centreline and increases the influence of neighbouring objects.

One more disadvantage due to the scale space appears for tangent objects. The boundaries between tangent objects fade on higher scales. But they have a tendency to get blurred for the smallest objects even on the lowest scales.

Another serious problems in the application of the linear Gaussian second derivative filters of ridge detection is their response to other features, such as edges or sheets in 3D. Already in early works, there has been suggested to use a nonlinear operator that combines the responses of two edge-detectors on both sides of a ridge [START_REF] Koller | From data to information: Segmentation, description and analysis of the cerebral vascularity[END_REF].

[LCB + 97] used an edge-indicator to suppress the response to edges. [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF] proposed to compute positions of ridges in terms of extrema of one operator and scales in terms of extrema of another operator.

In [START_REF] Majer | On the influence of scale selection on feature detection for the case of linelike structures[END_REF] they have used the parameter γ by Lindeberg [START_REF] Lindeberg | On Scale Selection for Differential Operators[END_REF][START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF] to give a priority to ridges and not edges. Figure 27 illustrates that in the variable scale setting the response to edges of a second order Gaussian derivative operator can be turned off simply by a suitable choice of γ parameter. To better illustrate this point the shortest ridges of the grass are not displayed in Fig. 27c. 27b displays all the detected ridges. Moreover, although these ridges cross each other in the projection of Fig. 27, they in fact occur at different scales.

Other way to treat edges is by integrating first-order derivatives responses, as in [START_REF] Cai | Multi-resolution vessel segmentation using normalized cuts in retinal images[END_REF][START_REF] Bouattour | Vessel enhancement in 2d angiographic images[END_REF], where it is suggested that first derivatives detect better thinner vessels than second ones. So they apply those, then after the response analysis, for the candidate vessels they perform a hysteresis thresholding and apply a tracking strategy in order to get all the other small vessels.

bifurcations Another issue for derivatives-based methods are vessel bifurcations, which, at some scales, do not show any principal direction, and at others, they can induce high curvature values along more than one principle direction. This is especially important in methods like tracking, directional filtering and centreline extraction. The Hessian-based model filters tend to attenuate junctions since junctions are characterized as blob-like structures. Junction suppression leads to the discontinuity of the vessel network, which is of course undesirable.

The first geometric bifurcation model was first introduced by [START_REF] Murray | The physiological principle of minimum work: I. The vascular system and the cost of blood volume[END_REF] with relationships between branching angles and vessel widths based on physiological optimality conditions. Other precise models were proposed by [START_REF] Antiga | Robust and objective decomposition and mapping of bifurcating vessels[END_REF][START_REF] Flaaris | Method for modelling cerebral blood vessels and their bifurcations using circular, homogeneous, generalised cylinders[END_REF], which require a pre-segmented centreline and surface mesh for their optimization.

More recent simpler models were introduced in works by Agam [START_REF] Agam | Vessel tree reconstruction in thoracic CT scans with application to nodule detection[END_REF] with a filter model which is based on the correlation matrix of the regularized gradient vectors (first-order derivatives). In the inertia moments space, Hernandez [START_REF] Hernandez | Three-dimensional segmentation of brain aneurysms in CTA using non-parametric region-based information and implicit deformable models: Method and evaluation[END_REF] has defined bifurcations as a contrario case of a vessel model .

On the other hand, although the image gradient is more general in handling structures of different shapes, due to the presence of intensity inhomogeneity such as bias field or overlapping objects, the intensity difference between vessels and background regions are not consistent. The boundaries of the low contrast vessels do not give large gradient responses for those vessel boundaries.

efficiency The multi-scale methods represent substantial computational costs, which is due to the iterative convolution of each pixel of the image. In order to compute eigenvalues of the Hessian matrix, the matrix itself represents four-times original image in 2D and nine-times in 3D. However, considering that it is symmetric matrix, there exists a relatively simple analytical solution to make the code much faster.

In [START_REF] Orłowski | Efficient computation of Hessian-based enhancement filters for tubular structures in 3D images[END_REF], it has been proposed to avoid the computation of the eigenvalues for as much as half of the voxels by checking the sign of the trace of the Hessian. But such closed-form solutions are prone to numerical instabilities, so Orlowski [START_REF] Orłowski | Efficient computation of Hessian-based enhancement filters for tubular structures in 3D images[END_REF] suggest that the eigenvalues for the remaining voxels be computed either by use of an iterative algorithm optimized for 3 × 3 symmetric real matrices, or by use of a hybrid algorithm that switches from non-iterative to iterative solution when risks of numerical instability occur, as suggested in [START_REF] Kopp | Efficient numerical diagonalization of hermitian 3 x 3 matrices[END_REF].

This problem is then solved by a trade-off between the quality of the filtering and gain in efficiency.

Model-based filtering

In general, vessel models can be used as prior for the final segmentation. Here they are described in the order of the increasing complexity. We start with intensity profiles and simple geometric models useful.

Intensity models

Intensity models, which are among the simplest ones, strongly depend on the considered imaging modality. They can integrate brightness, contrast and gradient prior, but also imaging properties, like intensity ranges or intensity variation based on location, or even noise distribution [START_REF] Agam | Probabilistic modeling based vessel enhancement in thoracic CT scans[END_REF] (see also the paragraph Statistical analysis, for a discussion on noise modelling).

In [START_REF] Wörz | A new 3D parametric intensity model for accurate segmentation and quantification of human vessels[END_REF] a cylindrical parametric intensity model is directly fit to the image intensities through an incremental process based on a Kalman filter for estimating the vessels radii, while in [QBD + 09], local neighborhood intensities are considered in a spherical polar coordinate system in order to capture the common properties for the different types of vascular points.

A natural integration into this kind of models is background description [SMS + 07, TdTF + 07]. See Fig. 5 for the illustration of the method. While being simple, intensity models are highly dependent on the nature of images. Therefore, they have to be tuned for all kinds of circumstances, like artefacts or other image distortions, as well as to compensate for image variability.

Geometry models

The assumption that vessels are elongated thin objects, globally similar to tubes has been used for the design of several geometric models, such as generalized cylinders, superellipsoids or Gaussian lines [KMA + 00, TBD + 01, TdTF + 07, SHS04, FHP08]. Bar-like profiles [KGSD95, BRT03, WR04, WGR09] account for the fact that bigger vessels are flatter than perfect cylinders. The result is the maximum of the convolution of linear filter with an image at multiple angles and it can distinguish between non-and tubular objects.

Based on second-order derivatives, several models incorporating geometrical properties have been developed. In [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], an ideal cylinder is proposed in order to enhance vessels within a measure called vesselness, while in [SNA + 97] a more general model incorporates elliptical shapes, also enriched in [START_REF] Lin | Stenosis detection using a new shape space for second order 3d-variations[END_REF] which also includes a conus-like descriptor, which can be used to detect vessel stenosis (see Fig. 6).

The bifurcation issue has also been considered, for instance in [START_REF] Antiga | Robust and objective decomposition and mapping of bifurcating vessels[END_REF] where a bifurcation models is proposed and optimised based on vessel centerline information. An alternative solution, described in [START_REF] Agam | Probabilistic modeling based vessel enhancement in thoracic CT scans[END_REF] is based on optimized probabilistic models derived from eigenvalue analysis of the structure tensor.

Geometry models are powerful tools for describing vessels and aiding to their further extraction within tracking schemes or by deformation. However, these methods, assume image regularities that are present in high-quality images but not necessarily in noisier ones, nor in pathological cases. Furthermore, they often require careful parameter tuning, which may change from one data-set to the next. Probabilistic and/or statistical approaches can contribute to decision-making whether pixel belongs or not to a vascular structure.

Deformable models

Deformable models aim at fitting a geometric hypersurface (e.g., a 2D surface in a 3D image), by moving it and modifying its shape from an initial model, under the guidance of several (generally antagonist) forces: external ("data-driven") ones, related to the image content, and internal ("model-driven") ones, devoted to preserve correct geometry properties (e.g., regularity). Such models have been intensively used in the field of image analysis [START_REF] Singh | Deformable Models in Medical Image Analysis[END_REF], especially due to their advantages: arbitrary shape representation, topological adaptivity, sub-pixel precision, etc.

Among the most classical methods, snakes (often used in 2D in order to segment vessel cross-sections), have been considered, e.g., in [START_REF] Mcinerney | Medical image segmentation using topologically adaptable surfaces[END_REF], or in [START_REF] Hernandez | Three-dimensional segmentation of brain aneurysms in CTA using non-parametric region-based information and implicit deformable models: Method and evaluation[END_REF], where two (1D and 2D) snakes are used for both segmentation and stenosis quantification.

Level-sets constitute another classical type of deformable models, and rely on an Eulerian version of contour evolution with partial derivative equations. The contour is integrated as the zero-level of a higher dimension function (level-set). In [LFG + 01], an original level-set based scheme proposes to deform an initial boundary estimate toward the vascular structures in the image using a codimension-two regularization force, based on the vessel centerlines instead of the vessel surface (see Fig. 7 for an illustration of the method).

In [START_REF] Yan | MRA image segmentation with capillary active contour[END_REF] a similar method called capillary active contours is introduced which adapts the evolving surface into very thin branches of blood vessels and obtains more accurate segmentation in comparison with [LFG + 01] (see Fig. 8).

An other level-set based method [MVvL + 06] proposes to estimates the background and vessel intensity distributions based on the intensity histogram, to more efficiently steer the level set onto the vessel boundaries. Several efforts have been conducted to improve deformable models in the quite specific case of elongated structures. In this context, [START_REF] Vasilevskiy | Flux maximizing geometric flows[END_REF] proposes to use flux maximization as an alternative curvature-based regularization which make surface normals evolve according to gradient vector field. The key idea was to evolve a curve or a surface under constraints by incorporating not only the magnitude but also the direction of an appropriate vector field. [START_REF] Descoteaux | Geometric flows for segmenting vasculature in MRI: Theory and validation[END_REF] has used in the similar manner Hessian values flow.

In [START_REF] Law | Weighted local variance based edge detection and its application to vascular segmentation in magnetic resonance angiography[END_REF], local variances are measured with first-order derivatives and are propagated according to their strengths and directions with an optimally oriented flux reporting more accurate and stable responses and higher robustness to disturbances from adjacent structures in comparison with Hessian-based measures.

The major advantage of deformable models methods is that they are sensitive to weak edges and robust to noisy structures. However, the intensity variation inside vascular structures can generate significant intensity gradient with this undesired discontinuity stopping the contour evolution at these regions. Due to this local minima, the initial forces should be described with such precision that the final object borders are not far from the initial ones. While the evolution of the deformation can be a costly process.

thin objects filtering and segmentation methods But, by integrating vessel features and forces in powerful optimization schemes help overcoming these problems.

Statistical approaches

Vessel segmentation based on statistical approaches generally relies on specific assumptions related to the intensity distribution of the vascular/non-vascular signals in MRA data (only very few statistical methods have been devoted to CTA, see, e.g., [START_REF] Florin | Particle filters, a quasi-Monte-Carlo-solution for segmentation of coronaries[END_REF], which proposes a particle-filtering strategy for coronary arteries segmentation), and especially physical models of blood flow. If the number and the nature of these distributions is supposed to be known correctly, it is then possible to determine their respective parameters (and in particular the mean intensity characterising the associated structures), via a standard Expectation-Maximisation (EM) technique [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

In (TOF or PC) MRA, two or three distributions are generally considered, for the blood, and the other anatomical structures and the background, respectively. They led, in particular to the definition of Gaussian-Gaussian-uniform [START_REF] Wilson | Determining X-ray projections for coil treatments of intracranial aneurysms[END_REF] and normal-Rayleigh-2×normal [SHFHM06] mixtures for TOF MRA, and Maxwell-Gaussian [CNSB01], Maxwell-Gaussian-uniform [START_REF] Chung | Fusing speed and phase information for vascular segmentation of phase contrast MR angiograms[END_REF] mixtures for PC-MRA. In [START_REF] Chung | Vascular segmentation of phase contrast magnetic resonance angiograms based on statistical mixture modeling and local phase coherence[END_REF], a hybrid model, enables to choose between these two kinds of mixtures. Alternatively, to these "constrained" mixture choices, it has been proposed in [START_REF] El-Baz | Automatic cerebrovascular segmentation by accurate probabilistic modeling of TOF-MRA images[END_REF] to consider a linear combination of discrete Gaussians with alternate signs, involved in a modified EM, which enables to adaptively deal with both laminar and turbulent (pathological) blood flow [EBFG + 06].

In the primarily considered strategies, the determination of the vascular intensity led to a straightforward segmentation by thresholding of the image (sometimes enriched by a hierarchical analysis of the image by octree decomposition [START_REF] Wilson | Determining X-ray projections for coil treatments of intracranial aneurysms[END_REF]). From an algorithmic point of view, segmentation improvements were also performed in considering of spatial information (i.e., statistical dependence) between neighbour voxels, by integrating Markov random fields (MRF) [GHP + 95] in a post-classification correction step [START_REF] Sabry Hassouna | Cerebrovascular segmentation from TOF using stochastic models[END_REF]. In other works, speed and phase information provided by PC-MRA were fused and involved in a maximum a posteriori-MRF framework to enhance vessel segmentation [START_REF] Chung | Fusing speed and phase information for vascular segmentation of phase contrast MR angiograms[END_REF][START_REF] Chung | Vascular segmentation of phase contrast magnetic resonance angiograms based on statistical mixture modeling and local phase coherence[END_REF].

Statistical methods globally inherit from the strengths and weaknesses of the EM algorithm. First, they generally require to establish hypotheses on the signal distribution. Moreover, they involve several parameters, for instance, weight, mean and standard deviation, of the distributions. The initialisation of the segmentation process then requires a special attention. Indeed, the convergence may possibly depend on the quality of the initial distribution settings (sometimes automatically determined based on heuristic rules [START_REF] Wilson | Determining X-ray projections for coil treatments of intracranial aneurysms[END_REF][START_REF] Chung | Fusing speed and phase information for vascular segmentation of phase contrast MR angiograms[END_REF]). As for any optimisation strategy, the termination also requires to decide whether the process has correctly converged or not (which is sometimes empirically determined, for instance by a maximal number of iterations [START_REF] Wilson | Determining X-ray projections for coil treatments of intracranial aneurysms[END_REF]). Finally, since the segmentation process is strongly based on photometric properties (the results often consist of global or local thresholdings), higher-level knowledge such as geometric assumptions are hardly considered, and require postprocessing steps based on a statistical framework [START_REF] Sabry Hassouna | Cerebrovascular segmentation from TOF using stochastic models[END_REF], or, more efficiently the collaboration of alternative image processing techniques (see examples in Section 2.3).

Minimal path techniques

Based on extremal intensity and connectedness criteria, the detection of a vessel segment (or more precisely of its medial axis) can be expressed as the determination of a minimal cost path in a weighted graph modeling voxels, their neighbourhood relations and their intensity.

Vessel segmentation based on such strategies can rely on standard minimal path finding techniques [START_REF] Dijkstra | A note on two problems in connection with graphs[END_REF] (i.e., on "global" minimisation strategies, while methods categorised in the next Tracking section will rely on "local" (step by step) minimisation strategies). This is, for instance, the case in [START_REF] Olabarriaga | Minimum cost path algorithm for coronary artery central axis tracking in CT images[END_REF].

Alternatively to classic path-finding methods, fast-marching strategies [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] have been considered. These methods are both related to the level-sets methodology (see the Deformable models paragraph) and minimal path-finding ones (they remain, in particuler, consistent with the continuous formulation of the minimal-path research). By opposition to fully discrete path-finding, they enable in particular to determine paths with a sub-voxel accuracy [START_REF] Avants | An adaptive minimal path generation technique for vessel tracking in CTA/CE-MRA volume images[END_REF].

The methods based on path-finding are globally well-fitted for the detection of vessel medial axes, especially in the case of small vessels which justifies in particular their frequent use in coronary detection. (For larger vessels, the optimal path may diverge from the medial axis, leading to eccentric results, this issue then requiring a specific care [START_REF] Li | Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centerlines[END_REF].) However, it has to be noticed that efforts have also been conducted to develop segmentation methods enabling to extract both vessel axes and vessel walls [BC09, LY07], which express the whole vascular volume segmentation as the minimisation of a path in a space enriched with a supplementary "scale" dimension corresponding to the vessel radius.

Despite attempts for segmenting whole vascular trees [START_REF] Yim | Gray-scale skeletonization of small vessels in magnetic resonance angiography[END_REF], such methods generally remain devoted to the segmentation of vessel segments, thus requiring to interactively provide at least an initial point, and also a final one [OBN03b, WFV + 02]. In this case of use, they may be robust to noise, and signal decrease (or shirt signal loss) along the vessel, especially in case of stenoses. These methods being based on monotonic and/or finite algorithmic processes, their termination is garanteed, and their theoretical algorithmic cost is generally low. Practically, the computational cost may however be high, and in this context, the proposal of initial and final points can potentially enable its reduction by computing paths from both points simultaneously [START_REF] Olabarriaga | Minimum cost path algorithm for coronary artery central axis tracking in CT images[END_REF].

thin objects filtering and segmentation methods

Tracking

The basic idea of the tracking methods is in following the object of interest along its principal direction. It starts from given seeds and direction and it is an iterative process of further steps prediction and correction. Usually, vessel tracking is combined with finding its centerline. A centerline can be estimated by looking at the 2D orthogonal plane in the principal vessel direction. Both centerline and direction can be estimated with geometrical and intensity model-based methods.

In [START_REF] Boldak | An improved model-based vessel tracking algorithm with application to computed tomography angiography[END_REF] tracking was based on geometrical moments for model fitting and local directions for vessel tracking, as well as enhancing along the vessel. If one bifurcation is detected, some seed points are extracted which are saved in a list to be further taken as initial point of the tracking process. The tracking scheme is illustrated in Fig. 9.

Many tracking techniques follow one branch at a time, relying on manual reseeding to extract a complete tree. To automatize bifurcation handling, some authors have proposed to perform the segmentation of the lumen locally and rely on topological knowledge.

In [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF] segmentation of each vessel begins from a seed point tracking the image intensity ridge representing the vessel skeleton in 3D based on Hessian values, using dynamic scaling and automated calculation of vessel width at each skeleton point. Since the method extracts only objects approximately circular in cross-section, it resists inclusion of non-vessel objects. An advantage of the approach is that it is capable of defining even very small vessels.

With geometric models the previous direction can be, in the simplest way, directly reused in the next step [TdTF + 07]. Please see Fig. 10 for a schematic illustration.

Otherwise, directions can be filtered, like in [FDCR01, WNV00] a set of sliding volumes is used to detect and follow branching vessels by automatical iterative tracking of the vessel centreline modelled by second-order B-spline.

Often the next tracking point is corrected according to the best fitness of a sphere [COH05, HHOPJ + 06] or medialness [START_REF] Wink | Fast delineation and visualization of vessels in 3-D angiographic images[END_REF]. Algorithm by [START_REF] Carrillo | Extraction of 3d vascular tree skeletons based on the analysis of connected components evolution[END_REF] applies local adaptive threshold and tracks the branches detecting the bifurcations by analyzing the binary connected components on the surface of the sphere.

Carrillo [START_REF] Carrillo | Recursive tracking of vascular tree axes in 3D medical images[END_REF] has reused the sliding spherical volume along the vessels. While this process is guided by the minimization of a new criterion consistently combining intensity variances with the spatial inertia moments for corresponding clusters in classification process. Not relying on contrasts assumptions (as clustering supposes intra-class intensitities to be similar), a strong a priori on the vascular shape is introduced via an adaptive geometrical model. This operation enhances cylindrical structures in the image.

Correction can be performed with 2D active contours with an accurate segmentation of the cross-section [START_REF] Lee | Automatic segmentation of 3d micro-ct coronary vascular images[END_REF]. Cross-sectional recentering can be done with local optimization of 3D models as in [TdTF + 07, WR08].

Kalman filtering was used as an other prediction and correction scheme in [WR04, WR08], which is theoretically very efficient under Gaussianity and linearity assumptions.

The integration of models generally improves robustness. A problem that may occur with such iterative approaches is that if the method proceeds in the wrong direction at one point, there is no mechanism to correct for this. The robust termination criteria should also be offered in order to achieve complete tree segmentation and not include false positives.

Furthermore, these approaches often depend on many parameters and can become complex due to the estimation of the next direction vector, making them more prone to failure especially for strongly curved structures, non-tubular structures (e.g., bifurcations), or vessels with pathological conditions.

One more way to further improve robustness is to use multi-hypotheses frameworks such as stochastic particle filters in order to deal with non-linear, non-Gaussian processes (see 2.2.6 for more details).

The fact that tracking, like some region-growing and deformable models, focuses on the object of interest, accelerates the processing.

thin objects filtering and segmentation methods

Morphological methods

In broad terms, the morphological filters that are of interest for vessel filtering and segmentation include openings, closings and thinnings. Openings filter out small bright objects over a dark background while closings remove small dark objects over a light background. They are dual of each other in the sense that a closing performs like an opening if the contrast is inverted. We also need to know that some composition of operators are interesting. Namely, a supremum of openings is still an opening, and an infimum of closings is still a closing. Finally, thinnings are like openings, except they are not increasing, meaning they don't maintain order.

structuring element based approaches Filtering thin objects with morphology can be achieved using appropriate structuring elements. Typically, thin structuring elements include segments and paths, combined over families. To account for arbitrary orientation, one can use families of oriented segments and compute a supremum of openings or an infimum of closings. This is described in [START_REF] Soille | Directional morphological filtering[END_REF]. To account for noise or disconnection, families of incomplete segments can be used instead, yielding so-called rank-max openings, which are just as efficient and also described in the same reference.

Paths are elongated structuring elements, but that are not necessarily locally straight. Even though the size of families of paths grow exponentially with their length, there exists a recursive decomposition that makes the use of such families tractable [START_REF] Heijmans | Path openings and closings[END_REF]. As with segments, it is useful to account for some discontinuities using so-called incomplete paths. As with segments, there exists an efficient implementation [START_REF] Talbot | Efficient complete and incomplete paths openings and closings[END_REF]. In fine, path and segment operations are comparable in speed.

In [START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation[END_REF], authors use segment-based operators and curvature evaluation for 2D blood vessel segmentation. In [START_REF] Tankyevych | Direction-adaptive grey-level morphology. Application to 3D vascular brain imaging[END_REF], the author shows that path and segment morphological operators significantly outperform linear and steerable filters for the segmentation of thin 2D structures, even in the presence of heavy noise. Segment-based 3D morphological operators were used in [START_REF] Dokládal | 3-d extraction of fibres from microtomographic images of fibre-reinforced composite materials[END_REF] with good performances. Paths operators have been extended to 3D in [START_REF] Hendriks | Constrained and dimensionality-independent path openings[END_REF], and show to outperform all other morphological filters for thin object segmentation in 3D, both for efficiency and performance. connected component approaches Connected operators are also supremum of openings or infimum of closings, but use families of structuring elements that are so large that it makes little sense to present them in this way. Instead, the concept of connectivity is used [Vin93a, SS93, SS95]. The simplest of those is the area opening or closing. Informally, the area opening suppresses objects that are smaller in area than a given size λ. It extends readily to arbitrary lattices, and corresponds to a supremum of openings with a very large family of structuring elements: all the connected sets that have an area smaller than λ. In the continuum this family is not countable, but in the discrete case it is still very large. Fortunately it is not implemented in this way. A very efficient way to implement this operator is via the component tree [SG00, MW02, NC06]. Many operators exist, that were derived from the area opening, and are well suited to thin objects filtering and segmentation. For instance, in [WME + 07], connected operators were used to segment blood vessels in the eye fundus and detect microaneurisms. In [START_REF] Michael | Shape preserving filament enhancement filtering[END_REF], a scale-independent elongation criterion was introduced to find vascular structures.

Some useful connected operators are thinnings rather than openings, as they make it possible to use more complex criteria for object selection, for instance using elongation measures, that are not necessarily increasing. In [START_REF] Westenberg | Volumetric attribute filtering and interactive visualization using the max-tree representation[END_REF], various connected criteria are uses for elongated object 3D volume rendering, including blood vessels.

hit-or-miss transforms Hit-or-miss transforms repeatedly use pairs of structuring elements (SEs) to select objects of interest, rather than single SEs. In [NPR07, BRB + 10], authors used such operators for 3D vessel segmentation, including brain and heart vessels. See Fig. 11 for an example of a segmentation with hit-or-miss transform according to [START_REF] Naegel | Grey-level hit-or-miss transforms -Part II: Application to angiographic image processing[END_REF].

combinations of methods

Despite the huge amount of methodological contributions dedicated to 3D vessel segmentation, proposed during the last twenty years, the results provided by such segmentation methods generally remain perfectible.

The handling of under-segmentation (especially in the case of small vessels, whose size is close to the image resolution, of signal decrease, or of partial volume effect) and over-segmentation (especially in the case of neighbouring with other anatomical structures, or of high intensity artefacts), the robustness to image degradations (low signal-to-noise ratio), the ergonomy (automation, or easy interaction), the low computational cost, the guarantee of termination and convergence, accuracy of the result (for instance, the ability to provide results at a higher resolution than the image one) are desirable properties for such methods. Unfortunately, none is generally exempt from drawbacks, even in the frequent (and justified) case where the method is devoted to a quite specific task, vascular structure, and/or image modality.

As nearly all the main strategies of image processing have been -not fully satisfactorilyinvestigated to propose solutions to this issue, a reasonable trend during the last years has consisted in designing hybrid segmentation methods obtained by crossing methodologies. An alternative way to overcome this issue is to inject more guiding knowledge in the segmentation processes, which justifies -among other reasons-the generation of anatomical vascular models, as discussed in Section 2.2.4. These strategies aim, in particular, at taking advantage of (distinct and complementary) advantages of different segmentation techniques.

A synthetic overview of such hybrid methods is proposed hereafter. A more detailed example of such a method is then presented.

Principal strategies

Vessel segmentation hybrid methods present a range of possible new paths, through which a number of discussed above challenges can be overcome.

derivatives-based

In [START_REF] Chen | Quantifying 3-D vascular structures in MRA images using hybrid pde and geometric deformable models[END_REF] multi-scale filtering based on the Hessian matrix is used to effectively enhance vessel structures with various diameters. The level set method is then applied to automatically segment vessels enhanced by the filtering with a speed function derived from enhanced MRA images. The segmented vessel surface is triangulated using 3D Delaunay triangulation and the resulting surface is used as a parametric deformable model. Energy minimization is then performed within a variational setting with a first-order internal energy; the external energy is derived from 3-D image gradients.

Works by Bemmel and Descoteux [BSVN03, DCS08] use Frangi's [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] vesselness combined with the maximum geometric flow by Vassilevskiy [START_REF] Vasilevskiy | Flux maximizing geometric flows[END_REF].

In[GLM + 08] a level-set-based geometric regularization method is proposed which has the ability to estimate the local orientation of the evolving front and utilize it as shape induced information for anisotropic propagation. It is shown that preserving anisotropic fronts can improve elongations of the extracted structures, while minimizing the risk of leakage. For an evolving front using its shape-offset level-set representation, a novel energy functional is defined.

Lesage [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF] has integrated Flux formulated as the measure of inward gradient flux through a local circular cross-section into a tracking strategy reducing false positive responses in situations such as step-edges.

One more hybrid method is by Frangi [FNH + 99] used vessel filtering [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] with two B-spline snakes to first detect the vessel axis (3D contour) and then the vessel wall (3D surface).

In [START_REF] Lee | Automatic segmentation of 3d micro-ct coronary vascular images[END_REF] there has been used a tracking strategy similar to [START_REF] Flasque | Acquisition, segmentation and tracking of the cerebral vascular tree on 3D magnetic resonance angiography images[END_REF] in combination with active contours while obtaining initial seeds with the Hessian matrix values as in [START_REF] Wink | Fast delineation and visualization of vessels in 3-D angiographic images[END_REF] tested on a rat's heart micro-CT and resulting in a full heart vasculature reconstruction. A limitation of this approach is that during tracking, the vessel axis is modeled by a sequence of piecewise linear elements as opposed to a b-spline or polynomial Hermite finite elements, which could offer continuity in gradient over the whole segment.

In [START_REF] Friman | Template-based multiple hypotheses tracking of small vessels[END_REF] a multiple hypothesis tracking (MHT) framework for 3D vessel segmentation was used with new Gaussian vessel profile and a statistically motivated criterion for assessing the model fit. This makes it a deterministic alternative to the recently published Bayesian particle filters [SMS + 07, FPW05]. Whereas the particle filters stochastically sample their way forward, the MHT samples its way with systematic sweeps. Where the particle filter furnishes a probability distribution of possible vessel paths, the MHT uses the tracking search tree. An advantage of the MHT over the particle filter is its computational complexity. Significant gain in robustness is obtained with the increase of the depth of the search tree.

Santamaria-Pang [SPCSK07] used the Hessian matrix together with Support Vector Machines (learning method) in order to segment neurone dendrites managing to achieve good result even for irregular shapes and image specificities. Blondel [START_REF] Blondel | 3D tomographic reconstruction of coronary arteries using a precomputed 4D motion field[END_REF] has proposed a method to compute 3D tomographic reconstructions of coronary arteries moving under cardiac contraction. It is based on a vessel detectors by Sato [SNA + 97] and Krissian [KMA + 00] where a binary mask of detected vessels in 2D is created. Then, with the help of morphological operators the artificial subtraction of the vessels is performed from the original image.

Manniesing [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] combined a smooth vessel filter based on geometrical analysis of the Hessian eigensystem with a non-linear anisotropic diffusion scheme. The amount and orientation of diffusion depend on the local vessel-likeliness. Vessel enhancing diffusion (Vessel Enhancing Diffusion (VED)) is show in a diameter study on phantom data to have least affects on the accuracy of diameter measurements compared to some other classical methods. It is shown to improve level-set based segmentation of the smaller vessels.

In [GAF + 09] they use Support Vector Machines (Support-vector machine (SVM)) that operate on the Hessian and quadrature eigenvalues for discriminating between filament and non-filament voxels. They show that using fourth order steerable features allow to detect filaments more accurately than second-order methods, notably to Frangi's vesselness model. See Fig. 13 for the result curves. Interestingly, despite expectations on SVM learning performance, in the best accuracy positions (somewhere in the middle) of these ROC curves -where the true positive values are around 80% -the named methods perform similarly.

Multiscale morphology with Gabor wavelets (providing vessel size and direction) filters were used in [START_REF] Sun | Enhancement of vascular angiogram by multiscale morphology[END_REF][START_REF] Sun | Morphological enhancement of vascular angiogram with multiscale detected by gabor filters[END_REF]. Taking advantage of the Gabor wavelet that is capable of tuning to specific frequencies allowing estimation of the vessel dimension.

statistics-based Hernandez [START_REF] Hernandez | Non-parametric geodesic active regions: Method and evaluation for cerebral aneurysms segmentation in 3DRA and CTA[END_REF] proposed an automatic method for segmentation of cerebral vascular structures applied to the segmentation of cerebral aneurysm. It uses a geometric deformable model that evolves minimizing an energy functional that incorporates statistical region-based information estimated in a high-order multiscale feature space with a non-parametric model. Gooya [GLM + 07] used region statistical measures and orientation information from ramp-like edges, which are fused within an energy minimization scheme that is based on a new interpretation of edge concept. A region driven abjection term simulates the edge strength effect is directly obtained from this minimization strategy.

conclusion and discussion

The above-described methods present a wide range of choice for a particular problem solving. The choice depends on the application and its constraints. For example, the most important features to decide upon are: (1) if a method should be automated or not; (2) if the efficiency is of high priority; (2) the precision and other supplementary information availability.

On the other hand, most of the methods have the same weaknesses in automatically distinguishing between vessel bifurcation and vessel fusion generated as an artefact in the image, as well as in the ability to separate between adjacent structures. Another weakness of current automated vascular reconstruction is in removing some features, such as presence of the objects that are not supposed to be included in the simulation or, otherwise, in dealing with incomplete data when a part of the vessel is not seen.

One of the ways to try to overcome limitations of some methods is by combining different complementary methods. This has been shown a way to achieve the best performance. The choice of each step depends on the computational constraints and on requirements of the next step.

Another way for achieving better results in vascular image analysis is by adaptivity. Tracking methods are good examples of such strategies. Their iterative nature permits to exploit the image data while evolving the analysis process. However, such schemes has also limitations in cases where they might be stopped by an obstacle without moving any further.

From such a variety of methods for thin objects detections, filtering and segmentation, we are interested in methods that can be automated, that are multi-scale, are able to distinguish thin objects and other structures and detect their directions. Moreover, we are interested in combining the scale-space methods that have an approach from localto-global with the purely local though well-defined but less exploited mathematical morphology methods. Our proposition with such a method is described further.

...the object becomes virtual, the image actual...

-Gilles Deleuze and Claire Parnet

For the screening and diagnostic of cerebral vascular diseases, as well as for planning incisions to a lesion, there is a strong need for efficient visualization. However, displaying large 3D data on 2D screens is problematic. The display of such data can benefit from the knowledge of quantitative shape features such as centrelines of vessels. These can be obtained from a curvilinear skeleton representation. If connectivity and topology are preserved, and if geometrical constraints such as smoothness and centering are satisfied, these methods serve not only for visualization, but as well for analysis of the topology, and estimations of length, width, orientation, curvature, and torsion.

introduction

The global external form of an object, or morphometry, can represent a robust quantification of complex structures. However, such representation treats the initial structure as one object, missing vessel-specific properties such as bifurcations, individual segments, or parent-child relations. Various measures, such as of topology, geometrical changes in diameter, length, branch point density or tortuousness are important in vessels disease studies. The vessel centrelines can serve as a basis for its description, and especially with a possibility of quantitative measures. Their main feature and advantage is that they process data locally while preserving global properties.

Extraction of the vessel medial line is useful in several vascular applications: evaluation of pathology severity and for measurements extraction such as lumen diameter, wall thickness, etc. These applications require axes to be aligned with the local orientation of the vessel.

In the following sections, the skeletonisation together with associated computational methods are reviewed, decomposition schemes and other technical details will be presented. In our work, we used only discrete skeletonisation schemes but we also present several works on vessel skeletonisation in continuous space.

Furthermore, vascular tree extraction is a first step toward an automated search of possible pathologies. It also helps physicians to pay attention at the vascular structures rather than aspect, noise, etc and therefore makes such huge datasets analysis more efficient. The vascular tree construction methods are presented in this chapter.

Since they were first introduced, curve-skeletons have been useful in many areas, e. g.image processing, visualization, animation, etc.

Efficient interaction between human observers and imaging systems is more important today than when screen-film systems dominated. Advanced imaging technologies have resulted in thousands of images per case that the radiologist must handle, which has amounted in a significant increase in the time required to view a case. Surveys state that clinicians work over-load has grown in comparison with the past [BHS + 06]. Image reconstruction and analysis, computer-aided detection and diagnosis, multi-modality comparison and integration, and a host of other software tools are needed to help clinicians make sense of the image data and render the best diagnostic decision.

In the clinical context, up till now, for displaying large 3D data on 2D screens, two major graphical approaches have been largely employed: maximum intensity projection (MIP) or similar and surface rendering. Volume rendering even on specialized hardware is very resource intensive. We also present in this chapter how the skeletons can aid visualization of large vascular data.

short presentation of digital topology

When considering the discrete space Z 3 , each point p in this space is called a voxel defined by its three coordinates (p x , p y , p z ). A voxel can be considered as a cube containing of 6 faces, 12 edges and 8 corners. Two voxels p and q ∈ Z 3 are 6-adjacent if they have a common face; 18-adjacent is they have a common face or edge, and 26-adjacent if they have a common edge, face or corner. The set of n-adjacent voxels (n = {6, 18, 26}), can be accordingly called 6-, 18or 26-neighbourhood of p, denoted by N n (p). The illustration of such neighbourhoods is presented in Figure 14.

Here, we use the definition of the following 3D neighbourhoods for a point x in an arbitrarily small neighbourhood V: 6-neighbourhood:

N 6 (x) = V 1 1 (x) and N * 6 (x) = N 6 (x) \ {x} 26-neighbourhood: N 26 (x) = V ∞ 1 (x) and N * 26 (x) = N 26 (x) \ {x} 18-neighbourhood: N 18 (x) = V 2 1 (x) ∩ V 1 ∞ (x)
and N * 18 (x) = N 18 (x) \ {x} An n-path is a sequence of voxels p 1 , ...p j with p i n-adjacent to p i+1 . While the adjancency is a symmetric and irreflexive relation. An n-connected component is then a set of voxels such that any two such voxels are connectd by an n-path included in that component.

A 3D binary digital image is then described as I = (Z 3 , m, n, O), where O ⊆ Z 3 is the set of voxels representing the object, while Z 3 \ O denotes the background voxels. The pair (m, n) specifies the object and background connectivity respectively. Often, with the aim to avoid topological paradoxes of objects being both connected and disconnected, different values must be chosen for m and n. Popular choices are (26, 6) and (18, 6).

A cavity is a background space surrounded by the object. A tunnel is a hole in the centre of torus.

Points classification in 3D

We use the same classification of 3D points as in [BM92]:

1. Each point is labeled with a topological type using the computation of two connected components numbers in a small neighbourhood.

2. Because some points (junction points) are not detected with the twonumbers, a less local approach is used for extracting them.

Let us consider an object X in the real space R 3 , let x ∈ X, and let V(X) be an arbitrarily small neighbourhood of x. Let us consider the number C R 3 which is the number of connected components in X ∩ (V(x) \ {x}) adjacent to x. These numbers may be used as topological descriptors of x.

We use this following number of connected components:

C = NC a [X ∩ N * 26 (x)] is the number of 26-connected components of X ∩ N * 26 (x) 26-adjacent to x.
All points in the 26-neighbourhood are 26-adjacent to x, therefore C is the number of 26-connected components of X ∩ N * 26 (x). It is not necessary to check the adjacency to x. The we obtain a classification of each point of the object using this number C an in Table 1. However, this classification depends only on the 26-neighbourhood of each point and some junction points belonging to a set of junction points which is not of unit-width are not detected. Then we use the following proposition for curves according to [START_REF] Bertrand | A new topological classification of points in 3d images[END_REF]:

we only need to count the number of neighbours of each curve point (type D), if this number is greater than two, the point is a missed curves junction point (type E).

This way, the connected component of the background is counted as a 18-neighbourhood N * 18 (x). Classification done with such a smaller neighbourhood ensures that we have finer details of the object.

skeletonization

The skeleton is a concept aimed at reducing the quantity of information of an object while retaining the same topological information as the object they were extracted from (in 3D, same number of junction points, tunnels, cavities and connected components). In the continuous framework, the ultimate skeleton of an n-dimensional object in an n-dimensional space is at most an (n -1)-dimensional object. Intuitively, in 2D, a skeleton of a shape X is a curvilinear representation of X, having the same topology as X. In 3D, the ultimate skeleton of an object will be at most a bi-dimensional object (no volume, only surfaces, curves, and points).

This property is very interesting in object analysis, especially in large image data. In general, the computer process of skeletonisation is in some ways similar to the processing the human brain performs when seeing an object -decomposing it into associated parts. In the mid/late 50s, it was realized that computers may be able to recognize patterns, and that reducing a thick pattern in a curvilinear object was a good strategy in reducing the amount of information needed to be processed by the machine. The skeletons were originally defined by Blum ([Blu62]) as the meeting point of the flames of a grassfire (see Fig. 15). Imagine an object as if it were a prairie grass paddock; if one sets on fire the contour of the object, then the meeting points of the firefronts would constitute the skeleton of the object. In the continuous framework, this definition is equivalent to saying that the skeleton is the set of points which are centres of maximal balls (balls included in the object, and not strictly included in any other ball) ( [START_REF] Calabi | A study of the skeleton of plane figures[END_REF]).

Figure 16: Illustration of a maximal ball in the Euclidean case. The ball B max is maximal because no other ball can contain it and simultaneously be included in the set. In contrast, the ball B 2 is not maximal (source: [START_REF]Mathematical morphology: from theory to applications[END_REF])

In 1969, Hilditch gave four properties that a skeleton in a bidimensional space should possess ( [START_REF] Hilditch | Linear skeletons from square cupboards[END_REF]). Adapted to the general case of n-dimensional skeletons, these properties are that a skeleton should be : 1) homotopic to the original object, 2) thin (in an n-dimensional space, the skeleton of an object should be at most (n -1)dimensional), 3) centered in the original object, and 4) skeletonising a skeleton should not change anything . This last property is only valid when dealing with ultimate skeletons.

In the continuous framework, the set of centres of maximal balls, called the medial axis, hold these properties. An illustration of medial axis can be viewed in Figure 17. In the discrete framework Z n , the discrete medial axis does not always hold two of these properties: it is not always homotopic to the original object, and it is not always thin. Shown on Fig. 18 is a counter-example. In Fig. 18(a) is a figure of a bear, which is homotopic to a disk. However, as shown on Fig. 18(b), its discrete medial axis with a 4-connected unit ball is neither connected not thin. On Fig. 18(c), its ultimate skeleton is not reduced to a single point, as it should be.

Here, we consider only the discrete skeletonization schemes, however, there are a few works on continuous skeletonization methods, meaning work where the formulation of the skeletonization process is purely continuous. In [START_REF] Deschamps | Fast extraction of minimal paths in 3D images and applications to virtual endoscopy[END_REF] centerline extraction of segmented tubular objects is accomplished by evolving monotonic fronts, where the cost function is a distance function from the edges of the binary object of interest. Similarly, in [START_REF] Hassouna | Accurate tracking of monotonically advancing fronts[END_REF] centerline extraction is proposed from monotonic front evolution, where the centerline follows the maximal distance from the boundary of the binary object. In recent work [START_REF] Bouix | Flux driven automatic centerline extraction[END_REF] the average outward flux is used through a Jordan curve.

In this case, the gradient vector field of the Euclidean distance to the boundary was exploited to approximate the centerlines.

Many methods have been developed to perform skeletonisation of a discrete object, an overview of these methods is proposed in the following sections. According to [START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF], discrete skeletons can be computed using four types of methods: Voronoibased transformations [BA92, NSK + 97], distance-based transformations [BNB99, TM01], general-field methods [AC97, RT02] and thinning. Here, we focus only on thinning methods.

Moreover, when working with computers, input objects and output skeletons are digital objects, which means that they are made of voxels. In such case, decomposing a skeleton into curves and surfaces is not a trivial problem, and we will have a look at the various methods that were proposed to solve it.

Even though skeletons have the same topology as the objects they were computed from, they do not always preserve their visual aspect. The following sections will be about how it is possible to obtain skeletons preserving the visual aspect of an object, and how it is possible to decompose them into basic parts (curves and surfaces).

Simple points in 2D and 3D

In the world of thinning, the "currency" is the simple point [START_REF] George | Three-dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF]. Intuitively, a point is simple if it can be removed from an object without changing its topology. In the discrete framework, the topology of an object depends on its connectedness; for this reason, when considering a k-connected object, we will talk about k-simple points. The notion of simple point is central for homotopic thinning in the digital framework: a skeleton is obtained by removing iteratively simple points from an object.

According to [START_REF] Lam | Thinning methodologies -a comprehensive survey[END_REF], in the 60s, simple points were characterized based on connectedness: a point p is k-simple for an object X if the removal does not change the number of k-connected components of X nor the number of k-connected components of X (where k is the usual connectedness for X when k is chosen as connectedness for X) [START_REF] Duda | Graphical data processing research study and experimental investigation[END_REF][START_REF] Marcel | Hexagonal parallel pattern transformations[END_REF]. But this definition complicates computational framework: indeed, in order to test if a single point is simple, this definition requires scanning the whole object in order to enumerate its connected components. Fortunately, local characterization of deletable points in 2D began to appear in the mid 60s, based on crossing numbers [START_REF] Rutovitz | Pattern recognition[END_REF][START_REF] Hilditch | Linear skeletons from square cupboards[END_REF], connectivity numbers [START_REF] Yokoi | An analysis of topological properties of digitized binary pictures using local features[END_REF] and simplicity [START_REF] Rosenfeld | Connectivity in digital pictures[END_REF].

All these works established that, in order to decide whether a point is deletable or not, it is only necessary to look at the configuration of the point's neighbourhood (no need to count the number of connected components of the whole object). Consequently, in 2D, deciding if a point is simple can be performed in constant time.

Proposition 3.3.1 [START_REF] Rosenfeld | Digital topology[END_REF] Let X ⊂ Z 2 , and p ∈ X. We denote by N(p) the set of 8neighbours of p.

If X ∩ N(p) has the same number of components (in the sense of X) as X ∩ (N(p) ∪ {p}), then p is simple for X.

In 3D, the removal of a point may not only change the number of connected components of the object, but may also change the number of tunnels of the object. As in the 2D case, 3D simple points can be locally characterized [START_REF] George | Three-dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF]. Further work on 3D simple points have established only connectivity of X and X is sufficient in order to characterize 3D simple points [MB92, BM94, Ber96, aBCaCaDM94, SC94]. As in 2D, deciding if a point is simple can be done in constant time in 3D. In order to do so, Bertrand and Malandain introduce topological numbers T 6 and T 26 : Proposition 3.3.2 [START_REF] Bertrand | A new characterization of three-dimensional simple points[END_REF] Let X ⊆ Z 3 and x ∈ X, let T 26 (x, X) be the number of 26connected components of (X ∩ N * 26 (x)), and let T 6 (x, X) be the number of 6-connected components of (X ∩ N * 18 (x)). In 26-connectivity, x is simple for X iff T 26 (x, X) = 1 and T 6 (x, X) = 1. In 6-connectivity, x is simple for X iff T 6 (x, X) = 1 and T 26 (x, X) = 1

Studies of simple points in 4-dimensions have also been achieved, leading once more to a local characterization of such points [START_REF] Kong | Topology-preserving deletion of 1's from 2-, 3-and 4-dimensional binary images[END_REF][START_REF] Couprie | New characterizations of simple points in 2d, 3d, and 4d discrete spaces[END_REF]. Thanks to these works, characterization of simple points up to 3D can be done once in constant time, and linear-time algorithms exist for 4D. In dimensions strictly greater than 4, it is known that local neighborhood may contain intractable configurations such as Bing's house [START_REF] Rh | Some aspects of the topology of 3-manifolds related to the poincaré conjecture[END_REF], and so simple points cannot be used [START_REF] Couprie | New characterizations of simple points, minimal non-simple sets and p-simple points in 2d, 3d and 4d discrete spaces[END_REF][START_REF] Couprie | New characterizations of simple points in 2d, 3d, and 4d discrete spaces[END_REF].

Simple sets

In general, an object possesses more than one simple points. When a simple point is removed from an object, three events can take place: non simple points can become simple, simple points can become non simple, or nothing changes. Consequently, removing two or more simple points simultaneously from an object may lead to obtaining a set of non homotopic to the original object (a process commonly called "breaking the topology"). Parallel thinning -removing simultaneously many simple points -is possible, but must be performed under certain conditions. In the following, we will give an overview of the various breakthroughs performed in the theoretical aspects of simple sets.

Given an object X ⊂ R n , we say that D ⊂ X is k-simple for X if D consists only of k-simple points of X and that there exists a discrete thinning process allowing to transform X into X\D. In 2D, this definition is equivalent to saying that X and X\D have the same connectedness (with regards to the connectivity of X) and X and X ∪ D have the same connectedness (with regards to the connectivity of X).

minimal non-deletable sets Minimal non-deletable sets were introduced by Ronse [Ron88b] in order to characterize under which conditions simple points could not be removed simultaneously from a 2D object without changing the topology.

Minimal non-deletable sets define "forbidden features" that should not appear in a set of pixels in order for it to be simple. Minimal non-deletable sets were designed in order to prove that 2D parallel thinning algorithms were topology preserving (and, therefore, valid) by testing only a small number of configurations of points. A computer-based implementation of these tests was later proposed in [START_REF] Hall | Tests for connectivity preservation for parallel reduction operators[END_REF], a 3D implementation -in [START_REF] Kong | Problem of determining whether a parallel reduction operator for n-dimensional binary images always preserves topology[END_REF][START_REF] Ma | On topology preservation in 3d thinning[END_REF][START_REF] Kong | On topology preservation in 2-d and 3-d thinning[END_REF], and a 4D implementation was proposed in [START_REF] Gau | Minimal non-simple sets in 4d binary images[END_REF][START_REF] Kong | Minimal non-simple sets in 4-dimensional binary images with (8,80)adjacency[END_REF].

p-simple points In 1995, Bertrand introduced the P-simple points in order to characterize in 3D simple points that could be removed simultaneously [START_REF] Bertrand | On p-simple points[END_REF]. To do so, given X ⊂ Z 3 and x ∈ X, he sets the geodesic n-neighbourhood of order k of x inside X (n being equal to 6 or 26) as the set

Γ k n (x, X) = ∪{Γ n (y) ∩ Γ * 26 (x) ∩ X|y ∈ Γ k-1 n (x, X)}, with Γ 1 n (x, X) = Γ * n (x) ∩ X.
Definition [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF] Let X ∈ Z 3 , P ⊂ X, x ∈ P and n equals to 26 or 6. The point x is P n -simple if, for all S ⊂ (P\{x}), x is n-simple for X\S.

Let S n (P) be the set of all P n -simple points. A set D is P n -simple if D ⊂ S n (P).

Based on this definition, one can see that, given the definition of P-simple points, if a set D is P-simple, then X and X\D are homotopic. The P-simple points allow to define sets of points that can be removed at once from an object during homotopic thinning. Let G 6 (x, X) = Γ 2 6 (x, X) and G 26 (x, X) = Γ 1 26 (x, X), the next proposition allows to locally characterize P-simple points: Proposition 3.3.3 [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF] Let X ⊂ Z 3 , P ⊂ X, x ∈ P, n be equal to 6 or 26 and n be equal to 32 -n.

x is P n -simple iff            The number of n-connected components of G n (x, X\P) is equal to 1, and
The number of n-connected components of G n (x, X) is equal to 1, and For all y ∈ Γ * n (x) ∩ P, Γ * n (y) ∩ G n (x, X\P) is not void, and For all y ∈ Γ * n (x) ∩ P, Γ * n (y) ∩ G n (x, X) is not void

As with the minimal non-deletable sets, the P-simple points allow to check if existing parallel 3D thinning algorithms work. Indeed, in [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF], the author gives a method for checking, based on the P-simple points framework, the topological validity of thinning algorithms. Moreover, P-simple points were widely used in order to propose new parallel 3D thinning algorithms (an example of a new algorithm is given in [START_REF] Bertrand | Sufficient conditions for 3D parallel thinning algorithms[END_REF]).

Critical kernels

The recent critical kernels, introduced in [Ber07], present a new framework for performing parallel thinning in 2D, 3D and 4D [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF][START_REF] Bertrand | Two-dimensional parallel thinning algorithms based on critical kernels[END_REF]. Critical kernels were also used to propose a new definition of simple points in 2D -4D [START_REF] Couprie | New characterizations of simple points in 2d, 3d, and 4d discrete spaces[END_REF], and links between this framework, P-simple points and minimal non-deletable sets were established in [START_REF] Couprie | New characterizations of simple points, minimal non-simple sets and p-simple points in 2d, 3d and 4d discrete spaces[END_REF]. Critical kernels were also used to prove that some thinning algorithms were valid, while others were not correct [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF].

Although they have applications in the voxel framework, the critical kernels rely on the cubical complex framework that will be presented in the remainder.

Simple pairs and non simple points

In order to conclude this overview of simple points, let us quickly talk of recent developments showing that simple points are not the only interesting elements in homotopic thinning. In 2D, it is always possible to reduce a set homotopic to a disk to a single point by removal of simple points. In 3D and up this is no longer the case, as there exists configurations that have no simple points but are still homotopic to a sphere (or hypersphere in dimensions greater than 3). This means that skeletonization procedures that proceed by removing simple points can get stuck on intractable configuration like Bing's house. It is therefore necessary to move beyond the notion of simple points.

Recent work from Passat et al. ([PCB08]

) is based on critical kernels in order to exhibit a new configuration called simple pairs. In such pair of points, none of the point is simple but the pair itself is simple and can be removed without changing the topology of the object.

Other works have shown that, in 3D, some points are not simple (relying on the local characterization given previously) but can still be removed without changing the topology of the input [START_REF] Bertrand | A note on 3-d simple points and simpleequivalence[END_REF][START_REF] George | Three-dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF].

Thinning process

Homotopic thinning in the digital framework consists of removing simple points from an object, until either no more simple point can be found, or a satisfactory subset of voxels has been reached : this method will reduce a ring into a circle, or a ball into a single point. An illustration of this process is displayed in Figure 19.

Two main strategies are possible for removing simple points: sequential removal and parallel removal.

Sequential algorithms

Sequential removal of simple points can be achieved by detecting simple points in an object, and removing them one after the other, until no more simple point can be found. After removing a simple point, the new set of simple points of the object must be computed. Such basic strategy does not guarantee the result to be centred in the original object, and does not preserve the "visual characteristics" of the object during thinning. It is important, when designing a sequential thinning algorithm, to decide of a removal order of simple points, and of a strategy for preserving interesting visual features of the object.

In order to obtain a centered skeleton, one must define a precise order of removal for simple points. Usually, in order to get a centered skeleton, it is necessary to delete simple points "layer by layer", from the outer layer to the inner one. Many strategies have been proposed in 2D for deciding of a removal order (a very exhaustive survey of thinning methods in 2D before 1992 can be found in [START_REF] Lam | Thinning methodologies-a comprehensive survey[END_REF]): for example, in the 80s, it was proposed to follow an object's contour in order to find and remove simple points "layer by layer" [START_REF] Arcelli | Pattern thinning by contour tracing[END_REF][START_REF] Pavlidis | A thinning algorithm for discrete binary images[END_REF]. However, this thinning scheme (as many more) can hardly be generalized to 3D and up.

A widely used strategy to obtain a centred skeleton with a thinning process consists of computing a priority function on the object and removing the simple points of X according to the value of this function [START_REF] Davies | Thinning algorithms: a critique and a new methodology[END_REF]: at each step, the simple point that is removed is the one with the lowest value. The euclidean distance map is widely used as priority function in order to remove points "layer by layer" [Tal92, CCZ07, MFV98].

Distance transform (DT), first introduced by Blum [Blu67], has a vast application in skeletonization. The DT-based thinning method was proposed by Rosenfeld et al. [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]. Other works use discrete distances, such as the chamfer distance ([Pud98]), to decide a removal order. An example of a distance transform is illustrated in figure 20.

One can use a thinning scheme based on a priority function, called inhibitor set, which is a set of points of the input object which should be in the resulting skeleton. An inhibitor set allows to choose "anchor points" for the skeleton, and therefore preserves the visual aspect of the original object in the skeleton.

Algorithm 1 shows the basic thinning scheme based on a priority function. The set W, called inhibitor set, is a set of points of the input object which should be in the resulting skeleton, and D is the priority function used to decide an order of points removal (here, the lower priority means faster removal, so it is possible to use an euclidean distance map as a priority function).

However, when performing a thinning guided by an euclidean distance map, the points of the inhibitor set and the directions of thinning followed by the algorithm are not always "compatible". In [START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF], the authors use a thinning algorithm where the slope of the priority function is used to dynamically add points to the constraint set. In [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF], the authors propose to merge the slope calculation into the priority function, leading to a new priority function and a new thinning algorithm which works in 2D and 3D.

Parallel algorithms

As previously explained in section 3.3.1.1, removing simple points simultaneously from an object usually "breaks" the topology. However, as shown previously, various theories have been elaborated in order to characterize sets of simple points which can be removed at the same time. It should be noticed that when performing parallel thinning, priority functions are rarely used as points are naturally removed "layer by layer".

Rutovitz was the first to propose a (fully) parallel thinning algorithm in 1966 [START_REF] Rutovitz | Pattern recognition[END_REF]. However, it is well known that Rutovitz's algorithm does not always preserve topology, but there exist "patches" for correcting it ( [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF]).

It is in 1981 that Pavlidis published the first fully parallel thinning algorithm ([Pav81]), in 2D for 8-connected objects, that was later proved to preserve topology in [START_REF] Couprie | Note on fifteen 2d parallel thinning algorithms[END_REF]. In this algorithm, the author defines multiple pixels, corner pixels and contour pixels.

Given X ⊂ Z 2 and x ∈ X, the point x is a contour point of X if it has an element of X in its 4-neighbourhood.

According to Palágyi [START_REF] Palágyi | A 3d fully parallel surface-thinning algorithm[END_REF], parallel thinning algorithms can be divided into three categories. In the first group, directional thinning algorithms, the main loop is divided into sub-iterations, and the thinning operator (the considered configurations of simple points) is changed from one sub-iteration to another. In subfield-based thinning algorithms, the points of the object are decomposed into subsets, and at at a given iteration of the algorithm, only simple points in a given subset are studied. Finally, in fully parallel thinning algorithms, no sub-iteration takes place: the same thinning operations (which usually remove sets of simple points) are performed on the object at each iteration of the main loop.

Preserving the visual aspect

In order to preserve the visual aspect of the initial shape in the skeleton and in order to compensate for generation of unnecessary branches of thinning and distance transform-based methods, that are sensitive to noise, three strategies can be adopted.

The first strategy consists of choosing, before the thinning process starts, an inhibitor set whose points won't be removed during skeletonisation. The set points will therefore act as anchor points for the skeleton, and if these points are well chosen, the final skeleton will have a satisfactory visual aspect. There are many different ways of computing an inhibitor set, such as using a filtered medial axis [START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF] or the discrete bisector function [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF].

The second strategy is to detect, during the thinning process, interesting points which should be included in the skeleton. Once such points have been detected, they are added to an inhibitor set. Although this strategy is mostly used with parallel thinning algorithms, it has been also used in works using sequential thinning ([Tal92], [START_REF] Couprie | Discrete bisector function and euclidean skeleton in 2d and 3d[END_REF]).

The last strategy consists of performing skeleton pruning once the thinning algorithm is finished. This strategy is used after performing a thinning where some points were retained during the skeletonisation (using one of the two previous strategies explained before). Skeleton pruning takes place after a skeleton with too many points is obtained, and consists of removing some points of the skeleton without changing its topology. In [START_REF] Vincent | Efficient computation of various types of skeletons[END_REF], the authors give multiple criteria for filtering a 2D skeleton, such as the size of its elements (the curves composing the skeleton), or the quantity of information from the original shape contained in each element (called area by the authors).

From skeleton to graph

The ability to distinguish the different components of the curve-skeleton depends on the ability to detect the junction points, i.e., the points where two or more curves meet. From this decomposition, one can infer the corresponding part structure of the original object.

Some thinning algorithms directly classify the skeleton points as junctions, either during thinning or as a post-processing step [ [START_REF] Bertrand | A new topological classification of points in 3d images[END_REF]].

From the geometric algorithm class, level-set methods directly identify the joints as the centroids of level-sets. Joint locations depend on the function used to define the level sets.

Here, a segment is a piece of vasculature connecting two vessel junctions, where vessel endpoints were considered as junctions. This mapping of segments and junctions can be represented by a mathematical graph.

A graph in this context is defined as a tuple (V, E), where V is a set of vertices and E a set of pairs of vertices forming edges. In a directed graph, the two vertices of an edge are ordered, the first vertex called the source, the second being the target. The degree of a vertex is defined as the number of edges associated with that vertex. A vascular network may now be represented as segment graph by associating each junction and each vessel dead end with a vertex and each segment with an edge in the graph. In Fig. 21 an example of such a graph is illustrated.

vascular tree construction

A substantial body of work has been done in measuring and enumerating the topology of vascular trees. The progress in precision and efficiency can be especially admired in the coronary vessel analysis. Due to the amount and the form of data in these trees, techniques for quantifying structural anatomy and comprehension of physiological processes are essential.

Among some of the methods that performed vascular tree reconstruction, we can mention the work of Bullitt [BAL + 99] where the vascular trees are reconstructed in the graph-based manner. The disadvantage of the method is the fact that each branch should be initialised by user. If the user does not inspect the image data carefully, a faint but important vessel may be missed. A more automated approach is preferable.

For tubular structures such as vessels, Krissian et al. [KMA + 99] applied multi-scale analysis to raw intensity image to detect vessel structures of various sizes. They proposed a "medialness" measure at a given point and scale based on the eigenvalues and eigenvectors of the Hessian matrix of the image. An adaptive medialness measure was proposed to be the mean of first order derivative information at the fitted circle.

[AJ01] has achieved interesting results with an interactive method for obtaining 3D centerlines given 2D user inputs. The problem is solved as a constrained optimization problem using a discrete graph with a shortest-path algorithm. But the price of interactivity with 3D data is expensive.

An alternative strategy to track the branches would be a method, as in [START_REF] Carrillo | Extraction of 3d vascular tree skeletons based on the analysis of connected components evolution[END_REF], where vascular tree skeletons are detected the bifurcations by analyzing the binary connected components on the surface of a sphere that moves along the vessels. The method assumes that the vessels can locally be separated from the background by an appropriate adaptive threshold.

Some approaches surpass the need to deal with segmentations and skeletonizations for medial-line extraction, as in [YMS + 00, BB08] for curve skeleton extraction directly from grey value images that. The approach of Bauer [START_REF] Bauer | Extracting curve skeletons from gray value images for virtual endoscopy[END_REF] is single-scale and its efficiency is boosted by a Gradient Vector Flow (GVF) GPU implementation.

Vascular tree analysis

There are three main approaches for vascular network representation as a graph. In the first, the minimal paths are searched between the root and the ends of the vessel branches. This is often done with Dijkstra's algorithm within skeletonization schemes, such as topological thinning or distance transforms. The limitation of such methods is discrete-voxel accuracy. In the second, top-down, approach, vessels are subsequently tracked from the largest ones to the smallest, as in [FDCR01, NBB + 06].

And the third approach is to independently track detected vessel segments and then organize them in a network, as in [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF]. In comparison with the top-down approach, this requires additional work on vessels identification and their segments reconnection.

Hierarchy (the ability to create a family of curve-skeletons of increased complexity) is not achievable using thinning algorithms because when processing a voxel there are only two choices: keep it, or remove it. A skeleton is obtained only after the last iteration of the algorithm. For this purpose, a number of ordering methods have been section suggests, imparts error associated with commuting subtrees with the entire vascular tree. A method was designed to assess how these errors affect results of the Strahler classification. For example, Fig. 3 illustrates common artifacts seen in medical images of vasculature (denoted with ‡). The gray marker infers a bifurcation where incorrect ordering has no effect on the upstream ordering, as changing the first-order vessel to a terminal zero-order vessel results in the same upstream ordering. Consequently, the probability that vessels were ordered incorrectly due to this terminal is zero. The dark markers, however, infer bifurcations where not applying the iterative technique entirely changes the ordering upstream as well as the calculated statistics (i.e., changing these terminal vessels to zero order causes a restructuring of the upstream ordering).

To assess these errors one might change terminal orders of the vascular network and compute the effect. This process, however, would be computationally intensive, especially for large-scale vascu-lar trees. The framework developed provides a means for quantifying the sensitivity of a vascular tree to likely perturbations of terminal vessels, assuming the given ordering scheme is a close match to the correct ordering.

We define i ϭ [ I i S i D i ] T to be a vector of probable outcomes for a vessel segment Si, where the components are the probabilities of that segment increasing ( I i ), remaining ( S i ), or decreasing ( D i ) an order. The transference of probability through vessel segments (and consequently the overall effect) is mediated by the series of bifurcations, or operators, and Strahler ordering rules at each. These operators summarize all possible outcomes given for the daughter segments and modify the upstream segment probabilities accordingly. Operators, denoted as , have subscripts that label a specific branch point and the superscript T referring to the matrix transpose (for complete details on operators, refer to APPENDIX B).

When calculating the probability profile of the jth terminal vessel segment Sj, with radius rj that is of order i, we define the upper (rmax) bound by the same diameter defined criterion proposed by Jiang et al. 6 (19) shown in Eq. 1, where ri and i are the average radius and standard deviation of order i vessels.

r max ϭ ͑r i ϩ r iϪ1,iϩ1 ϩ iϪ1,i Ϫ i,iϩ1 ͒/2 (1)
Note that the lower bound for order iϩ1 is set as the upper bound for order i. Using a normal distribution to describe error in ordering terminal vessels (with mean rj and standard deviation ␣rrj, where ␣r is a constant), we have assumed that terminal probabilities can be assigned by the area under the curve delineated by the bounds in Eq. 1 as seen in Fig. 4.

To calculate i or the probability distribution of the ith network segment (Si), the method requires the definition of a pathway ti. Considering the subtree beneath Si, ti is an index of segments that defines a single path down the ith subtree to a terminal vessel. For example, if one such path from Si to a terminal was given to go through segments S1, S2, and S3, ti ϭ [S1 S2 S3]. Choice for defining ti is arbitrary but requires only movement downstream. A different path for the same segment (Si) results in a modification of operators (which are defined by a chosen branch or branch pair) but no 6 Jiang et al. ( 19) applied Eq. 1 to reduce the radii overlap between successive orders and modify the process of assigning Strahler orders. It is important to note that Eq. 1 is applied, in this study, in a different context. proposed to analyse vascular data. The most common one is the Strahler Ordering method proposed by Strahler [START_REF] Strahler | Quantitative analysis of watershed geomorphology[END_REF] and Horsfield [START_REF] Horsfield | Morphometry of the small pulmonary arteries in man[END_REF].

Strahler ordering begins by labelling all terminal vessels (those that do not have children). Following vessel segments upstream, a second-order segment is formed by the joining of two first-order segments, a third-order segment by the joining of two second-order segments, and so on. Illustration can be seen in Figure 23.

This scheme is totally based on the topology of the vascular structure. Strahler ordering provides insight into vessel distribution of connectivity. Moreover, it can help in understanding of certain vessel disease and visualization of the vascular data.

visualization

In clinical applications and for medial education, it is more essential to facilitate viewing of anatomic and pathologic structures. Another need is to improve the 2D data visualizaion in 3D information, so that there be a support connected and reliable information between the two.

Four volume visualization techniques have proven useful over the years for clinical angiographic practice: Multi-Planar Reconstruction (MPR), Maximum Intensity Projection (MIP), Surface Rendering (SR) and Volume Rendering (VR).

In Multi-Planar Reconstruction, from a 3D volume, 2D slices are shown on the screen. These slices are the intersections of the volume by planes, either orthogonal to the original acquisition frame, or oriented arbitrarily. In the latter case, interpolation of the original data is needed to produce the greyscale values of the intersected pixels. As MPR displays the original, possibly interpolated, grey-level values, it is always the visualization method that is the reference in case of doubt. It is also possible to display the intersection of the data volume with any surface: the technique is called Curved Planar Reformation (CPR). It is particularly useful in the context of angiography, where the intersecting surface can be chosen to include the centerline of the vessel of interest, thus showing a cut of the vessel lumen along the vessel course. Pushed even further, the technique can produce visual renderings of the lumen of a whole arterial tree in one image [START_REF] Kanitsar | Advanced curved planar reformation: Flattening of vascular structures[END_REF].

Maximum Intensity Projection is the most common visualization technique aside from MPR. It shows, in a given direction, the maximum gray-scale value along a ray that is shot in that direction intercepting the volume. When the direction does not coincide with the original volume frame, interpolation is needed to produce grayscale values along the ray. MIP has the following advantages and drawbacks:

1. It looks like a conventional, X-Ray angiogram, therefore, it allows experienced radiologists to interpret MIP images without additional training.

2. The distribution of graylevels in the MIP is close to the one observed from the original volume.

3. The objects of interest can be overlapped by others.

It overestimates the vessel diameter [SPG96].

In surface rendering, first the objects of interest must be identified in the slice images by segmentation, which is usually using a threshold, resulting in a binary volume. The boundary of the binary volume is then first-order approximated by a set of polygons, typically using the "marching cubes" method or one of its more topologically-correct variants [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF][START_REF] Daragon | Marching chains algorithm for Alexandroff-Khalimsky spaces[END_REF] . For the perception of the three-dimensional structure it is helpful if the viewpoint is dynamically moved.

Volume rendering does not require any preliminary segmentation step. Each voxel is mapped via a transfer function to optical quantities such as color and opacity. Then a ray is sent from each pixel in the image plane through the volume data. As it passes the volume, each voxel it encounters may contribute to the pixel value, depending on its intensity and depending on the opacity of all voxels closer to the image plane. In the end, a pixel's intensity integrates contributions from many voxels along the ray. Volume rendering is an extension of the MIP idea, which can take into account edge information by including a gradient weight in the transfer function for instance.

From an other side, in 3D it is not trivial to recognize the objects. They might seem smaller, or might be occluded by other objects. In medical images, the last notion is highly prevalent, as most of the structures are surrounded by others, therefore, their view can be occluded. Change (rotation, translation) of the point views and change of the opacity of the objects' surfaces can facilitate the task. However, the first strategy might not satisfy the user, moreover it may distract him. And for the second parameter (transparency), it is often not sufficient due to possible multiple occluding objects and the very nature of transparency representing less of information, therefore being less recognizable.

When it comes to the visualization of vascular structures, due to their narrowness, limited spatial resolution of medical image data, these structures are represented with only a few pixels per slice. Therefore, a straightforward surface mapping of vessels in 3D leads to noisy visualizations. This makes the study of their topology very complicated from such visual presentations. This is an important step in surgery or treatment planning. One such example is the understanding of the risk of one of the branches being cut or connected to an other. To assist in such tasks, more abstract techniques are desirable. As such, skeletons are currently a useful tool [BAL + 99, PP03].

One of the early benefits of curve-skeletons was for virtual navigation, viewing its centredness to generate collision-free paths through a scene or through an object. With a scene composed of 3D objects, the skeleton of the background gives a free path through the scene.

The centered, thin, homotopic to the original image skeleton, which is fitted inside the vessel, can be used to assist visualization and study of vascular trees.

The example of such vessel visualization can be seen in Figure 24. According to [START_REF] Preim | Smart 3d visualizations in clinical applications[END_REF], this method for has been developed and refined in fruitful discussions with radiologists and is a good comprise between precise and easy-to-interpret visualizations.

This method (as many other previously discussed in Chapter 2) is based on the assumption that the cross-section of vascular structures is circular. This assumption, of course, is a simplification. It is not appropriate when vascular diseases should be diagnosed. As it is illustrated, vessels skeletons can be fitted more or less close to the original data.

discussion

In this chapter, we presented the skeletonisation methods together with decomposition schemes, associated computational methods and other technical details.

We also showed how the vascular tree construction can help physicians to better view and understand the vascular structure via the skeletons. And we presented the nowadays clinical visualization methods.

Thinning algorithms can either directly produce a curve-skeleton (by using curvethinning templates) or further thin a surface skeleton to a 1D representation. Parallel thinning algorithms, which remove all simple points at once, may not be able to achieve 1D skeletons due to homotopy constraints. An illustrative 2D example is the case of a rectangle whose width is an even number of voxels: in the last step of the thinning process, the middle section of the skeleton will be a curve of width 2. Although all points of this curve are simple points, removing them would completely remove the middle section. At this stage, no other simple points can be removed and the skeleton is not 1D. However, when it comes to precision of skeleton methods for clinical applications, thinning and general field methods do not guarantee centering of the result. In the case of directional thinning, this would depend on the order in which the different directions are applied. In the case of general field methods, since they are taking into account a larger surface area than the two closest points, centering is usually compromised. This means that for the measurements, this has to be taken into account.

Methods using a distance field can better achieve centering because of the information included in the distance field.

Moreover, resolution affects any centering measurement in the discrete domain. Using the same example of a shape whose width is an even number of voxels, if this shape is reduced to a 1D skeleton, at the grid resolution, the curve-skeleton must be one voxel closer to one of the sides than to the other.

Another aspect of proving the usefulness of the vascular tree representation is in the comprehension between the medical and technical specialists. As this demands a great effort in profound studies of the unknown fields for both sides, and due to the importance of each result, the research in this field takes longer to be proved. This in itself is a motivation for further progress in this field. -Paul Klee

Part II M O R P H O -H E S S I A N F I LT E R A N D I T S A P P L I C AT I O N S

In this chapter, we present our morpho-Hessian filter used for the enhancement of curvilinear structures applied prior their segmentation, analysis and visualization.

Our primary goal is in detection of thin elongated objects and their directions. For this, we have chosen to use second order directional Gaussian derivatives in their Hessian form. Such methods provide us not only with the detection of thin objects and their principal directions, but also with the possibility to do it on different scales. Further, we are making use of several tubularity measures by Sato et al. [SNS + 98] and Frangi et al. [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF].

Then, we are interested in enhancement and reconnection of tubular objects. The spatially variant mathematical morphology (SVMM) can be helpful in this case. The theory of SVMM and its corresponding algorithms are formulated with the purpose of filtering an image differently at various positions. In the case of elongated objects, we wish to discover their local orientation and filter them locally along this direction, for instance, with an oriented segment. We then expose a simple version of the theory of SVMM together with an efficient algorithm of its operations.

In the following chapter, we combine these two filters and perform their performance evaluation first on a 3D simulated data and on clinical data.

feature detection with gaussian derivatives

In order to detect one object or a pattern, one might start by detecting some features of the image in order to locate the object of search. For this purpose, local operators called feature detectors can be used. They are local operators that are applied to an image in order to detect and label local responses as particularly informative. The design of feature detectors is a fundamental problem in image analysis. The possible gain from a good feature detector is to guide the visual system to a few positions in an image that are not only labelled as particularly informative but that also are particularly useful in the image interpretation.

One of the pioneering and famous techniques is by Canny [START_REF] Canny | Finding edges and lines in images[END_REF], where he seeks an optimal balance between localization of interesting features and noise responses.

One of the bigger challenges in image analysis is to construct operators that are translation-, rotation-and scale-invariant. Translation invariance can be satisfied by any morpho-hessian filter convolution kernel. Rotation invariance can be guaranteed either by rotation-invariant kernels or when the preferred direction is fixed relative to the image. Scale invariance can be satisfied by derivatives of Gaussian filter. Linear combinations of derivative of Gaussian filter kernels are the basic feature detectors within linear scale-space theory.

From the first to higher dimensions, derivatives may be computed along different directions. To calculate these differential operators of the image I(x) in a well-posed fashion, we use here service of linear scale-space theory [START_REF] Koenderink | The structure of images[END_REF]. According to this theory, differentiation is defined as a convolution with derivatives of Gaussians:

∂I(x, σ) ∂x = σ γ I(x) * ∂G(x, σ) ∂x , (4.1)
where the N-dimensional Gaussian is defined as:

G(x, σ) = 1 (2πσ 2 ) N\2 e - |x| 2 2σ 2 (4.2)
The second order derivative of Gaussian kernel at scale σ probes the contrast of the regions inside and outside the kernel of size σ in one direction (as in Figure 25).

The Taylor expansion of a function in a neighbourhood of x o is where ∇ is the gradient vector (F x , F y , F z ) T in 3D and H is the Hessian matrix:

F(x o + δ x o ,σ ) ≈ F(x o , σ) + δx T o ∇ x o ,σ + δx T o H x o ,σ δx o , (4.3
H =     F xx F xy F xz F yx F yy F yz F zx F zy F zz     . (4.4)
In 2D, the gradient vector is (F x , F y ) T and the Hessian matrix is:

H = F xx F xy F xy F yy , (4.5)
Eigenvalue analysis can be performed on the Hessian matrix in order to extract one or more principal directions of the local structure of the image. This gives the direction of the minimal curvature, the principal one in the tubular structure and a high curvature in the vessel cross-section plane.

Compared with the image gradient, whose response is independent of the shape and local structures of boundaries, the Hessian matrix can capture the shape characteristics of objects, such as tubes, planes, blob surfaces or noise.

By the principal value analysis, the Hessian matrix H can be decomposed into three principal values, λ 1 , λ 2 and λ 3 (

|λ 1 | |λ 2 | |λ 3 |
) and e 1 , e 2 and e 3 are their corresponding eigenvectors, respectively. When λ 1 is close to zero and much smaller than λ 2 and λ 3 , the principal direction of the (bright) vessel is given by e 1 , the eigenvector corresponding to λ 1 . Vectors e 2 and e 3 form a base for the orthogonal plane.

Here, we are interested in detection of thin elongated objects and their directions. We are considering bright objects on a darker background. Accordingly, we will ignore all objects whose λ 2 and λ 3 are negative at the same time.

To summarize, in the Table 2 the relations of Hessian eigenvalues are shown for different shapes and intensities for 2D and 3D images.

morpho-hessian filter

All three eigenvalues are important in detecting tubular objects. Combinations of those lead to expressions differentiating curvilinear from other profiles. In particular, if appropriately designed and applied at multiple scales, such combinations, often called vesselness function, should give the strongest response at one particular scale corresponding to the plate-, blob-like and/or tubular objects [SNS + 98, FNVV98, KMA + 00].

Here, we are considering vesselness functions by Sato [SNS + 98] and Frangi et al. [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF].

Frangi's vesselness function

For a 3D image I(x), x ∈ R 3 observed at a scale σ, Frangi's vesselness function is formulated as follows:

ν(x, σ) =    0 if λ 2 > 0 or λ 3 > 0, (1 -e -R 2 A 2α 2 ) • e - R 2 B 2β 2 • (1 -e -S 2 2c 2 ) otherwise, (4.6)
with

R A = |λ 2 | |λ 3 | , R B = |λ 1 | √ |λ 2 λ 3 | , S = H σ = j λ 2 j , (4.7)
in which R A differentiates between plate-and line-like objects, R B differentiates bloblike ones, and S accounts for the intensity difference between objects and background. Parameters α, β and c influence the sensitivity of the filter to the corresponding measures.

For 2D images, Frangi's vesselness can be expressed as:

ν(x, σ) =    0 if λ 2 > 0, e - R 2 B 2β 2 • (1 -e -S 2 2c 2 ) otherwise, (4.8)
and

R B = |λ 1 | |λ 2 | , S = H σ = j λ 2 j , (4.9)
The above expressions of vesselness are probability functions of line-likeliness of each pixel. The parameter c depends on the intensity profile of the image and a quarter of the maximum intensity of the tubular objects is recommended in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF].

Often, in order to obtain the final segmentation, the vesselness can be thresholded. Such operation alone, however, does not provide satisfactory segmentation results. As stated in the Equations 4.6,4.8, the filter can be applied at different scales, which can provide a result in a larger range of line sizes. The vesselness function is normalized by σ 2 [Lin96] and the maximal vesselness is selected for each pixel x as its best response: ν max (x) = max(ν(x, σ)) with σ min σ σ max σ max (x) is, then, the scale at which the vesselness has the best response for a pixel, which can be an estimate of an object's width. However, the filter performance depends on how well the scale is chosen, which will be discussed in the following.

In Figure 26 an illustration of the vesselness performance for a tubular disconnected object and second order detected directions at each pixel of the object can be observed. It is obvious, that the elongated object is detected together with its directions, but it is not enhanced.

Sato's vesselness

The filter defined by Sato et al. [SNS + 98] was used as follows:

S(x) =      σ 2 |λ 3 | λ 2 λ 3 ξ 1 + λ 1 |λ 2 | τ , λ 3 < λ 2 < λ 1 < 0 σ 2 |λ 3 | λ 2 λ 3 ξ 1 -ρ λ 1 |λ 2 | τ , λ 3 < λ 2 < 0 < λ 1 < λ 2 ρ (4.10)
where ξ 0 influences cross-section asymmetry, τ 0 controls the sensitivity to bloblike structures, 0 < ρ 1 controls sensitivity to the tubular object curvature, and σ 2 normalizes responses across scales.

Hysteresis thresholding

As it was proposed by Canny, an alternative method to a simple threshold, would be an adaptive thresholding, or more precisely, a hysteresis thresholding. The main idea is that some edges might not respond with the same strength as their neighbours due to noise. Therefore, setting a second, lower threshold and making sure it is connected to the stronger ones, might help to get a more complete segmentation.

More precisely, at any point x being part of the object and is above the threshold T max , such a point is immediately a part of the segmented image. Moreover, all of its connected component points are also in the output if they they are above the lower threshold T min . In this case, the probability of a false positive elongated object is reduced as the high threshold can be rather big with no risk of under-segmentation.

Scale-space

As previously stated, the Gaussian derivatives are applied at some scale. Therefore, the choice and number of scales is especially important for such derivative-based methods. Too small scales emphasize small details, while producing spurious responses to noise, tangent objects and other features. Too large scales deform original objects at detection and lose details. On the other hand, it is computationally expensive to process along all possible scales.

The problem of scale-selection cannot be solved by an invariance requirement. Invariance "conserves" information under some specified transformations of the data. Scale selection ignores information: in the sense that the particularly informative scales, or positions and scales, of some operator response do not contain the same information as the original data. Of course, the purpose of scale-selection is not in the first place to ignore information, but rather to distinguish between relevant and irrelevant. This, however, amounts to ignoring the irrelevant information and, therefore, scale-selection cannot be derived from an invariance requirement.

The first work that dealt with positions and scales at the same time was proposed by Lindeberg [START_REF] Lindeberg | On Scale Selection for Differential Operators[END_REF]. And he has suggested such a heuristic principle:

In the absence of other evidence, a scale level at which some (possibly nonlinear) combination of normalized derivatives assumes a local maximum can be treated as reflecting the characteristic length of a corresponding structure in the data.

This idea has inspired works on multi-scale filtering with early popular approaches used for vessel enhancement [KGSD95, LCB + 97, SAH + 98].

In 1998, Lindeberg proposed to use a parameter γ-normalized derivative operator:

δG(x, σ) δx γ-norm = σ γ 2 δG(x, σ) δx (4.11)
This operator ensures that the filter responses are comparable between different scales.

In [START_REF] Majer | On the influence of scale selection on feature detection for the case of linelike structures[END_REF] they have used the parameter γ by Lindeberg [START_REF] Lindeberg | On Scale Selection for Differential Operators[END_REF][START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF] to give a priority to ridges and not edges. Figure 27 illlustrates that in the variable scale setting the response to edges of a second order Gaussian derivative operator can be turned off simply by a suitable choice of γ parameter. To better illustrate this point the shortest ridges of the grass are not displayed in image 27c. 27b displays all the detected ridges. Moreover, although these ridges cross each other in the projection of Fig. 27, they in fact occur at different scales.

spatially-variant mathematical morphology operators

Spatially-variant (SV) basic mathematical morphology operators, in this chapter, are the four standard basic morphology operators: erosion, dilation, opening and closing, using structuring elements (SE) that are not invariant by translation, i.e. that differ according to the location of their origin. Such operators are useful because they can reflect local content in actual images, which is seldom stationary except in some statistical models. For instance, structuring elements may vary according to orientation or perspective. Many well-known operators, such as some algebraic openings and closings (e.g. area operators) are spatially-variant in nature, but do not derive from spatially-variant erosions and dilations, rather from compositions of spatially invariant openings and closings.

Spatially-variant mathematical morphology (SVMM) has been known for a long time (it is mentioned in Serra's 1982 book), and even recently has been the topic of a few publications [CCS94, BCCS08, BS08], and was significantly mentioned in a special session at ICIP 2009 on adaptive morphology [START_REF] Maragos | OVERVIEW OF ADAPTIVE MORPHOLOGY: TRENDS AND PERSPECTIVES[END_REF].

Many researchers and practitioners wish to be able to use adaptive structuring elements that vary according to the location in the image. Applications include adaptive filtering [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF], segmentation taking advantage of perspective information [START_REF] Dokládal | Grey-scale morphology with spatially-variant rectangles in linear time[END_REF] or local orientation [START_REF] Tankyevych | Spatially variant morpho-Hessian filter: efficient implementation and application[END_REF][START_REF] Verdú-Monedero | Spatially-variant directional mathematical morphology operators based on a diffused average squared gradient field[END_REF].

While using spatially-variant erosions and dilation may be relatively easy, the same cannot be said of openings and closings using compositions of spatially variant erosions and dilations. Indeed, the computation of adjunct operators for spatially-variant erosions and dilations is not trivial.

There are still misconceptions about the concepts. Most publications do not mention adjunctions at all (e.g. [START_REF] Bouaynaya | Theoretical foundations of spatially-variant mathematical morphology part ii: Gray-level images[END_REF], by no means an isolated example). Some state that using spatially variant but symmetric structuring elements is sufficient to ensure openings and closings with adjunction properties, others suggest that using some SV SE families can lead to adjunct operators with infinite extension. Very few publications mention algorithms and even fewer efficient ones. The purpose of this section is to introduce the notion of adjunction in the SVMM case, to describe and justify an efficient algorithm first proposed by [START_REF] Lerallut | Image filtering using morphological amoebas[END_REF], but not justified enough before. We wish to justify the approach of [LDM07] thoroughly, and argue that if a SV family of SE leads to computable erosions and dilations, then their adjunct is also computable. Finally, we propose a simple and effective algorithm to compute the adjuncts.

An illustrative example

Basic filters of the mathematical morphology framework are not suited for oriented thin objects, as they do not take into account their directions. An example can be seen in 28, where a disconnected vessel is illustrated. One of the ways to try to reconnect the object could be using the closing morphological operation with a 3 × 1 vertical line segment. However, being translation-invariant, this does not give a desirable result as the operation does not follow the principal object's direction. In contrast, a 3 × 1 spatially variant segment utilizing the local direction information is much more effective. Here, we are interested in propagating and filtering the object along its main direction.

Let us illustrate an other simple example on the fact that SVMM can lead to filters that are not equivalent to their translation-invariant counterparts.

Let us consider the case of the colour reconstruction in modern digital cameras. Most use a so-called colour filter array (CFA) in front of a single luminance sensor, with the result that only one of the colour components (either red, green or blue) is captured in any single location. Most use a variant of the so-called Bayer array, illustrated in Fig. 29. We notice that on this CFA, green is twice more common than red or blue. This is meant to reflect the physiological fact that most humans are more sensitive to green than either red or blue.

If we want to reconstruct a complete colour image, we need to interpolate the missing components at every pixel location. If we do this in the classical way, using Let (i, j) be the local 2D location, and G(i, j) the value in the green channel pixel at location (i, j), then if G(i, j) needs to be interpolated, and |G(i -1, j) -G(i + 1, j)| > |G(i, j -1) -G(i, j + 1)|, we say that the local gradient is horizontal, and vertical if not.

Then we define the following SVMM filter:

1. If G(i, j
) is known (i.e. does not need to be interpolated), we perform an identity operation (i.e. a closing by a 1 × 1 structuring element).

If G(i, j

) is unknown, and if the local gradient is vertical, we perform a horizontal closing with a 3 × 1 structuring element.

3.

If the local gradient is horizontal, we perform a vertical closing with a 1 × 3 structuring element. This operator behaves adaptively and interpolated the green pixel values in the direction of the smallest gradient. This idea can be extended to the other two channels as well. The results can be seen on an excerpt of the light-house picture (Fig. 30) in Fig 31.

We see here that the results are much better than with the aforementioned methods. This SVMM result could not easily be reproduced by standard translation-invariant operators. Here the key idea is to use the local orientation to define a local structuring element to use: either a single point if interpolation is not needed, or a vertical or horizontal short segment to interpolate in the direction where intensities vary least, i.e. where the local texture is the most constant.

The question now remains how to actually implement this spatially-variant closing.

Morphological operators and adjunctions

While the basic morphological erosions and dilation are interesting by themselves, they are even more useful when combined to create openings and closings. However not every dilation and erosion can combine to form an opening or a closing, only those forming adjunctions do.

The notion of adjunction is central to morphology, as it encompasses the dual nature of the theory. Morphological filters come in pairs: those that work on the foreground and those that work on the background.

Basic concepts

A good, recent introduction to mathematical morphology can be found in [START_REF]Mathematical Morphology: from theory to applications[END_REF], especially chapters 1 and 2. We reuse some of their notations here. Since all of the exposition in this section is classical, we show no proof, but we refer to fundamental texts on this work, most notably [START_REF] Serra | Image analysis and mathematical morphology[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology -Volume II : Theoretical Advances[END_REF][START_REF] Matheron | Random Sets and Integral Geometry[END_REF][START_REF] Heijmans | Morphological image operators[END_REF].

We place ourselves in a space E equipped with a partial ordering forming a complete lattice, i.e. where every subset S of E has an infimum in E denoted S, and a supremum in E denoted S.

For illustrative purposes, we will often consider the binary case, i.e. the case where E is the Euclidean space of dimension n denoted R n , where subsets are arbitrary parts of R n , and where the supremum coincide with the union and the infimum with the intersection.

The notion of adjunction is central to modern morphology. It is defined in the following way: adjunction Let L and M be two complete lattice (equal or distinct). Two operators δ : L → M and ε : M → L form an adjunction (ε, δ) if and only if:

∀x ∈ L, ∀y ∈ M, δ(x) y ⇐⇒ x ε(y)
This concept is enough to introduce the basic operators of morphology: erosions and dilations Here the operator δ commutes with the supremum operator and is called a dilation, while ε commutes with the infimum and is called an erosion. and also their compositions: openings and closings If there exists two operators δ, ε that form an adjunction, then the operator γ = δε is called an opening, and the operator ϕ = εδ is called a closing.

From these, we can deduce the following properties : Proposition 4.2.1 (uniqueness) if (ε, δ 1 ) and (ε, δ 2 ) are adjunctions between L, M, then

δ 1 = δ 2 . Similarly if (ε 1 , δ) and (ε 2 , δ) are adjunctions, then ε 1 = ε 2 .
and the basic properties of the openings and closings: Proposition 4.2.2 (Properties of openings and closings) ∀x ∈ L, 1. idempotence : γγ = γ and ϕϕ = ϕ 2. increasingness : if x y, then γ(x) γ(y) and ϕ(x) ϕ(y)

extensivity: x ϕ(x)

4. anti-extensivity: γ(x) x.

While the above is extremely useful algebraically, it is also useful in practice because the behaviour of morphological operators is well defined. It is therefore possible to define new operators such as gradients or top-hats exploiting differences between these operators.

Translation-invariant operators using structuring elements

While it is possible to implement some algebraic morphological operators, most often simple erosions and dilations are constructed from operations with structuring elements. These are "special" sets whose shape, size and orientation are known in advance. Most definitions of structuring element-based morphological operators use the translation-invariant (TI) case. binary case The binary case is illustrative, so we introduce it first. Let B be a subset of E, i.e. an element of P(E). We denote the translate of B by p ∈ E is the set

B p = {x + p|x ∈ B}.
Here p is a translation vector. The morphological dilation of a set X ∈ P(E) by B is denote by δ B (X) or X ⊕ B and is given by:

δ B (X) = X ⊕ B = b∈B X b = x∈X B x = {x + b|x ∈ X, b ∈ B}.
The resulting dilation is the union of the B p such that p belongs to X: δ B (X) = {B p |p ∈ X}. As a consequence, the dilation of X by B 'enlarges' X, hence the name of the transform. In the formula, X and B play symmetric roles. Note also that when B is untranslated, (i.e. B o ), it is located somewhere relative to the origin of the coordinate system. We usually associate this point with B itself and call it the origin of the structuring element. When B is translated, so is its origin.

The morphological erosion of X by B is defined by:

ε B (X) = X B = b∈B X -b = {p ∈ E|B p ⊆ X}.
The erosion of X by B is the locus of the points p such that B p is entirely included in X. An erosion 'shrinks' sets, hence its name.

Note that in the binary case, the adjunction property reduces to:

∀X, Y, B ∈ P(R n ), X ⊕ B ⊆ Y ⇐⇒ X ⊆ Y B
In other words, the dilation of X is included in Y if and only if X is included in the erosion of Y.

general translation-invariant case There are multiple ways to generalise the binary case to arbitrary complete lattices L, but one simple way is to keep considering so-called "flat" structuring elements (in effect essentially binary sets), and switch the union and intersections above to the general supremum and infimum ones. This is legitimate and can be easily interpreted in geometric terms, as documented for instance in [START_REF]Mathematical Morphology: from theory to applications[END_REF], chapter 1.

We consider the case where L ≡ M is a family of functions, or images, defined as a mapping from the support D to the set of values V. Let , . ε : L → L be a dilation and an erosion of f ∈ L given by:

δ B (f) = p∈B f p (4.12)
and the erosion:

ε B (f) = p∈B f -p (4.13)
Here p is a translation vector, f b denotes the translation of f by b computed as f b (x) = f(x -b). B stands for the structuring element. In standard mathematical morphology, B is translation invariant, defined as B ⊂ D.

Notice the -p in the erosion with respect to the dilation, meaning that in the erosion case, the symmetric structuring element with respect to the origin is in fact used. We may want to make this fact more apparent with the following equivalent definition for the erosion:

ε B (f) = p∈ B f p (4.14) with B = -B = {-p|p ∈ B}.
Here B is called the symmetric of B or more generally the transpose of B. We will see that with a slightly more general definition of the transpose, the above definition apply to the spatially-variant case.

compositions With the definitions given above, the morphological erosion and dilation form an adjunction. In other words, we can now define the morphological opening and closing.

In the binary case, the opening of set X by B is defined in this way:

γ B (X) = X • B = (X B) ⊕ B = {B p |p ∈ E et B p ⊆ X} .
and the closing of X by B is defined by:

ϕ B (X) = X • B = (X ⊕ B) B. (4.15)
In the general case, the opening is simply γ B = δ B ε B and the closing is

ϕ B = ε B δ B .

Spatially-variant operators

We now consider the relatively unusual spatially variant (SV) dilation and erosion.

In general, one can construct a simple spatially-variant (SV) morphological operator by computing a max or min filter using a structuring element (SE) that is not the same everywhere in the image. In the context of SVMM, B is often denoted as structuring function and is defined as B : D → P(D), P being the collection of all subsets of D. As above, B(x) denotes the structuring element B originating at point x.

For instance, one can use parametrized disks or oriented segments, and vary respectively the diameter or the orientation according to a scalar field. Erosions and dilations with these SEs pose no problem by themselves. However, for filtering, openings and closings are the more interesting operators, but adjunct respective dilations and erosions are not trivially computed, as we shall see shortly.

In the general SV case, the transpose of a SE B, still noted B is now defined by:

B(x) = {y | x ∈ B(y)}, (4.16)
where x and y are points, and B(y) is the spatially-variant structuring element originating at point y. In the translation-invariant (TI) case, we observe that the general definition collapses to the one given above, i.e, B = -B, the symmetric of B with respect to the origin. However, for the SV case this is not true. We also note that the transpose of the transpose of B is B itself, i.e. B = B.

As before, for the spatially-invariant case, we denote (δ, ε) an adjunct pair of operators, with the same definitions as before:

δ B (f) = p∈B f p ε B (f) = p∈ B f p ,
bearing in mind that the definition for B is the one in eq. 4.16. With this, the opening γ B and the closing ϕ B are defined in the usual way from the definition of the adjunction:

ϕ B (f) = ε B (δ B (f)); γ B (f) = δ B (ε B (f)).
(4.17)

Note that there also exists a complementary, SE-based SV erosion ε B(f) = -δ B(-f) and its adjunct δ B(f), and that in general ε B and ε B differ.

Equivalence between morphological operators and oriented graphs

Here we consider the discrete case, modelled by a graph, as this is the practical case that interests us most, and for which illustration is easiest.

A general, spatially-variant morphological dilation or erosion in the discrete case is equivalent to a graph. More specifically, let (V, Γ ) be an oriented graph, with V a set of vertices and Γ a set of oriented edges. Then at each vertex v, Γ (v) defines a neighborhood relation, or equivalently, a structuring element with origin v. If we compute the maximum of all the neighbours of v, then we are computing a dilation. If we compute a minimum, then we are computing an erosion.

Conversely, given an arbitrary spatially-variant dilation δ B(x) on a discrete set, with x a vertex of the set, and family of structuring element B(x), then this B(x) defines all the points directly connected from x to all the vertices of B(x) and thus an edge set Γ B (x). The combination of the edge set x and vertex set Γ B defines an oriented graph.

Special case of the single-neighbour structuring element

In the following, we illustrate our operators in a restricted case for clarity, more specifically the discrete case as above, where the structuring element consists only of a single neighbour. In other word, the graph (V, Γ B ) is always such that every x ∈ V has a single successor (or neighbour).

It is easy to go from this specific case to the general case by observing, for example that a dilation with an arbitrary structuring element B(x) at vertex x is the supremum of all the dilations with single-neighbour structuring elements composing B(x).

Adjunct of the spatially variant operators

In Figure 32, we consider the case of the dilation and its adjunct erosion. For simplicity we consider only the case of the single-neighbour structuring element, as we just described. Equivalently, the SE consists of a single pixel. The arrows in Fig. 32 indicate the SE, with the origin at the base of the arrow and the pixel SE at the point.

In Fig. 32(a) we have the SV dilation for three sample pixel locations. At location (1) the value in the image is taken at the point of the arrow and set at the base in the resulting image, and the same thing happens at locations (2) and (3). We only notice that for the latter two pixels, the arrows point to the same pixel location (which is of course allowed!).

Following the definition for B necessary for the computation of ε B , in Fig. 32(b), for the computation of the erosion at location (1 * ) we use the value at location [START_REF] Avizo | Image J, FiJi 3. ITK Snap[END_REF]. For the computation of the erosion at location (2 * ), we need to remember that both pixels at location (2) and (3) pointed to (2 * ), and so B(2 * ) consists of the two pixels locations {(2), (3)}. This is conceptually simple, but becomes quickly more complicated as the size of B increases. Computing openings and closings with SV adjunct erosions and dilations computed in the way described above requires the transpose to be redefined everywhere, which can be computationally expensive requiring an exhaustive search or dedicated data structures. This becomes prohibitive for large, nD data (n 2). In the next section we present a way to decrease the computation complexity.

Adjunct by conditional propagation

In the following, we describe the SV dilation δ and its adjunct erosion ε both using spatially-variant SEs.

The transpose of B, denoted B, is defined in Eq. 4.16 and used in Eq. 4.17. In addition to being computationally costly, it can, under some conditions, be of arbitrary extent depending on the B family, which becomes problematic for forming filters based on adjunctions of dilations and erosions in order to make a closing or an opening.

In the following we propose an implementation of the adjunct operator identical in complexity to the normal one.

SV-closing implementation

The inf/sup-of-functions in Eqs. 4.17 and 4.17 are usually computed sequentially in raster scan order according to this scheme:

[δ B (f)](x) = max b∈B(x) f(x -b) [ε B (f)](x) = min b∈ B(x) f(x -b) = min y|x∈B(y) f(x -y)
The dilation of Eq. 4.18 is computed in O(MN), where N = Card(D) , and M = Card(B(x)). In the following erosion, from Eq. 4.18, given some x, the set {y | x ∈ B(y)} is a priori unknown and is computed in O(N 2 ) by exhaustive search. However, relaxing the sequential order of computing, the adjunct erosion of Eq. 4.18 can be computed more efficiently, as shown below.

Figure 33: Adjunct in-place erosion. We define the operator ε as follows: as we scan the image, for instance in raster order, we propagate the values at the origin of the SE (the base of the arrow) to their point (the point of the arrow) in a buffer image initialized to the original image value. We compute the min operator in place as the propagation is performed. In the text we prove ε ≡ ε.

Theorem 4.2.3 (Conditional propagation for the adjunct erosion) The spatially-variant dilation and the spatially-variant adjunct erosion can be computed at the same cost.

Proof Using Eq. 4.16 and Eq. 4.18 and permuting in the latter the roles of x and y, we obtain:

[ε B (f)](y) = min x∈B(y) f(y -x) (4.18)
This means that for computing e B (f) we only need B and not B. In practical terms, assume ε B (f), i.e. the output, originally set everywhere equal to f. One can sequentially read the input f at each point x. Considering the structuring element B(x) of origin x, for all elements y of B(x), we update the output value [ε B (f)](y) by taking the min operator between the current input value at x and the current output value at y.

It is easy to show that both ways of computing the adjunct are equivalent. As we scan the input image at x and update the value in the output image at y, we are indeed computing a min operator between all the origins of B(x) such that B(x) intersects y. So, once the whole image has been scanned, if B(x) intersects y, then x is in B(y) and vice-versa. This theorem is illustrated in Figure 33, where the adjunct structuring elements B are pointing at the selected pixel y for the efficient adjunct erosion.

This operation is performed with the minimal value updated everywhere along B as it spans the image.

This theorem leads to Algorithm 2 based on a conditional propagation of the value at the origin of each B(x). In this algorithm, the for loop (lines 4-5) computes sequentially for all x ∈ D, the dilation δ B (f) implemented by definition of Eq. 4.18. Every δ B (f) is then stored in d (previously initialized (line 3)). Notice that there is no need to store notes (1) In the illustrated algorithm the closing starts with the SE-based dilation, but it would be possible to start with the complementary, adjunct dilation δ B, followed by the erosion ε B. Note that in this case, due to the property B = B and using the definition of the SV erosion given by Eq. 4.17, the latter erosion is a "normal", SE-based erosion and not implemented through the ε definition. Note that the two operators differ in general. (2) Forming an opening starting with the adjunct erosion requires intermediary storage. (3) SE can be of any shape, not just a line segment.

In Fig. 34 the result of the spatially-variant closing can be observed on a tubular object with the directions provided by multi-scale second order derivatives. It can be stated that the object's segments were not reconnected, which can be explained by the fact that the dilation is performed on the pixels with a principal direction information. Therefore, a regularized direction field is necessary to perform reconnection. This is described in the following section.

direction field regularization

As shown in Fig. 36, in order to propagate objects in the space with the spatiallyvariant morphological closing (or, in fact, with many other filtering methods as well), it is necessary to have a dense direction vector field. In our case, the directional information is necessary only as far as the dilation can reach.

As can be noticed in Fig. 26, the directions obtained by second-order derivatives get chaotic at the end of tubular object segments, which might cause problems in methods that use them in further procedures. For this, a regularized directional vector field is highly desirable. One could use the gradient vector flow [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF] or an average square gradient within a diffusion scheme for such purposes as in [START_REF] Verdú-Monedero | Spatially-variant directional mathematical morphology operators based on a diffused average squared gradient field[END_REF]. While here, we have come up with a regularizing procedure based on several simple morphological operations.

At first, we perform the grey-scale morphological dilation on the direction field guided by the thresholded vesselness result with a spatially-variant structuring element of fixed size and line shape. By using the thresholded vesselness result, we make sure that we use the tubular objects as markers for direction field propagation. After follows the adjunct erosion -still only of the direction field guided by the segmented vesselness result -operation as described in section 4.2.2.7. This ensures an idempotent result, which guarantees that the resulting filter obeys all morphological rules. See Fig. 35 for a schematic representation of the operation. The resulting regularized field of the tubular object and its spatially-variant closing after the field regularization can be seen in Fig. 36.

algorithm

In the end, the algorithm of the morpho-Hessian filter can be summarized as a combination of three main operations: 1) derivative-based directions and tubular-objects detection, 2) field regularization, and 3) spatially-variant closing.

At the first step, the Hessian matrix in its eigen form is produced for each voxel on multiple scales of the image. The choice and number of scales is chosen in order to detect thin objects with different radii. The vesselness value is calculated across these scales according to Eq. 4.6, 4.8, where at the end the maximum response is chosen for each voxel giving the principal direction obtained with the Hessian matrix.

In the second step of the filtering, the directional field propagation is performed as described above.

In the following step, using a family B(e 1 (x)) of centred segments as structuring elements oriented in the direction of e 1 , and of fixed length for the morphological closing operation with the aim of vessel reconnection. The SE-based dilation is followed by the adjunct erosion ε B(δ B ) as in section 4.2.2.7.

conclusion and discussion

In this chapter, we have presented a non-linear, mixed Hessian-based and morphological method in order to enhance linear objects within the scale-space theory. The proposed filter has all the properties of a morphological filter (opening or closing), which ensures it can be used within a morphological processing pipeline with no side effect. This morpho-Hessian filter performs like a reconnecting inverse diffusion filter in some ways.

We have also presented a concise theory of spatially-variant morphology with the main result being an efficient algorithm in O(MN) for spatially-variant openings and closings, where N is the size of the image, and M (M N) is the size of the structuring element. While SVMM is not in itself novel, we believe that the presentation of the algorithm in this thesis is simple and enlightening.

Overall, this approach may be seen as a productive combination of linear and nonlinear techniques. Future work includes studying alternatives to the Hessian tensor for improved noise robustness and direction detection.

As proposed in [START_REF] Dokládal | Grey-scale morphology with spatially-variant rectangles in linear time[END_REF], the size of the spatially-variant structuring elements could be resized according to the eigenvalues of the Hessian matrix.

In addition, it is envisaged to use a variable and more flexible structuring element for morphological directional operations, such as paths instead of segments.

Vascular network detection from angiographic data is an important task, justified by various applications: pathology detection, treatment planning, and more. In this context, 3D MRA provides a detailed visualisation of vascular networks up to the resolution of the generated images (i.e. 0.5mm).

However, the small size and the complexity of the vascular structures, coupled to noise, acquisition artefacts, and blood signal heterogeneity (especially signal discontinuity) make the analysis of such data a hard task, thus justifying intensive efforts devoted to its segmentation.

Here, we especially focus on cerebral angiographic image filtering. Our final purpose is in particular to segment vessels to further assist visualization of arteriovenous malformations of cerebral vasculature. As a primary step, we aim at filtering vessels while detecting the smaller ones and correctly reconnecting these to larger vessels. We then segment the vascular structure and perform its skeletonization for further network topology analysis and visualization assistance.

motivation

Blood interaction takes place at different scales, through which large-scale flow is connected to cellular and sub-cellular biology. Blood vessels deviate and bifurcate while connecting with each other into networks. In the brain, as in all organs, these networks belong to two principal classes: arterial and venous. The arterial networks are the dominant vascular structures consisting of three trees connected to the Circle of Willis: the left and right carotids, and the basilar artery. From the Circle of Willis, the vessels are formed in a number of sub-trees (their actual number may differ from one individual to another). Around the inside of the skull, various venous trees are present, which are connected to the arterial vessels through capillary beds.

As a rule, the cerebral vessels are complex, tortuous, vary widely in geometry, morphology and topology even in normal cases, and even more so in pathological ones. Blood flows in these trees in a directed fashion: in the arteries from the parent to the child. An arterial subtree is, therefore, a smaller tree originating from and containing any branch, its children and all of its children's children to the end of descent. Conversely, blood flows in venous trees in the opposite way, from child to parent.

In the case of the arterio-venous malformation (AVM), some arterial vessels are directly connected to some venous ones and locally bypass the capillary beds. In simple cases, there is one feeder and possibly more than one draining vein. In more complex case, multiple compartments can be present, including several separated feeders and draining veins. They can even form so-called hidden compartments that are invisible on angiography compartments. Hidden compartments can only be detected with serial super-selective digital subtraction angiography or serial high-resolution MR angiography [START_REF] Moftakhar | Cerebral arteriovenous malformations. Part 2: physiology[END_REF]. Generally, it is very difficult to mentally picture the 3D intra-cerebral vascular tree.

Moreover, unfortunately, conventional imaging techniques do not supply sufficient information. Digital subtraction angiograms (DSA) provide useful information, but due to projection overlap, it is difficult to determine the location, connections and the vessel parent-child relationship. Magnetic resonance 3D are useful for such information, but images are noisy and still do not provide sufficient connectivity information directly, due to lack of resolution.

Indeed, in clinical MRA data acquired using a standard 1.5T main field, the resolution is so low that even in non-pathological cases some veins and arteries of the head seem to be connected. This is an important problem in disease recognition. The medical specialists can distinguish pathological and normal cases using their anatomical expertise. In gadolinium-injection MR angiography, during the on-line visualization of the imaging, in normal cases, arterial vessels are visualized first, and only then, as the gadolinium injection passes together with the oxygen in the blood, the venous vessels are visualised. However, in the presence of AVM, one or more veins may be visualised before all the others. But as the process takes only fractions of a second in some cases, it might go unnoticed by a clinician, and later this might even more difficult to recognise among the huge vascular structure data. An example of the MRA image with the AVM case can be seen in Fig. 37.

For neurosurgeons, it is critical to visualize correctly AVMs in 3D. The localization of the subtree of the corresponding artery can help understanding the actual problem better. Specialised visualization procedures can help guiding tailored surgical and endovascular treatments. Computer methods of such pathologies perform multi-scale processing. In this application, we consider the macro-vascular network consisting of large arteries, down to a diameter of 0.5 mm, which are patient-specific and can be reconstructed from clinical imaging. An example of such networks is illustrated in Figure 1.

To help with AVMs diagnosis and treatment, an ideal software suite should include the following features:

1. To correctly represent all the cerebral vasculature.

2.

To be able to show the original image and their projection views.

3. To show the selected and/or all vessels related to the AVM nidus.

4.

To be able to label the arteries and veins, 5. To be able to accept interactive corrections for this, for instance :

-To be able to show children of any selected vessel branch.

-To be able to hide certain tree or subtree to improve visualization.

-etc 6. That is fast and integrable with the physicians acquisition sources, image and patients database and other software suites.

This work is not aimed at providing such an ideal software package, but it is motivated by the completion of some of its points in order to propose such a clinical tool in future.

The process starts with the detection of the vascular structure and its enhancement with the morpho-Hessian filter. Then, with the segmentation and skeletonization procedures, a computer model of the vessel trees is built which can be used in study of vascular disease and their visualization. The main purpose of this work is to aid clinicians in image studies by providing enhanced, more precise data and its visualization for medical decisions. This involves multi-disciplinary communication. So, a complementary goal of the work is to establish links and identify important and relevant medical problems.

In the first section we describe these procedures. Then, their performance is evaluated on synthetic vessel image while comparing them with other similar techniques. In the following section, the same methods are applied to clinical angiographic data.

Although, the comparison of the methods by Sato and Frangi is not our final purpose (as it has already been performed, i.e. [START_REF] Olabarriaga | Evaluation of hessian-based filters to enhance the axis of coronary arteries in ct images[END_REF]), we are evaluating their performance according to their parameters and in combination with the spatially-variant morphological filter.

vessel filtering

Our first stage is a detection and filtering step. We apply the morpho-Hessian filter described in the chapter 4 together with the algorithm described in section 4.4. The main goal of this step is the retrieval of the smallest low-intensity vessels and their reconnection.

As described before, this process is a non-linear three-dimensional spatial filtering technique based on second-order partial derivatives of Gaussian and mathematical morphology. The filter first distinguishes the vessel-like objects and their the local orientation with the second order derivatives grouped in the Hessian matrix form. It then performs a spatially-variant morphological closing (assuming bright vessels on a dark background) according to these local directions.

vessel segmentation

Following the morpho-Hessian filter -which performs like a reconnecting inverse diffusion filter in some ways -in order to segment the vascular structure we perform a region growing method.

Seeded region growing

Seeded region growing (SRG) [AB94, MJ97, BT99] is a well-established method for the segmentation of objects with implicit contours in grey-scale or colour images, based on a set of markers.

SRG takes two images as input: a control image I and a seed image S. I can be virtually anything but in this case we will take the example of a single discrete grey-level image, on a nearest-neighbour connected grid (e.g. the 6-connected graph in the 3D cubic grid). S contains a collection of labelled binary regions, S i ⊂ S, i ∈ 1, 2, 3, .... Regions with the same label do not have to be connected.

Let U be the set of pixels unassigned to any region but connected to at least one of them: U = x ∈ S, ∀i, x ∈ S i , N(x) ∩ ∪ i S i = ∅, where N(x) is the immediate neighborhood of x. At each step of the algorithm, all the points x ∈ U are examined in turn and a distance δ(x, S j ) is computed between x and all the S j it is neighboring to, i.e: S j , N(x) ∩ S j = ∅. The pixels y that possess the minimum distance are assigned to their neighboring region, and the process is repeated until all the pixels are associated with a region.

Typically, the distance δ can be defined as follows:

δ(x, S j ) = I(x) -mean y∈S j (I(y)) ,

where I(x) is the value in I at x. A simple measure for δ is preferred, as complex ones like the Mahalanobis distance may require non-local updating rules.

Several algorithms exist to implement SRG. The original Adams-Bischof algorithm (ABA) [START_REF] Adams | Seeded region growing[END_REF] uses a simple linked list of pixels ordered with respect to δ. Pixels found around the boundary of each seed are inserted into the linked list at a position dependent on δ. The priority of each pixel is never updated. Pixels are simply processed as they appear at the top of the list. Once processed, a pixel is simply removed from the list.

The Mehnert-Jackway algorithm (MJA) [START_REF] Mehnert | An improved seeded region algorithm[END_REF] is significantly more complicated. It uses a main ordered queue of LIFOs to determine pixel priority and a secondary LIFO queue for use during pixel labelling. All the pixels at the same level of priority are processed at the same time. Pixels with the same priority are added to a growing region in parallel. Pixels that are at the boundary of two regions with identical distance with either (ties) are treated separately. Finally the pixels that are being labelled and their unmarked or unlabelled neighbours have their priority recomputed, and they are re-introduced in the queue if it turns out to be different from their original priority. Thus pixels can be present more than once in the main ordered queue.

Finally, the Beare-Talbot algorithm [BT99] uses a series of separate stable ordered queue structures (so-called splay-queues [START_REF] Sleador | Self-adjusting binary search trees[END_REF]) for each seeds, in order to update pixel priority statistics as the regions grow, as well as to process ties fairly. The authors show that their algorithm retains the speed of ABA with none of the side effect of MJA, particularly on 3D images.

SRG bears some resemblance to the watershed algorithm, but does not require the computation of a gradient, which on a noisy image is useful. It also takes into account region statistics, which induces a degree of flexibility to the algorithm. SRG with the BTA follows the mathematical description above, but does not segment the image into regions with minimum distance from the statistics of the seeds. This latter problem is NP-hard, whereas SRG is a greedy algorithm with O(N. log(N)) complexity, which N the number of pixels in the image. However SRG is known to perform well in practice, due to its region-competitive characteristics.

In our work, we used SRG with the Beare-Talbot algorithm, as it was shown to perform best on 3D images.

vascular network analysis

Our vascular network analysis algorithm operates in two major steps: (1) thinning of the binary input vasculature for centerline extraction using the efficient parallel algorithm of [BC10]; (2) identification of junction points (vertices) and segment centerlines (edges).

At the first step we construct a skeleton of the whole cerebral vasculature. Here we use a "curvilinear" asymmetric 1D-isthmus-based skeleton as described in [START_REF] Bertrand | A new 3d parallel thinning scheme based on critical kernels[END_REF]. A "curvilinear" element may be characterized by the fact that, if we remove it, we "break" (locally) the object. The voxels of such a 1D isthmus are characterized in cubical complex space in order to be defined as 1D isthmus voxel. By using the asymmetric algorithm, this method guarantees a minimal skeleton. See Fig. 38 for an example of 1D isthmus skeleton and an asymmetric 1D isthmus skeleton.

The skeleton is then decomposed into simple and end points, junctions and segments (see 3.2.1). This way the complex vascular network is decomposed into individual vessel segments suitable for topological and quantitative analysis. Here, a segment is a piece of vasculature connecting two vessel junctions. This mapping of segments and junctions can be represented by a mathematical graph. However, due to their discrete nature, thinning algorithms do not produce smooth curve-skeletons. Boundary irregularities propagate all the way to the curve-skeleton during the thinning process. Although some algorithms may produce smoother curveskeletons than others, smoothing can be performed in a post-processing step, regardless of the extraction algorithm used to compute the initial curve-skeleton. We have chosen to convert the skeleton to a graph and later to perform a smoothing of this graph for a better visual aspect. AN example of a skeleton produced with the asymmetric 1D isthmus algorithm of an object and the graph-skeleton can be seen in Fig. 39.

For each segment neighbouring junction candidates were grouped into one junction and stored as vertex of the segment graph. For each junction and each segment starting point, the segment centreline is followed voxel by voxel until the next junction is reached. The segment along with the centreline is stored as new edge in the segment graph. Spurious leaf segments shorter than a user-defined minimum segment length are pruned, remaining segments are merged where appropriate.

In order to avoid topological holes (see section 3.2), which create cycles in the skeleton, and prior to skeletonization, the segmented image is hole-closed using the algorithm described in [START_REF] Aktouf | A three-dimensional holes closing algorithm[END_REF].

In order to compensate for the occasional disconnected noise and spurious skeleton branches, a certain number of simple points can be deleted [START_REF] Talbot | Euclidean skeletons and conditional bisectors[END_REF].

algorithm

In the end, the algorithm of the cerebral angiography image analysis can be summarized as a combination of the following operations:

methods evaluation

An important issue in medical image analysis is the validation of the proposed methods, which for most algorithms is based on at least one of the following items:

1. Comparing the results with specialists' expertise; the extracted centerlines are evaluated either visually or by calculating distances between tracked and userindicated centerlines, which are considered as the ground truth.

2. Constructing physical phantoms with known characteristics, taking their image and comparing the results of processed images with the known features; or 3. Using synthetic images.

While the first evaluation is desirable as expert opinion can carry a lot of weight, in reality a lack of objective criteria as well as inter-observer variability are highly likely. In addition, for 3D cerebral angiographic data, this represent a considerable amount of work. To our knowledge, this has not yet been performed. The second method requires machining resources as well as access to a full imaging system to evaluate the method. The inherent complexity and size constraints of cerebral vascular networks makes the use of this validation strategies practically infeasible. The only remaining option is to use 3D synthetic images of tubular structures that resemble the features of the real vessels.

In this work, the results were validated at first on synthetic data and then the same methods are applied on clinical data and discussed.

Evaluation on synthetic data

We used the synthetic dataset from [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF], which is a 100 × 100 × 100 voxels, isotropic dataset containing a tortuous, branching vessel-like object of varying radii (0.5 to 4 voxels), which does not simulate a specific anatomic structure. A rendering of the surface of the object is shown in Fig. 40. The object's cross-sectional intensities have a parabolic profile ranging from 150 at the object's edge to 200 at the object's midline, and the background intensity is 100, which makes it a typical profile of contrast MRA and TOF MRA for small vessels.

Two different types of noise are present in this data. First, the discrete sampling of the continuous object produces errors. Second, intensity noise added to the data, namely: Gaussian noise with standard deviations of k equals 10, 20, 40, and 80. The k = 20 data is representative of the noise level in MR and CT data. The k = 40 data more closely resembles the noise magnitude of ultrasound data. The k = 80 data was chosen to explore the methods' performance on the worst case data which is hardly possibly in clinical applications. Slices of these synthesized images with four levels of noise can be seen in Fig. 41. Hereon, we refer to these synthesized images as phantoms.

Based on the noise-less phantom image and its centerline, Receiver Operating Characteristics (ROC) curves [START_REF] Metz | Basic principles of roc analysis[END_REF] are computed to assess the detection performance of the methods. The ROC curve plots the probability of true positive (TP) detections versus the

Accuracy varies between values 0 -1.

All the methods were tested on the original phantom image and images with four levels of noise. In the case where the resulting image is gray-scale, it is thresholded from its minimum to maximum in small increment steps. In the case of binary result (e.g., segmentation, skeletonization), the ROC curve will contain only one point.

Filtering results

As a filtering method we have used the morpho-Hessian filter. However, in this section we compare only the performance of the vesselness measure by Frangi and Sato associated with this filtering. The thresholded filtering results are compared to the ground truth. We use both normal and hysteresis thresholding methods.

For experiments with Frangi's vesselness, its parameters have been set to: α = β = 0.25 and γ = 5.

In general, the lower scales are more sensitive to noise, while bigger ones tend to join tangent structures. Moreover, in [SWB + 00] it has been reported that the vesselness response curve differs by not more than 3 percent between σ = 1 and σ = 2 √ 2 using a set of discrete values

σ i ∈ 1, √ 2, 2, 2 √ 2.
For our filtering method, in order to evaluate the scale-space performance, we use four discrete scales between σ = 1 and σ = 4. We also perform a multi-scale filtering between scales σ = [1 -4] with four natural logarithmic incrementation steps. For Sato's vesselness, we have used the parameters α 1 = 0.5 and α 2 = 2. The same scale parameters were used.

The vesselness results is a probability function with values from 0 to 1. It was thresholded with all possible values within this range with a step size of 0.1. The performance results of two vesselness measures and their two thresholds are calculated on the ROC accuracies and averaged for each image. The plot of these results can be seen in Fig. 42. We notice that Sato's vesselness shows better results on all images. Moreover, the hysteresis thresholding gives better results for both vesselness functions while giving a higher improvement for Frangi's vesselness.

We have tested Frangi's vesselness both normalized by the largest eigenvalue and not normalized. These segmentation results are obtained with simple thresholding. In Fig. 43, Frangi's vesselness both normalized and not, as well as Sato's measure are shown for the image with Gaussian distribution 10. We can see that at the first scale, all the methods respond to the noise as well as to the phantom. It should be noticed that for the non-normalized vesselness response the bifurcation point is not well detected. It is important to notice that the bifurcation seems to be poorly detected because of the automatic threshold result chosen according to ROC curves, which is a compromise between true positive rate and false positive rate. This means that this bifurcation is detected well at some intensity value where there are many false positives.

At the second scale, the results are similar for normalized Frangi's and Sato's vesselness, while non-normalized Frangi is noisy and with the bifurcation is not well marked.

However, at the third and following scales, Frangi's non-normalized vesselness is the only one that detects the bifurcation even though it responds clearly to the noise and, as a result, the surface of the segmentation is more irregular. Therefore, for further tests we proceed only using Frangi's non-normalized vesselness response.

By using hysteresis thresholding for the three aforementioned vesselness measures, we get rid of the noisy responses while increasing the accuracy. The results can be seen in Figure When looking at the plot of the same results in relation to the average of scales (Fig. 45), more regular, almost increasing results are achieved by Sato's vesselness giving their maximal result with multi-scale σ = [1 -4]. Frangi's vesselness performs similarly to Sato's on scales 1 and 2, while giving worse results at higher and multiscales. Again, for both measures it can be noticed that the hysteresis threshold improves the results, especially at lower scales, which is due to the fact that at such scales there are more spurious and false responses which can be eliminated by such thresholding. However, here we can also observe the controversy between the accuracy numbers and the vistual results: while Frangi's vesselness detects an important part of the phantombifurcation -it enlarges the object and by doing so it achieves lower accuracy values.

Segmentation results

The images are segmented with the seeded region growing (see section 5.3.1) algorithm. Region growing is performed on the morpho-Hessian enhanced image in order to test the morpho-Hessian filter and to obtain the image better reconstructed in spite of heavy noise. For the morpho-Hessian filter, as structuring elements, we use a family B(e 1 (x)) of centered segments oriented in the direction of e 1 and of fixed length 11 voxels.

For the seeded region growing, the inside marker is obtained with the same thresholded vesselness result that is used in the field propagation step (see section 4.3). The outside marker is the inverse of the thresholded morpho-Hessian filter result. Threshold is chosen according to the original image allowing for a slight over-segmentation, which is augmented at the end by the fact that the morpho-Hessian filter increases intensities of the vessels. In such a way, we are ensuring no useful information is lost. Due to the result of the previous step, we have used the hysteresis thresholded results of the Sato's and Frangi's vesselness.

In order to compare the performance of our morpho-Hessian filter, we also perform region growing on the anisotropic diffusion enhanced image [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF]. The external marker is the inverse of the thresholded result of this filter. The internal marker is the Sato's vesselness thresholded result, the same as for the SV-closed region growing.

We also perform region growing on the original phantom images, the internal marker as in the previous case. And the external marker is the inverse of the threshold of the image.

So, we compare four different results of the region growing operation:

1. morpho-Hessian filtered based on Frangi's vesselness hysteresis thresholded (in the following referred to as "Frangi's vesselness").

2. Morpho-Hessian filtered based on Sato's vesselness hysteresis thresholded (further referred to as "Sato's vesselness").

Anisotropic diffusion filtered image.

4. Original image. In Fig. 46 the average ROC accuracies can be seen for region growing segmentation results with the four strategies described above. From this graph it can be noticed that all the methods except vesselness performed perfectly for the original phantom image. For the noise-corrupted images, with the levels 10 and 20 Frangi's vesselness with morpho-Hessian outperforms Sato's vesselness, and it yields only slightly better results on the higher-noised images. Moreover, it can be noticed that the region growing results on all the images except Frangi's vesselness with SV-closing look very similar.

The results of the region growing segmentation results of the morpho-Hessian enhanced images based on Frangi's best performance (σ = 2) on the original phantom image and the Sato's vesselness both mapped on the ground truth can be seen in Fig. 47. The over-segmentation results can be seen (in orange), which, once again, can be explained by the fact that Frangi's vesselness slightly enlarges the object and SV-closing according to its directions also slightly dilates the original image.

Average ROC accuracies of the segmentation results based on four different methods in relation to scales σ = {1, 2, 3, 4, [1 -4]} can be seen in Fig. 48. There we can seen that Frangi's vesselness based result has better ROC accuracies than other methods up to scale 3. It shows worse accuracies on the higher scale and multi-scale, which is again questionnable observing its multi-scale result on the image with k = 10 as in Fig. 49.

Skeletons

The results of the region-growing segmentation with four strategies were skeletonised according to the algorithm 5.4. According to these results, the skeleton of the morpho-Hessian filtered image based on the region growing with Frangi's vesselness markers performed the best in three images except the original one and with noise k = 80 (see Fig. 50). The worst performance for the original data can be explained by the fact that in comparison with other methods, Frangi's vesselness enlarges the object and might give incorrect directions which influences the performance of the SV-closing. Holes smaller than 3 pixels were filled. Simple points belonging to branches of length smaller than 8 were deleted in order to avoid spurious responses. However, by doing so, the length of each branch is also decreasing. In addition, if the junction point was not well detected, it will be missing even more points. This operation might be useful in practice, but it influences accuracy results (as can be seen in Fig. 51).

The comparison of skeletons in relation to scales can be seen in Fig. 52. From this plot, we can see that all the methods have performed similarly at each scale. The skeleton based on Frangi's vesselness outperforms 2, 3 and [1 -4]. The morpho-Hessian with vesselness outperforms the original image-and the anisotropic diffusion-based skeletons.

Fig. 53 shows the results of the skeletons based on the morpho-Hessian with Frangi's vesselness and anisotropic diffusion (segmented with region growing) for all the noiseadded images with k = {10, 20, 40, 80} and at respective scales σ = {1, 2, 3, [1 -4]} (this figure presents only few of the results). These images exhibit the fact that the two strategies perform similarly at scales 1 and 2. However, at scale 3 and [1 -4] the morpho-Hessian detects more of an object even though with the pruning step some pixels are deleted. In any case, the image with k = 80 is a very difficult case for which the methods perform surprisingly well. Here, we are interested in the enhancement of smaller vessels (in particular their reconnection), then segmentation of the whole vascular network, and its skeletonisation for further vascular network analysis. In particular, we are interested in developing an aid to clinicians with the visualization and analysis of arterio-venous malformations. Our particular interest with regards to reconnection is that disconnected networks do not lend themselves well to topology analysis.

To illustrate the performance and its pros and cons on the real vascular data, we apply the algorithm described in section 5.5 on an excerpt of the vascular structure of the smallest visible size in our image data. For the vesselness measures we have used the same parameters as for the synthetic images, except scales are σ = [0.5 -2], as the smallest vessel size in such images is 0.5 mm. The length of the structuring element is the same, And no pruning is used in the skeletonization step.

The results for all the methods are illustrated in Fig. 54. The maximal intensity projection (Figure 54a) and its surface-rendered image (Figure 54b) show that it is densely disrupted by noise. In Fig. 54d we can see that the directions are regularized in comparison with the initial Hessian detections (54c). However, we can notice that due to noise, irregularity of the surface (due to the resolution) and the digitalization, the directions are still not perfectly corresponding to the object surface rendering.

In Figure 54f, it can be observed that the vesselness diffusion filters out the noise while smoothing the surface of vessels. However, doing so, it also removes the smaller tubular objects. In contrast, morpho-Hessian filter (Fig. 54e) keeps the same rough vessel surface as the original surface rendering, but it also preserves and reconnects the smaller objects and eliminates the noise. In 54g (in blue) we can see the result obtained with the seeded region growing of the morpho-Hessian enhanced image and in 54f the result of the thresholded anisotropic diffusion. In comparison with the anisotropic diffusion, we can observe that we have performed the reconnection of the object and got rid of the noise. However, we have also thickened the vessels somewhat. While the anisotropic diffusion preserved much of the noise, it also preserved the topology of the object while noticeably smoothing the surface of the vessels.

When looking at the skeletonization results according to the original image (Fig. 54g), we can see that we joined some small parts of the vessels, due to scale-space issues, further filtering and segmentation. At the same time, the smaller branches were produced as an effect of the irregular surface which can be pruned but with this step we risk to lose some important details. The skeleton transformed into a graph can be seen in Fig. 54h.

We have compared the vesselness measures of Frangi and Sato for the whole-head image data. The result of the thresholding of the two can be seen in Fig. 55. Frangi's vesselness is less smooth than Sato's and produces more responses to noise. Still, we have selected it for further processing as we are interested in reconnecting all the smallest parts and vesselness can serve effectively as a marker.

Arteriovenous malformations segmentation and visualization

As we are interesting in assisting the clinical doctors in the assessment of AVMs, here, we apply the above procedures to the image with an example such malformation.

The example of AVM can be seen in Fig. 37, indicated by the arrow. However, even though this image topologically has all the properties of an AVM, it is in reality in this case a normal feature. The direct apparent connection between arteries and veins is caused by noise and lack of resolution in the IRM data. We stress that we received the dataset as a potential AVM with the feature identified by an arrow as in the figure. Nevertheless, for simplicity, we refer further to this as AVM. In addition, this example shows how difficult it is to distinguish pathological from non-pathological cases from the clinical data. We apply our algorithm as described in section 5.5 on the whole-head image data with the AVM. The image is analysed on scales σ = [0.5 -4]. Its segmentation result can be seen in Fig. 56 (on the left) with the skeletonization result (on the right in blue).

We have simulated the user-guided selection of AVM by choosing its connected component of the thresholded vesselness image. The selected AVM can be seen in Fig. 56 together with the skeleton of the whole head.

As our aim is analysis of the AVM, we have performed the grey-level reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] using the selected connected component of the AVM as marker, and the result of the whole-volume segmentation as mask with finite number of iterations. This has provided us with the bigger part of the AVM showing the connection to both identified veins and arteries. The result can be seen in Fig. 57 together with the maximal intensity projections (MIP) of the image.

We have obtained the skeleton of this larger AVM (Fig. 57b). Simulating the procedure of AVM analysis, we further enlarged the initial AVM to connect it to more vessels (Fig. 57c) and constructing its skeleton (Fig. 57d). The skeleton helps with viewing the connections between individual networks and understanding their topology. This was demonstrated to us when we ran the same procedure by mistake on the other side of the head and found the same structure, which should have been unsurprising, given that on this dataset the structure is in fact normal.

To illustrate the usefulness of such visualization methods in clinical applications, we demonstrate a real example of an AVM. An image of an operated AVM is demonstrated in Fig. 58a. The darker zone (indicated by the green arrow) is the removed brain tissue with the vascular structure. However, with time the pathology is persisting and arteries and veins are still connected. Fig. 58b shows the AVM where it was operated. In Fig. 58c there are pointed with arrows the places of connection of the AVM to veins and arteries.

conclusion and discussion

In this chapter, we have applied a combined morpho-Hessian filter for detecting, enhancing and reconnecting blood vessels. We used the seeded region growing method in order to segment the vascular structure. Then, we skeletonized this result with the method by Couprie and Bertrand in order to obtain a 1D representation of the vascular structure which can help in assessment of vascular disease, more particularly, arterio-venous malformations.

Our algorithm acts on several levels of the image data: voxel-and neighbourhoodbased filtering -giving the precision of the detection, segmentation based on regions information -semi-local representation of the image, and vascular tree -global topological and geometrical representation of the object. At the voxel level, only low-level appearance hypotheses are exploited with the second order derivatives and morphological closing operation. Assumptions of local tubularity are expressed naturally with the vesselness function, while seeded region growing takes in consideration the intensity information. The tree skeleton beholds bifurcations and centerline regularity.

The proposed algorithm was evaluated on synthetic image data representing a vascular bifurcating object with different levels of noise and the real clinical angiography data. The use of the morpho-Hessian filtering was compared via segmentation and skeletonization results with anisotropic diffusion filter [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] and initial image. These evaluations have shown the capacity of the proposed approach to reconnect the smaller parts of the vessel network with the larger ones while removing noise.

As previously said, most often, clinicians visualize 3D medical data slice by slice. This, however, does not facilitate the task of understanding the connections of vessels, which is sometimes very tortuous in 3D. The skeleton representation of vessels used in this work can help clinicians for vessel following from one slice to another, featuring the best of two worlds: their usual 2D environment and 1D representation of vessels presenting within themselves the topology information. By using the asymmetric 1D isthmus algorithm, we have guaranteed a thin skeleton. It has to be noticed that some parts of the approach could however be improved. First, the filter is based on the directions computed from the Hessian matrix, therefore it strongly relies on the quality of these computed directions. The direction computation could actually be improved by also considering the first-order derivative of the image (as in [START_REF] Bouattour | Vessel enhancement in 2d angiographic images[END_REF]) in addition to second derivatives.

Second, our scale-space approach tends to enlarge objects and to reconnect tangential objects, which especially causes problems in the detection of the smallest objects. This issue is difficult to address. Using only the smallest scales yields spurious responses to noise and responds to edges at larger scales. Only using the larger scales is also insufficient, as smaller objects are likely to disappear altogether. Using a maximal responses selection from smaller to larger vessels is not satisfying, as it acts like a sum of all the responses. A more elaborate analysis should be used, like γ-scale fitting, as proposed by Lindeberg [START_REF] Lindeberg | On Scale Selection for Differential Operators[END_REF]. Alternatively, measure based on both first and second order derivatives can be used.

The commonly-observed issue of incorrect vesselness measures at bifurcation does not cause us much problems as we perform direction field regularization and spatiallyvariant filtering. However, it il less precise for smaller objects.

As far as the efficiency issue is concerned, a closed-form solutions of Hessian matrix could be calculated instead of the complete ones, as in [START_REF] Orłowski | Efficient computation of Hessian-based enhancement filters for tubular structures in 3D images[END_REF].

Another step forward in facilitating the visualization of vascular trees that was envisioned in this work and left unfinished for the time being is the vascular tree hierarchy calculation with such schemes as Strahler [START_REF] Strahler | Quantitative analysis of watershed geomorphology[END_REF], as presented in section 3.4.1.

Such a hierarchy analysis can also help in handling scales with regard to the topology of the vascular tree. By knowing the level of the tree, experts know approximately what vessel size is considered.

Overall, our approach may be seen as a productive combination of linear and nonlinear techniques. We have shown that the proposed segmentation and skeletonization procedures can significantly aid vascular network visualization and topology study.

acknowledgments

The author would like to thank Nicolas Passat, Michel Couprie and John Chaussard for their contribution with materials, methods and generous insights. We wish as well to thank the following medical personnel from Colmar Hospital: A. Tournade, M. Musacchio, M. Lagneau, H. Oesterlé, who provided medical expertise and data.

Low dose X-ray images are used in interventional radiology therapies for observing insertion of catheters and guidewires in the vascular system of patients. Guidewire detection assists in such applications as interventional navigation, adaptive visualization enhancement, guidewire reconstruction and respiratory motion tracking. To detect guide-wire from such "2D plus time" images is an essential initial step prior to higher level techniques. The quality of the guidewire enhancement algorithm boosts the performance of the final application.

In this chapter, we compare the performance of several line enhancement algorithms for the application of guidewire detection. The purpose of this work is to select a small number of the most promising 2D line detection methods for medical applications.

guidewire detection methods

In X-ray images, guide-wires (GW) appear as thin, dark lines. Therefore, their detection can be done with methods commonly used for thin objects detection (as presented in chapter 2). The low contrast-to-noise ratio of such images is a common property, as shown in Fig. 59. On this representative image, one can see tissues of the heart superimposed to the diaphragm, the lungs, the spine, the ribs and the interventional devices.

While there exist numerous papers on the subject of line detection in noisy images, as discussed in chapter 2, there is relatively little literature devoted to the more specific topic of guidewire tracking, which is rather unexpected given the clinical importance of endovascular interventions.

In [BNM + 00] a tracking of the guide wire is performed with a spline, and then the spline position is optimised numerically using Powell's direction set method. The optimization is designed to deform the spline so that it has the minimal length, remains smooth, and matches the guidewire position in the image.

In a similar work, Slabaugh [START_REF] Slabaugh | Variational guidewire tracking using phase congruency[END_REF] derived analytic equations to evolve the control points of a spline in order for the spline to match the image data.

Kunz [START_REF] Kunz | Line detection in strongly noise-corrupted images[END_REF] has applied watershed-line method to a steerable filtering result on guide-wire and synthetic images showing improved SNR values after the filtering.

However, no quantitave evaluation of the guidewire detection has been performed except a visual one in the aforementioned works. This comparison of line-detection methods is similar to the work of Ayres [START_REF] Ayres | Performance analysis of oriented feature detectors[END_REF] where the steerable, line-kernel filter [START_REF] Freeman | The design and use of steerable filters[END_REF] and Gabor filters [Chu92] are compared. Bismuth [START_REF] Bismuth | A comparison of line enhancement techniques: applications to guide-wire detection and respiratory motion tracking[END_REF] has compared steerable filter, the Hessian based methods and rotated line filters.

In this work we compare two morphological filters especially designed for thin objects detection: rank-opening and path-opening filters, the steerable filter, the structure tensor. These filters are described in the following sections, the Frangi's vesselness filter (described in section 4.1), and our morpho-Hessian filter (as presented in chapter 4).

Steerable filters

Steerable filters are a subgroup of derivative-based methods. One of the best performing classes of such filters are Gaussian derivatives and their linear combinations [START_REF] Freeman | The design and use of steerable filters[END_REF]. They are efficient means to compute filters that can be rotated with a small computational cost.

These filters allow to calculate a filter response for an arbitrary orientation by linear combination of the responses of a set of basis filters. The angular selectivity of the filters is coupled to the number of basis filters, i.e. increased selectivity implies an increased number of filters and accordingly increased computational effort. So far, the question of how to reduce this effort e.g. by appropriate subsampling has not been addressed. Moreover, little is known about how to combine filter responses in order to detect thin straight or curved lines.

Steerable filters are comprise for applying a filter bank according to directions. This is not only more efficient, but is also more precise as you do not need to sample the orientations.

Let f be a function of R 2 to R defining impulse response of filter, Let f θ be th result of the rotation of f by an angle θ. The filter f is said to be steerable [START_REF] Freeman | The design and use of steerable filters[END_REF] if f θ for any θ and can be expressed as a linear combination of f θ i with 0 < i < N:

f θ = i ω i (θ).f θ i (6.1)
As a result convolving with f θ for various θ can be speeded up. Moreover the value of θ maximizing the response to f theta can be derived analitically.

A large class of filters can be steered, 2D directional gaussian derivatives as well as circular windows times a polynomial (in this case a separable steering basis can be built). Freeman defined the conditions for a filter to be steerable and proposed techniques to compute the basis filters and the coefficients. He also studied separable basis for steerable filters. His work encompasses the computation of second derivatives of Gaussian. The Hessian filter is a particular case of steerable filters.

Several teams designed and used steerable filters for elongated structure detection [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF]. Jacob [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] designed a family of oriented filters optimized according to Canny like criteria for the detection of ridges. The steerable filter approach makes sense when the number of basis filters is reasonable. However, when one attempts to increase the SNR of the filter, one must make the filter as elongated as possible, at the cost of a higher number of basis filters. Therefore the steerable approach can loose its efficiency when one wants to design very elongated filters. It is still important to keep in mind that a large class of steerable filters can be decomposed on separable basis that is favorable to efficent implementation.

Among the steerable filters, we used the work of Jacob [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF]. This implementation is based on the optimization of a Canny-like criteria. We use the fourth-order Gaussian derivatives with the parameters M = 4, µ = 1/4 and 8 separable filtres proposed in [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] as a more anisotropic filter (in comparisson for instance with secon-order derivatives) with the highest theoretical SNR. See an illustration of second-order and fourth-order derivative filter detectors in Fig. 60.

Parametric opening

An opening and closing morphological operators were proposed to calculate local orientation information for, respectively, bright and dark structures by Soille in 1998 [START_REF] Soille | Image structure orientation using mathematical morphology[END_REF].

Denoting by f a grey-scale image, γ the morphological opening, and L α,λ is a line segment defined along Bresenham lines of length λ and orientation α, the orientation of a bright object at a give pixel x can be defined as:

Dir λ (f)(x) = {α i | γ L α i ,λ (f)(x) γ L α j ,λ (f)(x), ∀α i = α j }.
(6.2)

Figure 60: Ridge detectors with second-order (left) and fourth-order (right) derivatives.

The bright-object directional signature at a given pixel can be obtained by plotting the normalized opened values versus the orientation of the line segment. This signature can then be used to detect crossing lines, flat zones, etc. Then, the following quantities are defined for each point x of the input image f:

Max λ (f)(x) = {γ L α i ,λ (f)(x) | ∀α i = α j , γ L α i ,λ (f)(x) γ L α j ,λ (f)(x)}, Min λ (f)(x) = {γ L α i ,λ (f)(x) | ∀α i = α j , γ L α i ,λ (f)(x) γ L α j ,λ (f)(x)}, Gdir λ (f)(x) = Max λ (f)(x) -Min λ (f)(x).
Morphological openings and closings are sensitive to noisy pixels. They remove all foreground pixels that cannot be covered by the structuring element. A less restrictive operation can be performed by accepting partial fits, i.e., at least k pixels of the structuring element that have to fit the foreground opening. It is equivalent to the union of morphological openings by the structuring elements B i which are obtained by removing λ -k pixels to a given structuring element B containing λ pixels. A parametric opening is an algebraic opening denoted as:

γ B,k = ∨ i {γ B i ∪ {p 1 , . . . , p λ-k } = B}, (6.3)
where λ > k. B is a line segment of k pixels.

Rank-max opening

A fast implementation is possible as the parametric opening corresponds to a rankmax opening [START_REF] Ronse | Extraction of narrow peaks and ridges in images by combination of local low rank and max filters: implementation and applications to clinical angiography[END_REF][START_REF] Heijmans | Morphological filters for dummies[END_REF], which is a point-wise minimum between the original image f and the dilation by a structuring element B of its rank filtering rho with a mask equal to B and a rank equal to k:

γ B,k = f ∧ δ B,ρ B,k (6.4)
When computing openings (or closings), the non-translation-invariant recursive implementation of min/max filters can be combined so as to result in a TI implementation. An example of parametric openings with a rank parameter controlling the noise sensitivity is illustrated in Figure 61.

A detailed discussion about rank-max filters and their relationships with other filters can be found in [START_REF] Soille | On morphological operators based on rank filters[END_REF].

The choice of the k value in the rank morphological filters should be driven by the amount of noise corrupting the local connectivity of the oriented image patterns. The orientation resolution is constrained by the length of these patterns and can be automatically determined from a multi-scale approach, the scale being defined by the length L of the line segments.

Path openings

Path openings and closings can isolate oriented and linear structures that are of length at least L min and that are respectively brighter or darker than their local neighbourhood. To perform this, these filters assign to each pixel the highest (lowest, respectively) grey level where a path L is formed. Within mathematical morphology theory, this task is based on the study of the connectivity of the pixels belonging to one one given path [START_REF] Heijmans | Path openings and closings[END_REF].

To detect thin elongated object, the morphological approach consists of using narrow structuring elements and applying them in all possible (or many) angles. Assuming the object has a bright intensity compared to darker surroundings, the Supremum of all openings obtained by this SE application the linear object is preserved and all other are removed. When using a rather long structuring element, the disconnected (seemingly) parts of an elongated object will be reconnected.

Let E be a discrete 2D image, a subset of Z 2 . We assume E has an adjancency relation x → y meaning that there is a directed edge going from x to y.

With this adjacency relation, the dilation can be defined as δ

({x}) = {y ∈ E, x → y}. The L-tuple a = (a 1 , a 2 , . . . , a L ) is called a δ-path of length L if a k 1 ∈ δ({a k }) for k = 1, 2, . . . , L -1.
Given a path a ∈ E, σ(a) denotes the set of its elements, i.e. : σ(a 1 , a 2 , . . . , a L ) = {a 1 , a 2 , . . . , a L }.

The set of all δ-paths of length L is denoted by Π L , and the set of δ-paths of length L contained in a subset X of E is denoted by Π L (X).

An opening operator α L (X) can be defined as the union of all δ-paths of length L contained in X: α L (X) = ∪{σ(a), a ∈ Π L (X)} L can be considered as a structuring element size with a flexible (to some extent) shape. The path opening α L is the Supremum of the morphological opening using these paths of size L as structuring elements.

Grey-level image decomposition

The binary operator defined above can be extended for grey-level images by replacing the union with a supremum. The linear-time recursive decomposition of the operator α L has been proposed by Heijmans [START_REF] Heijmans | Path openings and closings[END_REF].

By the threshold decomposition, the grey-scale image can be processed with binary morphological operators.

Here, binary images are redefined as functions of the form b : E → {false, true}. Then, given a grey-scale image ginG, a threshold operator T t : G → B with threshold t, and a binary opening γ B : B → B, there exists a grey-scale opening γ G : G → G such that for all thresholds t, T t • γ G = γ B • T t , where • is the composition operator.

This grey-scale opening may be constructed explicitly by "stacking" the results of the binary opening applied to each threshold of the original image. This stacking assigns openings by lines in the latter algorithm are all running in constant time irrespective of L, but for larger L more orientations need to be explored. Note also that the presented algorithm for the supremum of opening by lines is not the translation-invariant implementation, which would be slower still. The area opening algorithm seems to converge to a constant-time algorithm with low constant. The area parameter was simply L, although k × L with k small and constant (e.g: 3) could have been chosen without significantly affecting the result. Memory demands for these algorithms are all low except the recursive path opening implementation which requires an amount of memory proportional to LN , with N the number of pixels in the image.

We observe that the area opening is the fastest algorithm, but that the presented path opening algorithm comes second and significantly faster than the other two algorithms for most useful values of L. For grey-scale images, the path opening transforms for each threshold, the active pixels whose maximal path length λ [p] has decreased store a point (t, λ [p]) in a linked list. This list is a monotonically decreasing in t and monotonically decreasing in λ [p]. Once computed, it is possible to query this structure with any desired path length to extract the associated grey-scale path opening.

This algorithm can be computed in O(NlogL). The average memory required is O(NlogL).

An example of a path opening and, for comparison, an area opening is illustrated in Fig. 62 In order for path openings and closings to be useful in a context where features of interest are arbitrarily oriented, paths openings and closing should be constructed in the usual way through supremums and infimums (resp.), in the same way as with straight-line structuring elements. However, in order to achieve acceptable levels of isotropy with path operators, far fewer compositions are necessary.

Paths openings with increasing L are also increasingly sensitive to noise, as long paths are more likely to contain noisy pixels. In order to decrease the sensitivity to noise, it is useful to allow a few pixels to be ignored along the structuring paths. For this, the rank-max opening, described in 6.1.2.1, can be helpful. In practice, there does not exist at present a way to efficiently compute the median over a family of paths, and so a decomposition of paths openings like Eq. 6.4 hasn't yet been proposed, however, Appleton in [START_REF] Appleton | Efficient path openings and closings[END_REF] did propose so-called incomplete paths, which have the same effect. However, efficiency and memory requirements are proportional to k, the number of "ignored" pixels in the paths.

Structure tensor

Another popular feature detection tool in image processing is the structure tensor. Traditional structure tensor 6.5 [START_REF] Förstner | A fast operator for detection and precise location of distinct points, corners and centres of circular features[END_REF] is obtained by calculating the first order Gaussian derivatives and in 2D can be expressed as ∇f = ((∂ x G σ ) * f, (∂ y G σ ) * f) T . By taking the product U(x) = ∇f(x)∇f(x) T a 2 × 2 matrix for each pixel with one positive eigenvalue λ 1 (x) = ∇f(x) , with corresponding eigenvector ∇f(x)/ ∇f(x) and one zero eigenvalue. The tensor field U is then convolved with an isotropic Gaussian kernel with scale ρ, W = U * G ρ . The resulting tensor field containing richer information on edges and junctions is defined as:

G p * (∇U σ ∇U T σ ) = G p * (U σ ) x (U σ ) x (U σ ) x (U σ ) y (U σ ) y (U σ ) x (U σ ) y (U σ ) y (6.5)
The eigenvectors of this structural-analysis tensor at a position x show the two principal directions of the gradient vectors in the neighbourhood of this point, and the eigenvalues represent its deviation distribution in the gradient space. σ denotes the area of study the distribution of intensity gradients. By denoting eigenvalues with λ 1 > λ 2 , the principal direction of the line-like structure is indicated by the corresponding eigenvector e 1 . We evaluate line-likeliness as in [START_REF] Bigün | N-folded symmetries by complex moments in gabor space and their application to unsupervised texture segmentation[END_REF][START_REF] Deguchi | Detection and enhancement of line structures in an image by anisotropic diffusion[END_REF]:

S(x, y) = λ 1 -λ 2 λ 1 + λ 2 (6.6)
This measure function acts like a probability function, meaning that when S(x, y) ≈ 1, the corresponding pixels belongs to a line-like structure. Conversely, its probability of affiliation to such a structure decreases as S(x, y) → 0.

As for the second order derivative vesselness functions, structure tensor as well can be computed at multiple scales.

S max (x, y) = max(S(x, y, σ)) with σ min σ σ max

The final segmentation can be obtained by simple or hysteresis thresholding as in Section 4.1.3.

guidewire detection results

For diagnostic purposes, for the X-Ray guide-wire images, the critical importance is on accurate positioning of the guide wire and the catheter with respect to the vascular network, patient motion and the low signal to noise ratio of the images. The limited quality of the image is defined by the constraint of low dose used in fluoroscopy in order to minimize the radiation exposure of the patient and radiologist.

Because of the physics of the X-Rays, the clinical images are transparent. In these images, due to the transparency effect, the tissues of heart superimposed to the diaphragm, the lungs, the spine, the ribs and the interventional devices can be seen. Moreover, the images are dominated by various sources of quantum noise. In the darker regions of the image (where the signal is weak), the electronic noise is significant, which can be assumed Gaussian. Everywhere else, the dominant source of noise is the X-Ray photon shot noise, which is accurately modelled as a Poisson process. To simplify the noise pattern, an image transform was applied called a variance stabilization transform (similar to the Anscombe transform). After this step, the noise becomes Additive White Gaussian, which is stationary.

Here, we have used clinical X-ray images of guidewires with a resolution of 1000 × 1000 pixels. For the current study, we have chosen three images corresponding to three levels of increasing difficulty (shown in the first line of Fig. 67) for illustration.

Evaluation framework

In order to evaluate the performance of the above-describes methods, ground truth was provided by manual delineation by an independent human expert operator.

The ground truth for these images was segmented by the medical specialists and was divided into two parts: 1) the main one and the easier to detect, and 2) the secondary one, on the extremities of the guidewire, often present in the high-noise and contrastcorrupted parts. The main part of the guidewire is the most important and is already sufficiently difficult to segment, therefore we have chosen to concentrate on this section for our evaluation.

Based on the ground truth of these images, Receiver Operating Characteristics (ROC) curves [START_REF] Metz | Basic principles of roc analysis[END_REF] have been computed to assess the detection performance of the methods. The ROC curve plots the probability of true positive (TP) detections versus the probability of false positives (FP) against the ground truth as a function of the detection threshold. Pixels above certain threshold along the ground truth count as detections; others are false positives.

The ROC curve is constructed from false positive rate (FPR) and true positive rate (TPR) where each value of the curve is accuracy as in equation 5.1.

However, due to the thinness of guide wires and presence of many other similar objects, we have chosen to dilate the ground truth of the object and mask the nonrelevant part of the image for the preliminary evaluation of the methods. The ground truth was dilated with a near-isotropic polygonal structuring element of radius 10 in order to create the new ground truth, pixels of the resultant image present in this zone are marked as true positives. To create the mask, the ground truth is dilated with a polygonal SE of radius 40, which means that resulting pixels being in this zone will be considered as false positives, and pixels outside will be ignored. The schematic view of the ground truth is shown in Fig. 63.

All the methods were tested on three images: Image 1, Image 2, Image 3.

The path opening and rank-opening filters have been tested with structuring elements of length 17, 31 and 61. They have also been tested in their original formulation and integrating the noise sensitivity value k that should compensate the amount of noise corrupting the local connectivity of the oriented objects.

Both structure tensor and Frangi's vesselness was tested at different fixed scales σ = 1, 2, 3, 4 and the multi-scale filtering with σ ∈ [1 -4] showed the best performance. Frangi's vesselness sensitivity parameters were: α = 0.25, β = 0.25 and γ = 5.

We have also tested our morpho-Hessian filter with the same parameters as for Frangi's vesselness. However, as it is the enhancement filter, producing the image similar to the original one only with the enhanced intensities at the elongated objects, we had to apply another method in order to detect these objects. A natural choice was to use Frangi's vesselness detection scheme again with the same parameters.

The final results of these filters are grey-scale images with maximal values assigned to the detected elongated objects. Therefore, the detection results were evaluated after the thresholding of these results. In addition, the methods were tested with hysteresis thresholding in order to improve the results by the neighbouring responses. The threshold ranges are not the same from one method to another. q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0. Best results: 0.932 0.919 0.911 q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0. Best results: 0.959 0.939 0.935 q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.

Results

The selection of the best performing filters according to their accuracies as in Eq. 5.1 is shown in plot 64 1 . According to this graph, the largest bar represents the best result (max = 3), while the possible accuracy result for one image is between 0 and 1.

The best result with the ROC accuracy measure achieved with the rank-opening (SE length is 31 pxls, width is 7 pxls, noise tolerance r = 10%) can be seen in Figure 65 for Image 1. The best segmentation results for each method and image are shown in Figure 67.

In general, we can note that we have two types of challenges:

1. For each method its parameters are indicated, e.g. "Rank-opening SE = 31 7 , r = 0.1" means that the used structuring element is of length 31 and width 7 and r is the noise tolerance parameter. 116 guidewire detection q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0. -ability to detect the linear structure on a flat noisy background -ability to ignore shorter elongated structures, blobs and other noise. The morphological path opening and rank filters were successful on both counts. We can explain their success on these noisy images by their strongly elongated SE. Such long filters and integration of the noise sensitivity parameter allow to significantly reduce the influence of noise along the structures of interest and permitting the reconnection of the guidewire parts.

For longer SE (e.g. 31, 61), the noise tolerance parameter for path opening was set lower than for rank-opening for speed-up reasons.

Scale-space methods as vesselness and structure tensor cannot be tuned in the same manner, as the choice of the scale depends on the width of the object and not the length in order to measure the contrast across the object. Moreover, the image is blurred on this scale, so that a too large scale would erase the important parts. However, the Frangi's vesselness performs well at detecting the object, but produces many false positives (see its ROC curve in Fig. 66a)2 . However, those objects are also thin elongated structures. Moreover, it detects objects with higher curvatures than path opening methods. These properties can be of interest in other applications, where mostly all the objects with high curvature are prevalent.

Morpho-Hessian filter has not given good accuracy results (Fig. 66b), however, when observing the ROC curve, it has improved the TP rate in comparison with Frangi's performance. But it has also added some false negatives, which can be explained firstly by incorrectly detected objects directions (often orthogonal to the principal one), which in itself is due to the highly noise-corrupted nature of the image. Then, the morphological closing enlarged the object in width. Afterwards, by applying the second vesselness filter for object's detection even more smoothness was added to the image, which increased the FPR.

Finally, in our implementation, the steerable filters do not give good results at all, due to not very anisotropic filter and no additional blobness characteristic as in Frangi's vesselness.

On the tested images, hysteresis thresholding has improved only about 1 -2% for morphological methods and steerable filter in comparison with simple thresholding. Against the improvement of about 10 -20% for vesselness and structure tensor. This can be explained by the fact that the morphological opening detected only the longest elongated object, while other methods detected other present thin objects which are not related to the GW.

conclusion and discussion

Through ROC curve accuracy detection assessment we demonstrated that the morphological rank-and path-opening filters are the most suited for detecting low SNR and low curvature devices. On the other hand, the second order Gaussian derivative filters are more robust to strong curvature of the devices. While the structure tensor and steerable filter approaches showed less interesting detection capabilities.

Moreover, morpho-Hessian filter presents a limitation of its own in that the artefacts present in the original data are not suppressed. If these artefacts have an apparent structure similar to a plane or line, they will be enhanced as well. This is the case with transparency and motion artefacts, both of which create ghost artefacts. There, this method is not well suited for the enhancement of such guidewire images as the directions detected by the Hessian filter were disturbed by the presence of noise and due to the inherently local nature of the Hessian filter. The results can be improved if the directions are detected with some other filter, i.e. as described in section 6.1.2.

For steerable filters, beyond the simple idea of having an elongated filter to reach high SNR, one can go further in the idea of combining responses of fast filters in a local neighborhood. Ideally, one would like to be able to combine responses of filters to re-inforce each other, for instance along all possible "admissible" curves and retain at each pixel the most relevant one. Local minimal path, path opening and tensor voting all tend to do this.

There is a wide scope for further work on this topic, as it has much applicability and relevance. As real-time constraints are imposing that only efficient filters be tested, in the future more complex filters should be of interest, in particular the local shortest paths of Vincent [START_REF] Vincent | Minimal path algorithms for the robust detection of linear features in images[END_REF]. Here, we were mainly interested in detection of one part of the guidewire, but segmentation of the whole of it is another interesting and challenging issue. We have performed experiments based on the morpho-Hessian enhancement results together with the fast-marching method [START_REF] Sethian | Fast marching methods[END_REF] and skeletonization, where we managed to segment the whole image without much of the noise response. This approach is planned to be extensively tested for an article in progress on our real-image database containing of 25 images with 100 datasets temporal image sequences per image and a database of synthetic images simulating the guidewire and noise.

acknowledgements

This work is a contribution for an article in progress together with the following authors: Tankyevych O., Bismuth V., Talbot H., and Dokladal P. I would especially like to thank Vincent Bismuth for the provided images, materials and steerable filter result that were used in this thesis.

C O N C L U S I O N A N D P E R S P E C T I V E S

To emphasize only the beautiful seems to me to be like a mathematical system that only concerns itself with positive numbers.

-Paul Klee

contributions

In this thesis, we discussed the thin objects analysis challenges and their detection methods. More particularly, applied for vascular image analysis medical applications with the main purpose of aiding clinicians in image studies by providing enhanced, more precise data and its visualization for medical decisions. This work has involved multi-disciplinary communication, which helped us to establish stronger links and gain better understanding of important and relevant medical problems.

A first contribution of this thesis is the review of the current methods for thin objects analysis. We have discussed the methods of filtering and segmentation of such objects categorised in mathematical frameworks. We have as well outlined some of the new successful combined methods and suggested that the hybrid strategies are the most promising ones. Such a study helped us in our further choice of methods for implementation and tests.

In the following chapter (Chapter 4), we introduced the hybrid morpho-Hessian method based on the second-order derivatives Gaussian filter and spatially-variant morphological operations for detecting, enhancing and reconnecting such objects along with a segmentation and skeletonization strategies.

In Chapter 5, we have then applied this method for the analysis of cerebral vasculature and arterio-venous malformations from 3D magnetic resonance angiography imaging.

Our algorithm acts at several levels of the image data: voxel-and neighbourhoodbased filtering -giving the precision of the detection, segmentation based on regions information -semi-local representation of the image, and vascular tree -global topological and geometrical representation of the object. At the voxel level, only low-level appearance hypotheses are assumed and exploited with the second order derivatives and morphological closing operation. Assumptions of local tubularity are expressed naturally with the vesselness function, while seeded region growing takes in consideration the intensity information. The tree skeleton contains and represents bifurcations and centerline regularity.

The proposed algorithm was evaluated on synthetic image data representing a vascular bifurcating object with different levels of noise, as well as real clinical angiography data. We have also compared two vesselness functions performances: Sato's and Frangi's. The use of the morpho-Hessian filtering was compared via segmentation and skeletonization results with anisotropic diffusion filter [START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF] and initial image. These evaluations have shown the capacity of the proposed approach to reconnect the smaller parts of the vessel network with the larger ones while removing the noise.

In addition, in Chapter 6, we performed thin-objects detection methods evaluation for an application of guidewire detection from highly noised 2D X-ray interventional images. We compared two morphological filters especially designed for thin objects detection: rank-opening and path-opening filters, together with well-established, more classical approaches: the steerable filter and the structure tensor. Through ROC curve accuracy detection assessment we demonstrated that the morphological rank-and pathopening filters are the most suited for detecting low SNR and low curvature devices. On the other hand, the second order Gaussian derivative filters are more robust to strong curvature of the devices. While the structure tensor and steerable filter approaches showed less interesting detection capabilities.

While this work has not provided a complete software package, it has nevertheless allowed us to propose performing and promising procedures for a future clinical tool. Overall, this approach may be seen as a productive combination of linear and nonlinear techniques. We have shown that the proposed segmentation and skeletonization procedures can aid vascular network visualization and topology study.

During this work, we made careful methodological choices in order to obtain efficient and novel filtering and segmentation methods. As described in chapter 8, we have implemented or used the tools with the consideration of its reusability and efficiency. We have developed a specialised testing environment and made extensive tests on the synthetic and clinical data to obtain better knowledge on the performance of the used methods according to quantitative and qualitative evaluations.

perspectives

In the respective discussions of each chapter, we discussed various theoretical and practical possibilities for the improvement of our approaches.

First, for the morpho-Hessian filter. As it is based on the directions computed from the Hessian matrix, its performance strongly relies on the quality of the computed directions. The direction computation could actually be improved by also considering the first-order derivative of the image (as in [START_REF] Bouattour | Vessel enhancement in 2d angiographic images[END_REF]). More recent simpler models were introduced in works by Agam [START_REF] Agam | Vessel tree reconstruction in thoracic CT scans with application to nodule detection[END_REF] with a filter model which is based on the correlation matrix of the regularized gradient vectors (first-order derivatives).

Then, the scale-space approach. We have observed that it tends to enlarge objects and to reconnect tangential objects, which especially causes problems in the detection of the smallest objects. This issue is difficult to address by using only the smallest scales, as Fig. 6. Model-based vessel visualizations with different settings for the smoothing of the vessel diameter. The vessel diameter is color-coded. In the left image, no smoothing is applied which results in obvious discontinuities in the periphery. These discontinuities are due to the limited spatial resolution of the underlying data.

such a system selects an object in a list via its name and the system should provide feedback emphasizing this object. Another obvious reason for the need of 3d emphasis techniques is due to the integration of 2d and 3d views (recall Sect. 4). After an object is selected in a 2d view it should be highlighted in the 3d view to support the connection between the two visualizations. The emphasis of objects in 3d visualizations is difficult by its very nature. From the current viewing position objects might be too small to be recognizable or they might be occluded by other objects. In medical visualizations these problems are prevalent: objects are often concave, are at least partly occluded by others. Therefore simple emphasis techniques such as the use of a special colour for highlighted objects or blinking do not work well. Emphasis techniques therefore should ensure the visibility of the involved objects (see [18]). In principle, two strategies are possible to achieve visibility: the camera position might be changed to make the desired object visible or visualization parameters of occluding objects might be adapted to allow to look through them. The first strategy, in general, cannot be recommended. A radical change of the camera position, not initiated by the user, is often not comfortable because the user has to interpret this change and perhaps dislikes the chosen perspective. Unnatural viewing directions may result. The change of visualization parameters, the second strategy, can be implemented by using transparency: occluding objects are rendered semitransparently to reveal objects behind. This is a viable approach, however showing apparent drawbacks. If several objects are in front of the object to be emphasized all occluding objects must be strongly transparent with the result that they are almost unrecognizable (see Fig. 7). As an alternative, fast silhouette generation algorithms might be used to enhance the visualization of the transparently rendered objects. For the sake of brevity, emphasis techniques could only be touched here. There is a variety of emphasis techniques suitable for medical visualization, including those based on the shadow generation for selected objects (see [15] for a review on such techniques). this gives spurious responses to noise and responds to edges of larger scales. Using only the larger scales is also insufficient, as we will likely lose smaller objects. Maximal responses selection from smaller to larger vessels is not satisfying either, as it acts like a sum of all the responses. A more elaborate analysis should be used, like γ-scale fitting, as proposed by Lindeberg [START_REF] Lindeberg | On Scale Selection for Differential Operators[END_REF], or perhaps measure both based on first and second order derivatives.

For the efficiency issue, the closed-form solutions of Hessian matrix could be calculated instead of the complete ones, as in [START_REF] Orłowski | Efficient computation of Hessian-based enhancement filters for tubular structures in 3D images[END_REF].

When it comes to segmentation, many other techniques discussed in Chapter 2 could be tested. One of the ways that interests us for future research is a user-guided segmentation, for instance using tracking. This is especially interesting as the unavailability of manual intervention can also introduce computational error.

The methods based on the assumption that the cross-section of vascular structures is circular are, of course, too simplified. Moreover, they are not appropriate in vascular diseases diagnosis. However, for many therapy planning tasks it is crucial to understand the spatial relation between pathological structures and adjacent vascular structures. So, it is desirable to be able to choose the visulization fitting the real image data or being more idealized, as it is illustrated in Figure 68. This kind of fitting could also compensate for the irregular surface resuts that our morpho-Hessian shows in comparison with the diffusion filter.

Another step forward can be in facilitating the visualization of vascular trees that was envisioned in this work and unaccomplished at the time of writing, is the vascular tree hierarchy calculation with such schemes as Strahler [START_REF] Strahler | Quantitative analysis of watershed geomorphology[END_REF], as presented in section 3.4.1.

Such a hierarchy analysis can also help in the handling of scales with regard to the topology of the vascular tree. By knowing the level of the tree, experts know approximately what size is expected.

To exploit the skeleton information better, we can extract measurements such as vessel diameter, length, density, torutosity, etc. that could be useful for example in evaluation of pathology severity. Generally, automated vessel size measurements are preferred over visual interpretation because of the inter-and intra-observer variability.

Another way to consider computer-aided analysis of 3D angiographic images is to provide to the human experts a base of high-level anatomical knowledge which can possibly be involved in more specific analysis procedures such as vessel labeling. Such knowledge can in particular be embedded in vascular atlases which are devoted to model qualitative and/or quantitative information related to vessels.

The digital image working environment in the clinical context becomes more complex. As physical and psychological problems such as carpal-tunnel syndrome, visual and physical fatigue become more prevalent, imaging systems need to be evaluated not only with respect to diagnostic accuracy, but also toward the totality of perceptual, cognitive, and environmental factors that contribute to the diagnostic decision-making process [START_REF] Krupinski | Anniversary Paper: Evaluation of medical imaging systems[END_REF].

For the guide-wire detection application, there is a wide scope for further work on this topic, as it has much applicability and relevance. As real-time constraints are imposing that only efficient filters be tested, in the future more complex filters should be of interest, in particular the local shortest paths of Vincent [START_REF] Vincent | Minimal path algorithms for the robust detection of linear features in images[END_REF]. Our approach is planned to be extensively tested for an article in progress on our real-image database containing of 25 images with 100 datasets temporal image sequences per image and a database of synthetic images simulating the guidewire and noise.

1 . 1 Figure 1 3D

 111 Figure 1 3D model of major vessels and bifurcations of the human arterial tree reconstructed from a set of CT, DSA CT and MRA images. Colors represent different parts of the model. Left: Aorta and adjacent arteries. Right top: Cranial arterial network. Right bottom: Carotid artery (source: [GAM + 09]). . . . . . . . . . . . . . Figure 2 (a) Illustration of the normal arterial-capillary-venous relationship. (b) Drawing of the failure of normal capillary genesis, a fistular connection between the developing artery and vein occurs, which causes an arterial shunting of blood with arterialization of the venous outflow (source: [MHMM09]). (c) Illustration of cerebral AVM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3 A fluoroscopic C-arm machine that generates x-rays from one side and photographs them on the other side. . . . . . . . . . . . . . . .

12

 12 Figure 12Results of segmentation of vessels by[START_REF] Chen | Quantifying 3-D vascular structures in MRA images using hybrid pde and geometric deformable models[END_REF]. (a) Reconstructed images for the abdominal arteriography images. (b) Labeled result obtained by directly applying the level set method. (c) MIP of enhanced vessel tree. (d) Interim and final results of the evolution of the propagating surface. (e) The final corrected vessel surface. (f) An enlarged view of part of the final vessel surface triangulated with 3-D Delaunay triangulation. (g) The intersection contours of sagittal, coronal, and axial image planes with the levelset segmented result (white) and the final result obtained after the application of a geometric deformable model segmentation technique (gray). (Illustration is from [CA04]). . . . . . . . . . . . . Figure 13 ROC curve for Frangi method and algorithm of [GAF + 09] (source: [GAF + 09]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 14 3D

 14 Figure143D neighborhoods: 6-connectivity, 18-connectivity and 26-connectivity (from left to right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure15Illustration of the concept of firefront. Fire propagates from the border of the object toward its interior, in an isotropic manner.The point C is a quench point and is the center of a maximal ball (source:[START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure16Illustration of a maximal ball in the Euclidean case. The ball B max is maximal because no other ball can contain it and simultaneously be included in the set. In contrast, the ball B 2 is not maximal (source:[START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure17Illustration of medial axes (source:[START_REF] Blum | A Transformation for Extracting New Descriptors of Shape[END_REF]). . . . . . . . . . . . . . Figure18Discrete medial axis: (a) the silhouette of a bear; (b) its discrete medial axis, with the 4-connected neighborhood as unit ball; and (c) the ultimate erosion with the same unit ball . . . . . . . . . . . Figure19The thinning process on an example 2D shape. Boundary points are marked "B" at the beginning of each iteration and then removed if they are simple (source:[START_REF] Nd Cornea | Curve-skeleton properties, applications, and algorithms[END_REF]). . . . . . . . . . . . . . Figure20Example of a distance transform (right) and its original image (left). Figure21A graph based on the skeleton. . . . . . . . . . . . . . . . . . . . . . Figure22Extraction of junctions and segments from binary image data. The original vascular structure (left) is thinned to its centerline (center), voxels with one or more than 2 neighbors are marked in red as junction points. Segments are identified by following the centerline voxels from one junction to the next one (right). . . . . . Figure23An illustration vascular tree using Strahler Ordering (left), where ( †) indicate examples of trifurcation ordering rules, and ( ‡) indicate junctions upstream of vessel pruning (source: [NBB + 06]). . .

  Figure143D neighborhoods: 6-connectivity, 18-connectivity and 26-connectivity (from left to right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure15Illustration of the concept of firefront. Fire propagates from the border of the object toward its interior, in an isotropic manner.The point C is a quench point and is the center of a maximal ball (source:[START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure16Illustration of a maximal ball in the Euclidean case. The ball B max is maximal because no other ball can contain it and simultaneously be included in the set. In contrast, the ball B 2 is not maximal (source:[START_REF]Mathematical morphology: from theory to applications[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure17Illustration of medial axes (source:[START_REF] Blum | A Transformation for Extracting New Descriptors of Shape[END_REF]). . . . . . . . . . . . . . Figure18Discrete medial axis: (a) the silhouette of a bear; (b) its discrete medial axis, with the 4-connected neighborhood as unit ball; and (c) the ultimate erosion with the same unit ball . . . . . . . . . . . Figure19The thinning process on an example 2D shape. Boundary points are marked "B" at the beginning of each iteration and then removed if they are simple (source:[START_REF] Nd Cornea | Curve-skeleton properties, applications, and algorithms[END_REF]). . . . . . . . . . . . . . Figure20Example of a distance transform (right) and its original image (left). Figure21A graph based on the skeleton. . . . . . . . . . . . . . . . . . . . . . Figure22Extraction of junctions and segments from binary image data. The original vascular structure (left) is thinned to its centerline (center), voxels with one or more than 2 neighbors are marked in red as junction points. Segments are identified by following the centerline voxels from one junction to the next one (right). . . . . . Figure23An illustration vascular tree using Strahler Ordering (left), where ( †) indicate examples of trifurcation ordering rules, and ( ‡) indicate junctions upstream of vessel pruning (source: [NBB + 06]). . .

[ 1 -

 1 Figure 37Slice of a MRA data of the head (left) and a close-up of the AVM. . . Figure38An object and its symmetric (left) and asymmetric 1D isthmus skeleton (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure39An object and its asymmetric 1D isthmus skeleton (left) and its graph based on the skeleton (right) . . . . . . . . . . . . . . . . . . . Figure40Rendering of the synthetic vessel-like object with the added Gaussian noise of κ = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 41 Slices of phantom images with four levels of Gaussian noise (from left to right): 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . Figure 42 Comparison of average ROC accuracies for Frangi's or Sato's vesselness and their respective hysterysis thresholding for original phantom image and with four levels of Gaussian noise: 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 43 From left to right: Morpho Hessian with Frangi's normalized and non-normalized vesselness or Sato's measure on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4]. . . . . . . . . . . . . . Figure 44 From left to right: Frangi's and Sato's vesselness on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4]. . . . . . . . . . . . . . Figure 45 Comparison of average ROC accuracies for Frangi's or Sato's vesselness according to results at four discrete scales and one multi-scale (σ: 1, 2, 3, 4, and [1 -4]). . . . . . . . . . . . . . . . . . . Figure 46 Comparison of average ROC accuracies for segmentation results with region growing on SV-enhanced and Frangi's and Sato's vesselness markers, anisotropic diffusion, and original image. Comparison is performed on the original phantom image and the phantom corrupted with four levels of Gaussian noise: 10, 20, 40, and 80. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 47 Region growing segmentation results on (a)original phantom data with Frangi's and (b) Sato's vesselness (right) seeds (orange) mapped with the ground truth (white). . . . . . . . . . . . . . . . Figure 48 Comparison of average ROC accuracies of region growing segmentation results based on four methods according to scales σ = {1, 2, 3, 4, [1 -4]}. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 49 Region growing segmentation results (orange) superimposed with the ground-truth image (white).(a): ground-truth. (b): morpho-Hessian enhanced. (c): anisotropic diffusion enhanced. . . . . . . . Figure 50 Comparison of average ROC accuracies of skeletons results based on four methods according to synthetic images. . . . . . . . . . . . Figure 51 Skeletons of the image with noise k = 80 at σ = 2. (a): Simple skeleton. (b): Prunned skeleton. (c): The two superimposed: white -(a), orange -(b). . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: 3D model of major vessels and bifurcations of the human arterial tree reconstructed from a set of CT, DSA CT and MRA images. Colors represent different parts of the model. Left: Aorta and adjacent arteries. Right top: Cranial arterial network. Right bottom: Carotid artery (source: [GAM + 09]).

Figure 2 :

 2 Figure 2: (a) Illustration of the normal arterial-capillary-venous relationship. (b) Drawing of the failure of normal capillary genesis, a fistular connection between the developing artery and vein occurs, which causes an arterial shunting of blood with arterialization of the venous outflow (source: [MHMM09]). (c) Illustration of cerebral AVM.

Figure 3 :

 3 Figure 3: A fluoroscopic C-arm machine that generates x-rays from one side and photographs them on the other side.

  Skeleton and watershed-like approaches (geometrical ridges) 3. Orientation-scale-space [MCL + 07] 4. Edge detection[START_REF] Canny | Finding edges and lines in images[END_REF]

2. 1

 1 thin objects filtering

Figure 4 :2. 1 . 1

 411 Figure 4: A thin, noisy brain blood vessel. There is no low-gradient zone in this object, disconnections are due to noise.

Figure 5 :

 5 Figure 5: Volume rendering of MRA volume, from left to right: original image, ground truth, entropy-based measure from [QBD + 09], Hessian-based measure. (Illustration is from [TdTF + 07]).

Figure 6 :

 6 Figure 6: Geomteric vessel models: (a) Hessian eigenvalues' shape-space according to [LD01]. (b) Superellipsoid model by [TdTF + 07].(Illustrations are from according articles.)

Figure 7 :

 7 Figure 7: Schema of segmentation with level-sets. (a) 2D case (left -evolving curve, right -curve evolution with level sets). (b) Codimension-two curve evolving on tubular isolevel set Γ of C. Right: the tangent to C at p (the normal plane), the external vector → d and its projection onto the normal plane.(Illustration is from [LFG + 99]).

Figure 8 :

 8 Figure 8: MRA segmentation results using: CURVES (left), capillary active contours (right).(Illustration is from [YK05]).

Figure 9 :

 9 Figure 9: Tracking scheme from [BRT03]. (Illustration is from the corresponding publication.

Figure 10 :

 10 Figure 10: A sequence of superellipsoidal fits traversing vessel. (Illustration is from [TdTF + 07]).

Figure 11 :

 11 Figure 11: Brain vessels segmentation based on grey-level hit-or-miss transform. Illustration is from [NPR07].

Figure 12 :

 12 Figure 12: Results of segmentation of vessels by[CA04]. (a) Reconstructed images for the abdominal arteriography images. (b) Labeled result obtained by directly applying the level set method. (c) MIP of enhanced vessel tree. (d) Interim and final results of the evolution of the propagating surface. (e) The final corrected vessel surface. (f) An enlarged view of part of the final vessel surface triangulated with 3-D Delaunay triangulation. (g) The intersection contours of sagittal, coronal, and axial image planes with the level-set segmented result (white) and the final result obtained after the application of a geometric deformable model segmentation technique (gray).(Illustration is from[START_REF] Chen | Quantifying 3-D vascular structures in MRA images using hybrid pde and geometric deformable models[END_REF]).

Figure 13 :

 13 Figure 13: ROC curve for Frangi method and algorithm of [GAF + 09] (source: [GAF + 09]).

Figure 14 : 1 Type

 141 Figure 14: 3D neighborhoods: 6-connectivity, 18-connectivity and 26-connectivity (from left to right).

Figure 15 :

 15 Figure 15: Illustration of the concept of firefront. Fire propagates from the border of the object toward its interior, in an isotropic manner. The point C is a quench point and is the center of a maximal ball (source: [NT10a])

Figure 17 :

 17 Figure 17: Illustration of medial axes (source: [Blu67]).

Figure 18 :

 18 Figure 18: Discrete medial axis: (a) the silhouette of a bear; (b) its discrete medial axis, with the 4-connected neighborhood as unit ball; and (c) the ultimate erosion with the same unit ball

Figure 19 :

 19 Figure19: The thinning process on an example 2D shape. Boundary points are marked "B" at the beginning of each iteration and then removed if they are simple (source:[START_REF] Nd Cornea | Curve-skeleton properties, applications, and algorithms[END_REF]).

Figure 20 :

 20 Figure 20: Example of a distance transform (right) and its original image (left).

Algorithm 1 : 1 A3

 11 Basic Thinning(X, W, D, k) Data: A k-connected shape X, a priority function D and a subset W ⊆ X Result: A skeleton of X while there exists a k-simple point in X\W do = {y ∈ X\W|y is k-simple for X} 2 B = {x ∈ A| for all y ∈ A, D(x) D(y)} Let

Figure 21 :

 21 Figure 21: A graph based on the skeleton.

Figure 22 :

 22 Figure 22: Extraction of junctions and segments from binary image data. The original vascular structure (left) is thinned to its centerline (center), voxels with one or more than 2 neighbors are marked in red as junction points. Segments are identified by following the centerline voxels from one junction to the next one (right).

Fig. 3 .Figure 23 :

 323 Fig. 3. An example vascular tree using Strahler Ordering, where ( †) indicate examples of trifurcation ordering rules, and ( ‡) indicate junctions upstream of vessel pruning. AJP-Heart Circ Physiol • VOL 291 • JULY 2006 • www.ajpheart.org

Figure 24 :

 24 Figure 24: Hybrid combination of the direct volume rendering of bone structures and surface rendering of intra-hepatic structures (liver, vascular systems, tumours in the right view). The bone structures serve as anatomical reference for the intra-hepatic structures. (Source: [PP03]).

  ... not to reproduce what we can already see, but to make visible what we cannot.

Figure 25 :

 25 Figure 25: Left: the second-order derivative of a Gaussian kernel probes inside-outside contrast of the rage (-σ, σ) (in this example σ = 1). Right: the second-order describes the local principal directions of curvature. (source: [FNVV98]).

Figure 26 :

 26 Figure 26: a) Original image of a tubular object. b) Second derivative detected object and its pixel directions. c) Frangi's vesselness function thresholded for a surface rendering.

Figure 27 :

 27 Figure 27: a) Original image. b) Second derivative scale-space ridges (center). c) Ridges depicted by their boundaries (at scale of ridge). Ridges were computed for γ = 1.25 and scales range from 1.5 to 16 in steps of 0.5. (Illustration is from [Maj04].)

Figure 28 :

 28 Figure 28: Image of disconnected vessel branches (a), result of a 3 × 1 vertical spatially-invariant closing result (b), and a 3 × 1 spatially-variant closing (c).

Figure 29 :Figure 30 :

 2930 Figure 29: The Bayer colour filter array. In modern digital cameras, only one of the colour component is captured at every pixel location.

Figure 31 :

 31 Figure 31: Colour interpolation using a SVMM closing. (a) Original image data ; (b) Bayer sampling simulating a digital camera sensor ; (c) green channel interpolation using an SVMM closing ; (d) full colour interpolation using and SVMM closing.

Figure 32 :

 32 Figure 32: (a) Standard single-neighbour dilation. (b) Standard adjunct erosion by considering the transpose of the structuring element.

Algorithm 2 :Figure 34 :

 234 Figure 34: Result of the spatially variant closing on the tubular object with the directions provided by multi-scale second order derivatives.

Figure 35 :

 35 Figure 35: Vector field regularzation.

Figure 36 :

 36 Figure 36: Regularized vector field of the tubular object and the resulting spatially-variant closing.

Figure 37 :

 37 Figure 37: Slice of a MRA data of the head (left) and a close-up of the AVM.

Figure 17 :Figure 18 :Figure 38 :

 171838 Figure 17: An object in G 3 with, in red, (a) the minimal symmetric K -skeleton produced by MK 3 , (b) a minimal asymmetric K -skeleton produced by AMK 3 .

Figure 39 :

 39 Figure 39: An object and its asymmetric 1D isthmus skeleton (left) and its graph based on the skeleton (right) .

1 .

 1 derivative-based directions and tubular-objects detection ; 2. field regularization ; 3. spatially-variant closing ; 4. seeded region growing segmentation ; 5. hole-closing ; 6. skeletonization with asymmetric 1D isthmus curvilinear method ; 7. 3D skeleton points classification ; 8. Vascular network graph construction ; and 9. Skeleton pruning.

Figure 40 :

 40 Figure 40: Rendering of the synthetic vessel-like object with the added Gaussian noise of κ = 10.

Figure 41 :

 41 Figure 41: Slices of phantom images with four levels of Gaussian noise (from left to right): 10, 20, 40, and 80.

Figure 42 :

 42 Figure 42: Comparison of average ROC accuracies for Frangi's or Sato's vesselness and their respective hysterysis thresholding for original phantom image and with four levels of Gaussian noise: 10, 20, 40, and 80.

Figure 43 :

 43 Figure 43: From left to right: Morpho Hessian with Frangi's normalized and non-normalized vesselness or Sato's measure on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4].

Figure 44 :

 44 Figure 44: From left to right: Frangi's and Sato's vesselness on the phantom image with k = 10. 1st row: σ = 1. 2nd row: σ = 2. 3rd row: σ = 3. 4th row: σ = 4. 5th row: σ = [1 -4].

Figure 45 :

 45 Figure 45: Comparison of average ROC accuracies for Frangi's or Sato's vesselness according to results at four discrete scales and one multi-scale (σ: 1, 2, 3, 4, and [1 -4]).

Figure 46 :

 46 Figure 46: Comparison of average ROC accuracies for segmentation results with region growing on SV-enhanced and Frangi's and Sato's vesselness markers, anisotropic diffusion, and original image. Comparison is performed on the original phantom image and the phantom corrupted with four levels of Gaussian noise: 10, 20, 40, and 80.

Figure 47 :

 47 Figure 47: Region growing segmentation results on (a)original phantom data with Frangi's and (b) Sato's vesselness (right) seeds (orange) mapped with the ground truth (white).

Figure 48 :

 48 Figure 48: Comparison of average ROC accuracies of region growing segmentation results based on four methods according to scales σ = {1, 2, 3, 4, [1 -4]}.

Figure 49 :

 49 Figure 49: Region growing segmentation results (orange) superimposed with the ground-truth image (white). (a): ground-truth. (b): morpho-Hessian enhanced. (c): anisotropic diffusion enhanced.

Figure 50 :Figure 51 :

 5051 Figure 50: Comparison of average ROC accuracies of skeletons results based on four methods according to synthetic images.

Figure 52 :

 52 Figure 52: Comparison of average ROC accuracies of skeletons results based on four methods according to corresponding scales.

Figure 53 :

 53 Figure 53: 1st line: Phatom skeleton -ground-truth. 2nd line: skeletons of morpho-Hessian with Frangi's vesselness results. From left to right images with noise level κ = 10, 20, 40, 80 on scales σ = {1, 2, 3, [1 -4]}, respectively. 3rd line: skeletons of anisotropic diffusion results, for the same images and at the same scales as on the previous line.

Figure 54 :

 54 Figure 54: (a) Sample of an original MRA showing small disconnected vessels at the resolution limit. (b) Surface rendering of the smallest disconnected vessels. (c) Surface rendering of the original image (higher threshold than in (b)) and Hessian-based directions. (d) Regularized vector field. (e) Geodesic reconstruction of the thresholded Frangi's vesselness into the morpho-Hessian filtered image. (f) Anisotropic diffusion result as in [MVN06]. (g) Region growing result and its skeleton. (h) Surface rendering of the original image and the skeleton of region growing.

Figure 55 :

 55 Figure 55: Frangi's vesselness (yellow) and Sato's with the same thresholds.

Figure 56 :

 56 Figure 56: Whole head volume segmentation (on the left). Its skeleton with the AVM (yellow).

Figure 57 :

 57 Figure 57: (a) MIP and the selected AVM. (b) AVM with its skeleton. (c) AVM with more vessel connections and (d) their skeleltons.

Figure 59 :

 59 Figure 59: Example of a guidewire image (cropped and contrast-enhanced).

Figure 61 :

 61 Figure 61: Illustration of noise sensitivity of the local orientation operator. The top image is the original image. The bottom left image is obtained using the local orientation operator with morphological openings. The bottom right image is obtained using a parametric opening with a rank parameter of 20%, which makes it less noise-sensitive. (Source: [Soi98].)

Figure 3 .

 3 Figure3. Electron micrograph of glass fibres: to detect the small thin fibres in the bottom of the image, a white top-hat is useful but noisy. When this top-hat image is filtered by an area opening some compact noise remain while a path opening yields a better result.

Figure 62 :

 62 Figure 62: Electron micrograph of glass fibres: to detect the small thin fibres in the bottom of the image, a white top-hat is useful but noisy. When this top-hat image is filtered by an area opening some compact noise remains while a path opening gives a better result. (Source: [AT05].)

Figure 63 :

 63 Figure63: Ground truth (GT): dark green -provided GT, light green -dilated (diameter of the balls structuring element is 7 pixels) GT used as a true positive here; blue -dilated GT (size 20 pixels) used as a negative. The rest of the image is ignored.

6.2 guidewire detection results 115 Figure 64 :

 11564 Figure 64: Graph of results for three guidewire images.

Figure 65 :

 65 Figure 65: Graph of the best results obtained with rank-opening for three guidewire images.(Points in red are the best results.)

Figure 66 :

 66 Figure 66: ROC curves of accuracies of hysteresis thresholded results of a) Frangi's vesselness , b) Frangi's vesselness on morpho-Hessian filtered image. (Points in red are the best results.)

Figure 67 :

 67 Figure 67: 1st line: the original guide-wire images (cropped and contrast-enhanced) and their corresponding best segmentation results. 2nd line: rank-opening, 3rd line: path opening, 4th line: Frangi's vesselness with hysteresis thresholding, 5th line: structure tensor, 6th line: morpho-Hessian filter, and 7th line: steerable filter.

Figure 67 :

 67 Figure 67: 1st line: the original guide-wire images (cropped and contrast-enhanced) and their corresponding best segmentation results. 2nd line: rank-opening, 3rd line: path opening, 4th line: Frangi's vesselness with hysteresis thresholding, 5th line: structure tensor, 6th line: morpho-Hessian filter, and 7th line: steerable filter.

Figure 68 :

 68 Figure 68: Model-based vessel visualizations with different settings for the smoothing of the vessel. The vessel diameter is color-coded. In the left image, no smoothing is applied which results in obvious discontinuities in the periphery, which are due to the limited spatial resolution of the underlying data (source: [PP03]).

Table 1

 1 Graph of results for three guidewire images. . . . . . . . . . . . . . Figure65Graph of the best results obtained with rank-opening for three guidewire images. (Points in red are the best results.) . . . . . . . Topological classification of 3D points according to the values of C. 41

	Figure 63	Ground truth (Ground Truth (GT)): dark green -provided GT,
		light green -dilated (diameter of the balls structuring element is 7
		pixels) GT used as a true positive here; blue -dilated GT (size 20
		pixels) used as a negative. The rest of the image is ignored. . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 61 Illustration of noise sensitivity of the local orientation operator. The top image is the original image. The bottom left image is obtained using the local orientation operator with morphological openings. The bottom right image is obtained using a parametric opening with a rank parameter of 20%, which makes it less noisesensitive. (Source: [Soi98].) . . . . . . . . . . . . . . . . . . . . . . . Figure 62 Electron micrograph of glass fibres: to detect the small thin fibres in the bottom of the image, a white top-hat is useful but noisy. When this top-hat image is filtered by an area opening some compact noise remains while a path opening gives a better result. (Source: [AT05].) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . line: rank-opening, 3rd line: path opening, 4th line: Frangi's vesselness with hysteresis thresholding, 5th line: structure tensor, 6th line: morpho-Hessian filter, and 7th line: steerable filter. . . . . Figure 68 Model-based vessel visualizations with different settings for the smoothing of the vessel. The vessel diameter is color-coded. In the left image, no smoothing is applied which results in obvious discontinuities in the periphery, which are due to the limited spatial resolution of the underlying data (source: [PP03]). . . . . . L I S T O F TA B L E S

Table 2

 2 Possible eigenvalue responses and their signs after the eigenanalysis corresponding to different shapes and color intensities. . . . 61

xvii 1

Table 1 :

 1 Topological classification of 3D points according to the values of C.

  )

			4.1 feature detection with gaussian derivatives	61
	2D			3D		Object's profile
	λ 1	λ 2	λ 1	λ 2	λ 3
	Big+	Big+	Big+	Big+ Big+	blob-like (dark)
	Big-	Big-	Big-	Big-	Big-blob-like (bright)
	Small Big+ Small Big+ Big+	tubular (dark)
	Small Big-	Small Big-	Big-	tubular (bright)
			Small Small Big+	surface (dark)
			Small Small Big-	surface (bright)

Table 2 :

 2 Possible eigenvalue responses and their signs after the eigenanalysis corresponding to different shapes and color intensities.

(a) (b) (c) Figure 58: (a) MIP and the selected AVM. (b) AVM with its skeleton. (c) AVM with more vessel connections and (d) their skeleltons.

Plots in this figure have more points than in Fig.65due to the fact that the thresholding range for vesselness is bigger and this is a hysteresis thresholding.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Part III A P P E N D I X

D E V E L O P M E N T A N D T E S T E N V I R O N M E N T

During this thesis work there have been used many university libraries, public and commercial tools.

The derivative-based and morphological parts of the morpho-Hessian filter were written with the ITK library based on the available implementations of filters of Frangi's and Sato's multi-scale filters. Also, the code of the anisotropic diffusion filter [START_REF] Enquobahrie | Vessel enhancing diffusion filter[END_REF] was used within this library.

The regularization of the direction field, the region growing, the structure tensor were written by Hugues Talbot within the library Voir.

The skeletonization procedures and numerous other morphological and arithmetic image analysis functions were performed with the library Pink (author -Michel Couprie). Other skeletonization procedures based on cubical complexes were used with the library of John Chaussard.

The test environment is written with the Bash code with the generous and thorough help of Hugues Talbot.

The results of steerable filters were obtained by Vincent Bismuth.

The automatic performance evaluation of the method is written with R language also with help of Hugues Talbot.

The specialised visualization modules for Avizo, Amira were kindly provided by Erwan and Nicolas Combaret.

Other free and commercial programs used for the development and visualization programs are: P U B L I C AT I O N S Some ideas and figures have appeared previously in the following publications: