posteriori probability (MAP) solution.

Recall that the posterior probability can be stated, according to the Bayes rule, in terms of prior, likelihood and observation probabilities as

Summary

We develop generic methods for the efficient semiautomatic multiobject segmentation of high-resolution D medical images. We formulate segmentation as a Bayesian labeling problem defined in terms of a prior distribution, image-derived evidence and user constraints and solve it by minimizing an energy with multilabel Graph Cut algorithms which produce provably good approximate solutions and in some cases global optima.

The most widely used prior in Graph Cut segmentation is the Potts model. It implies the presence of a single boundary in the image, which is untrue in multiobject segmentation. We investigate extensions of this model to account for multiple boundaries in multiobject segmentation which allow the minimization of resulting energies by Graph Cut algorithms. Our main contribution is a structural prior model, called the vicinity prior, which captures the spatial configuration of objects. It is formulated as shortest-path pairwise constraints on a graph model of interobject adjacency relations.

Graph Cut algorithms use grid graphs which become complex on large input making optimization impractical. We reduce memory budget and optimization runtime by clustering image voxels prior to segmentation by an adaptive centroidal Voronoi tessellation which strikes a good balance between the conflicting goals of cluster compactness and object boundary adherence and help to place segmentation boundaries precisely. We formulate segmentation according to the graph of clusters such that its energy is independent of clustering resolution. This constitutes our second contribution.

Qualitative and quantitative evaluation and comparison with the Potts model on synthetic, simulated and real medical images confirm that the vicinity prior model achieves significant improvements in the correct segmentation of distinct objects having identical intensities, the accurate placement of structure boundaries and the robustness of segmentation with respect to clustering resolution. Comparative evaluation of our clustering method with competing ones confirms its benefits in terms of runtime and quality of produced partitions. Compared to voxel segmentation, the clustering step improves both overall runtime and memory footprint of the segmentation process up to an order of magnitude virtually without compromising the segmentation quality.

Structural prior models can be derived from other relationships as well, such as distance, area of common surface and spatial directionality. Graph models describing such relations are fairly easy to define at an arbitrary level of detail from simple specifications of relationships on pairs of objects, and can be generated from existing anatomical models in medical applications. Furthermore, relations can be learned from training data or several variants of a model thus capturing its variabilities. i

Présentation des travaux

La pratique clinique et la recherche médicale génèrent aujourd'hui un grand nombre d'images D de haute résolution, spécifiques aux patients, rendant difficile les tâches d'accès à ces données, de transfert, d'analyse et de visualisation, particulièrement dans les environnements informatiques distribués qui témoignent d'une utilisation des terminaux mobiles pour la visualisation interactive des données anatomiques.

Traditionnellement adaptées à un seul organe ou une pathologie particulière, les méthodes de traitement d'images récentes ont rendu possible l'analyse et la segmentation simultanées de plusieurs structures anatomiques. Par ailleurs, les technologies de traitement et de visualisation à distance, pour les grands volumes de données, ont mûri et ouvert la porte à la navigation sémantique de l'anatomie du corps entier. Ainsi, la création de modèles anatomiques et physiologiques «patient-spécifiques» devient crédible.

La segmentation d'objets (d'organes) dans les images est une condition préalable à leur analyse, manipulation et visualisation efficaces. Dans les applications médicales elle est l'épine dorsale de nombreuses procédures cliniques allant du diagnostic aux soins post-opératoires. Elle est également essentielle à d'autres domaines de la médecine impliquant l'étude et la représentation de l'anatomie. Cependant, l'analyse et la segmentation automatiques des images médicales est difficile pour plusieurs raisons. Différents artefacts, liés à la modalité d'imagerie et aux conditions d'acquisition, dégradent la qualité des images et réduisent leur contenu informationnel. En outre, de nombreuses structures distinctes manifestent des niveaux d'intensité (niveau de gris dans l'image) identiques et se fondent dans les tissus voisins. Ces problèmes, couplés avec la variabilité des tissus à travers les populations de patients, rendent irréalisable la conception d'algorithmes de segmentation d'usage générique. Alors que certaines tâches élémentaires, dans des conditions bien définies, sont déléguées aux ordinateurs pour un traitement automatique, les experts médicaux ont fréquemment recours à des méthodes de segmentation essentiellement manuelles dans la construction des modèles D «patient-spécifiques» qui demandent des investissements majeurs de temps et d'effort.

Afin d'alléger le fardeau de la délimitation manuelle et comme une alternative aux méthodes automatiques, souvent irréalistes, les méthodes de segmentation semiautomatiques cherchent un équilibre entre les méthodes manuelles et automatiques. Elles exploitent les facultés d'un opérateur humain à reconnaître des objets de façon iii plus fiable qu'un ordinateur, et les capacités d'un ordinateur à les délimiter de façon plus efficace qu'un opérateur humain. De nombreuses méthodes reposant sur différents paradigmes d'interaction ont été conçues pour permettre à un utilisateur d'initialiser, de piloter un processus de segmentation et d'affiner son résultat avec une interaction nécessitant un effort raisonnable. L'intervention de l'utilisateur fournit aussi des indices, pour obtenir des informations a priori, qu'un algorithme de segmentation peut ensuite exploiter, tels que les statistiques d'intensité, l'emplacement et la disposition spatiale des objets ciblés.

Nous développons une méthode générique pour la segmentation multi-objet semiautomatique efficace d'images médicales D de haute résolution qui pourrait être facilement adaptée à un problème particulier par une simple reformulation du modèle d'a priori que nous proposons et une appropriation de l'algorithme aux besoins de l'application en termes de temps d'exécution et de capacité de stockage par une simplification adaptative de l'image à segmenter.

Un nombre d'approches de segmentation multi-objets variationnelles et combinatoires ont eu succès dans la dernière décennie. Notre approche se base sur les méthodes de coupure de graphe qui s'inscrivent dans le dernier groupe. Ces méthodes sont connues pour leurs formulations probabilistes discrètes flexibles qui se généralisent naturellement aux problèmes multi-objet et pour leurs mécanismes d'optimisation numériquement stables qui, pour des fonctionnelles d'énergie à minimiser sous-modulaires, produisent des solutions locales de bonne qualité, à une distance connue de l'optimum global, et, dans certains cas, la solution correspondant à l'optimum global. Nous formulons la segmentation comme un problème d'étiquetage bayésien, défini en termes de distribution a priori, décrit par un champ aléatoire de Markov, de vraisemblance dérivée de l'image et des contraintes d'utilisateur. Nous le résolvons en minimisant une énergie à l'aide des algorithmes d'optimisation multi-étiquette par coupure de graphe.

De nombreuses formes d'a priori sur l'emplacement, la forme et la disposition spatiale des objets ont été introduites dans les approches de segmentation par coupure de graphe pour différentes applications afin d'assurer la cohérence spatiale de l'étiquetage par rapport à un modèle de référence. Le modèle d'a priori le plus souvent utilisé est celui de Potts dont le mécanisme de pénalisation à deux niveaux implique la présence d'une seule frontière dans l'image; une hypothèse fausse en segmentation multi-objet. Nous étudions des extensions de ce modèle qui tiennent compte des multiples frontières figurant en segmentation multi-objet et qui, en outre, permettent la minimisation des énergies résultantes à l'aide des algorithmes d'optimisation par coupure de graphe. Notre contribution principale est un a priori de structure, appelé «modèle de proximité». Sous-modulaire par définition, il introduit plusieurs niveaux de pénalisation et représente la configuration spatiale des objets en segmentation multi-objet. Il est formulé à partir des contraintes de paires du plus court chemin définies sur un modèle de graphe décrivant les relations d'adjacence entre les objets.

Les algorithmes d'optimisation par coupure de graphe utilisent des structures de données sous forme de grille qui deviennent complexes sur des grands volumes d'images D rendant l'optimisation impraticable. Nous observons que l'image en pleine résoiv lution est souvent redondante et réduisons l'empreinte mémoire et le temps d'exécution de l'algorithme en partitionnant l'image avant de la segmenter par une tesselation de Voronoï barycentrique adaptative. Les partitions de l'image ainsi créées établissent un bon équilibre entre les objectifs concurrents de la compacité des régions et du respect des frontières des objets, permettant en outre de contrôler le compromis, et contribuent à la mise en place fine des frontières de la segmentation ultérieure. Nous formulons la segmentation par coupure de graphe selon le graphe des régions de Voronoï tel que son énergie est indépendante de la résolution des partitions. Ceci constitue notre deuxième contribution.

Des évaluations qualitatives et quantitatives, sur des images médicales synthétiques, simulées et réelles de tailles et de contenus différents, confirment les avantages de notre modèle d'a priori de proximité sur celui de Potts en segmentation multi-objet par coupure de graphe comptant jusqu'à quelques dizaines d'objets. Les améliorations apportées par notre modèle sont significatives, en particulier dans la segmentation correcte des objets distincts ayant des intensités similaires, dans le positionnement précis des frontières des structures ainsi que dans la robustesse de segmentation par rapport à la résolution de la partition.

L'évaluation comparative de notre méthode de partition avec ses concurrentes confirme ses avantages en termes de temps d'exécution et de qualité des partitions produites. Par comparaison avec l'approche de segmentation appliquée directement sur les voxels des images, l'étape de partition améliore à la fois le temps d'exécution global et l'empreinte mémoire du processus de segmentation jusqu'à un ordre de grandeur sans compromettre la qualité de la segmentation en pratique.

Les travaux futurs portent sur des objectifs à court et à long terme. Des variantes parallélisées de l'algorithme d'optimisation sous-jacent peuvent directement être utilisées afin d'accélérer encore le processus de segmentation défini sur des partitions d'images. L'implantation expérimentale mise en oeuvre afin de développer et valider la méthode peut être élaborée dans une application interactive intégrant le geste de l'utilisateur et l'algorithme de segmentation en temps réel. À plus long terme, de nouveaux modèles d'a priori de structure peuvent être définis à partir d'autres relations, notamment la distance, l'aire de la frontière commune et l'orientation spatiale. Les modèles de graphe décrivant ces relations sont faciles à définir à un niveau arbitraire de détail à partir de simples spécifications des relations sur des paires d'objets. De manière plus intéressante, ces graphes peuvent être générés automatiquement, soit à partir des modèles anatomiques existants pour les applications médicales, ou par des méthodes d'apprentissage en déduisant les relations sous-jacentes à partir des données étiquetées ainsi représentant ses variabilités possibles. v First and foremost, I would like to thank my thesis advisors Messrs. Rémy Prost and Michel Desvignes for their guidance and trust throughout the past three years and half, and especially my supervisor and close collaborator Mr Sébastien Valette without the assistance, the meticulous attention to detail and the encouragement of whom this work would have never come to fruition. The insight and the invaluable feedback these gentlemen generously provided me with over the course of elaborating the present treatise undoubtedly shaped its scientific content and vigor, keeping it relevant to the times and current research trends. It should not go without notice that these three men were also instrumental in helping me meet and fulfill many deadlines and ridiculous administrative chores the life of a PhD candidate is never free of, often borrowing time from weekends and holiday leaves.

I extend my thanks and gratitude to the funding organization La Région Rhône-Alpes (France) and particularly to Mr Desvignes for generously and patiently allowing me the time to complete the thesis manuscript after recently recruiting me into his research group. I am looking forward to new research endeavors with him.

Finally, I address very special words of gratitude to my better half, Lilie. It is unavoidably cliché to thank beloved ones and family members; their support and inspiration can never be fully appreciated and returned. While the writer agonizes over his oeuvre, they suffer through his absence. Creation is beautiful agony, the wait is not. By nature, writing is a selfish process, and the writer only hopes to achieve salvation by the consciousness that his work would offer some inspiration and insight to readers. vii

M

edicine owes much of its developments to advances in anatomy throughout the history of modern science. Although the study of anatomy begins at least as early as  BC in ancient Egypt, the treatise of Vesalius [Ves] was the first published work on anatomy to overthrow Galenic and Arabian medical doctrines which prevailed until the Renaissance. Prioritizing dissection, Vesalius established a scientific method introducing what is today called an anatomical view of the body which sees the human internal essentially as a corporeal structure filled with organs arranged in a three-dimensional space. Vesalius would not have been able to accomplish his seminal work without the advances in art during the Renaissance and the technical developments of printing press at the same time. The latter technology in particular allowed him to produce highly detailed illustrations of human anatomy accurately depicting his findings. In the th century, the finalization and the systematization of descriptive human anatomy developed during the previous centuries, along with growing contributions from new sources of knowledge such as histology and developmental biology, culminated in Gray's Anatomy [Gra], an influential anatomical atlas which continues to be revised and republished to the present day in print and digital media [Sta]. It is considered one of the most comprehensive and detailed references on the subject of human anatomy.

With the revolution of medical sciences in the th century, anatomical research has taken advantage of technological developments and expanding knowledge in new disciplines, such as endocrinology, evolutionary and molecular biology, to achieve a thorougher understanding of the human body's organs and structures. Importantly, medical imaging devices, such as magnetic resonance imaging (MRI) and computed tomography (CT) scanners, have enabled researchers to inspect the anatomy of living as well as of dead organisms after being limited to dissection studies on cadavers for centuries. Advances in medical imaging technologies simultaneously with the development of computing platforms and computer graphics technologies together with standardization efforts, such as DICOM [Pia], moved many medical disciplines to  the digital realm, while anatomical research and education increasingly relied on digital photography and three-dimensional computer-generated representations. A major milestone at the turn of the century, perhaps of equal significance to Vesalius' contribution for its time, was the release of the Visible Human dataset by the U.S. National Library of Medicine [SASW], which is an open-access digital image library of the adult male and female anatomy comprising high-resolution cryosectional color photographs of unprecedented detail, as well as the corresponding MR and CT images for the same subjects. This dataset has formed the basis of many biomedical and engineering applications and medical education platforms  which have all exploited computer and visualization technologies for the examination of anatomy in its natural three-dimensional setting. In particular, three-dimensional representations reconstructed from images in the Visible Human dataset have made it possible to view human anatomy from perspectives unavailable in dissections and thus have helped to identify and correct several inaccuracies in two-dimensional illustrations of conventional anatomical atlases [VIM]. Today, as the macroscopic aspects of human anatomy have been largely cataloged, the research focus is on the development, the evolution, and the function of anatomical features, particularly on microscopic scales.

 For recent bibliography, refer to the Visible Human project web site [U.S]  . From the archetype to the subject-specific

The importance of anatomical understanding in modern medical practice cannot be stressed enough. A thorough working knowledge of gross anatomy, and sometimes of microscopic anatomy as well, is required for most medical workers; physicians, dentists, physiotherapists, nurses, paramedics, radiographers and especially surgeons, not only to support diagnosis but also in preventive care, therapy and particularly in preoperative planning, surgical intervention and postoperative care. Such clinical procedures usually involve an imaging step by the appropriate modality, depending on the tissue or structure type, which often constitutes the only means to discover a disease and proceed with diagnosis, particularly in the case of cancerous tumors. Obviously the accuracy of diagnosis and the effectiveness of surgical intervention depend crucially on the quality of the acquired image, the sufficiency of its coverage of the pathology and, more importantly, the relevance and the precision of information that a practitioner can extract from the image. Although, in general, a trained practitioner is considered of superior image analysis and pathology identification abilities, relegating some limited tasks to computers is of great practical interest; for example, the location and the isolation of a tumour from affected structures in the image for measurement and the outlining of moving structures for tracking purposes are routinely performed tasks on present-day imaging devices. Such tasks are carried out by "learned" software programs entailing, in almost all cases, structure recognition and delineation mechanisms, that is, image segmentation algorithms.

In addition to their critical importance in clinical practice, medical imaging technologies are exploited by other disciplines of medical science as well, notably in the study of anatomy and physiology. In many applications in these fields, medical images of different modalities are frequently analyzed for structure segmentation in preparing anatomical and physiological models for the diverse purposes of research, medical education and visualization. A noteworthy example is the creation of anthropomorphic phantoms for applications in health physics and nuclear medicine, such as whole-body and individual-organ radiation dosimetry for radiation protection and radiotherapy. Computational procedures of such studies usually involve Monte Carlo simulations to track radiation interactions and energy deposition in the body. Needless to say, the accuracy of such simulations largely depend on the precision and the realism of underlying anatomical representations. The construction of whole-body anatomical models from real images naturally entails a nontrivial task of identification and segmentation of, for the applications just mentioned, at least several tens of structures from a series of whole-body images, in some cases originating from different modalities to account for different anatomical structures. Unsurprisingly, the Visible Human image dataset has formed the basis of several such digital phantoms, notably the voxel-based VIP-Man [XCB] and Voxel-Man [PHP + ] phantoms and the NURBS-based Cardiac-Torso (NCAT) phantom [Seg]. All of these phantoms were constructed from the cryosectional color photographs in the dataset by trained anatomists mostly by a painstaking process of manual contouring and image volume editing complemented by elemen-tary computer-assisted techniques such as thresholding and color-space classification. However, unlike the dataset from which they were created, none of these phantoms are available to the scientific community on an open-access basis.

While in the past cadavers were the definitive source of anatomical information, it seems ironic today that with the availability of open-access digital anatomical datasets, such as the Visible Human image library, and the effective means to process and visualize them, that is, powerful computing platforms and computer graphics technologies, open-access three-dimensional anatomical atlases and models are lacking. The sheer volume and anatomical coverage of image data to process, coupled with difficulties associated with the analysis and the segmentation of such images by automatic unattended methods prompts researchers and engineers to resort to manual methods in the construction of three-dimensional whole-body models which require significant investments in time and labor making it often infeasible to give away its outcome.

Clinical practice today, especially whole-body CT and MRI scanning, often generates large numbers of high-resolution images, which makes tasks of efficient data access, transfer, analysis and visualization challenging, especially in today's distributed computing environments which have seen growing use of handheld terminals for interactive data access and visualization of anatomy. With increasing availability and sophistication of industrial-grade medical imaging technologies, unprecedented amounts of patient-specific anatomical data is now at the disposal of medical practitioners. Traditionally single-object or pathology oriented, recent image processing methods, discussed later in this chapter, have made the simultaneous analysis and segmentation of multiple anatomical structures increasingly possible. Moreover, remote processing and visualization technologies for large datasets, such as [JJAM], have matured and opened the door for a number of popular applications, notably the semantic navigation of full-body human anatomy [BCK + ]. Thus, the creation of patient-specific anatomical models, at least systemic or disease-specific, is becoming of great interest.

. Why is it difficult to segment medical images?

Despite more than  years of progress, image segmentation remains one of the major challenges and among most active research topics in image processing and computer vision [FP]. As we have seen, segmentation of objects in images is a prerequisite to their effective analysis, manipulation and visualization. In medical applications in particular, segmentation of anatomical structures in medical images is the back-bone of many clinical procedures from diagnosis to surgical planning and intervention, and is crucial to other fields of medical science involving the study and the representation of human anatomy. The segmentation of high-resolution medical images for the purpose of creating application-specific reduced representations of targeted structures, such as surface and volume meshes, is also interesting from the viewpoint of efficient data storage, dissemination and three-dimensional visualization. CT and MR images have intrinsic characteristics that render its analysis and segmentation, especially in automatic unattended mode, particularly challenging [BSJB]. In theory, they are regarded as piecewise-constant intensity maps over a number of tissue classes. In practice, however, the piecewise-constance property is degraded considerably by a number of factors. Even with significant technological advances, current image acquisition systems often yield low signal-to-noise ratio data, consequently produced images are commonly degraded by various noise sources. The limited resolution of imaging systems often results in the dual effects of partial volume and spillover leading to structure volume under-and overestimation respectively. These, in addition to other image artifacts due, for example, to patient motion or limited acquisition time, reduce the prominence of intensity edges in images creating diffuse or poorly defined ones which make it difficult to identify true structure boundaries. MR images may, in addition, suffer from spatial distortion of structure intensity due to magnetic field inhomogeneity and gradient field nonlinearity. Good contrast between tissue class intensities, which greatly facilitates segmentation, is usually achieved by using contrast agents depending on the imaging modality, in the lack of which contrast enhancement postprocessing methods are resorted to, which often give less satisfactory results. However, regardless of the imaging modality and related artifacts, many anatomically and functionally distinct structures, especially those corresponding to soft tissues, have similar intensity levels in images and, furthermore, blend into surrounding tissues which have intensities close to their own. It is impossible for an automatic process to identify and segment such structures on the basis of intensity information only, hence most advanced segmentation methods exploit some form of prior information on structure location, spatial extent or interrelations to achieve greater robustness and precision.





. The choice of a segmentation approach

Segmentation, as the reader can tell, consists of the two closely-interrelated tasks of structure recognition and delineation. Recognition involves identifying the location of the targeted object and distinguishing it from others that are present in the image, and delineation involves designating its precise spatial extent. According to the degree of automation in these two tasks, existing segmentation methods cover a wide spectrum of approaches from completely manual tracing or painting of object boundaries or regions to fully automatic recognition and delineation of objects. Automatic methods have had considerable success in many well-defined application contexts, usually coinciding with fixed image acquisition protocols and little variation in object pose, shape and location in the image. Due to the availability of image datasets acquired in such conditions by routine clinical or laboratory procedures, most automatic segmentation approaches rely on image and object features learned from presegmented datasets. This makes them prone to be bound to the characteristics of training sets and therefore to have difficulty in accounting for previously unencountered but possible large variations in object pose, location and shape, typically arising in pathological cases, and the photometric features of images acquired by different protocols in different conditions. Semiautomatic segmentation methods are usually resorted to in order to lighten the burden of manual delineation and as an alternative to fully-automatic methods when the latter are inapplicable, for example due to the singularity of the studied case or the lack of training sets. Semiautomatic methods achieve a balance between manual and automatic methods exploiting the synergy between a trained expert's superior object recognition abilities and a computer's faster delineation capacities. Many methods relying on different interaction paradigms have been devised to allow a user to initialize, steer a segmentation process and refine its result with a reasonable amount of input. User interaction and feedback provide clues on what he intends to segment which can be exploited by a segmentation algorithm in a number of ways, such as to constrain segmentation boundary placement or the assignment of a pixel to some object, and to collect intensity statistics of objects and derive prior information on their interrelations.



User interaction is usually carried out through a traditional human-computer interface using a pointing device to mark certain image pixels in order to indicate an object or a boundary type. Often called "seeds", such pixel markings are consequently interpreted by the segmentation algorithm as delineation guides or pixel assignment constraints. Commonly encountered user-specified constraints fall in three categories: object boundary, object region and image region constraints. A popular object-boundary constraint formulation for semiautomatic segmentation was introduced in D by [MB, FUS + ] and generalized to surface constraints in D by [FU]. In these approaches, the user roughly traces the boundary of each object indicating a number of seed points through which each segmentation boundary should pass. The segmentation algorithm is tightly coupled with user gesture and finds the best boundary as the minimum-cost path between the current and the last specified seed points on the image graph according to some image gradient-based local cost. If a computed path deviates from the desired one, due to a weak intensity edge for example, the user specifies additional seeds. The intricacy of this method is obvious; the user is required to place seeds accurately on the object boundary. By contrast, approaches relying on object-region constraints require that the user places seeds anywhere inside targeted objects. It is the segmentation algorithm's duty to ensure that marked pixels are part of corresponding objects in the result. Assuming that objects consist of fairly homogeneous regions, such constraints are quite robust to placement precision and, compared to object-boundary seeds, their specification requires much less time and effort and can be carried out in D in a straightforward fashion by placing seeds in D slices of the image volume corresponding to a chosen cross-sectional view. The interpretation and the use of such seeds often depends on the algorithm. For example, [VK] define seeds as the starting point of their cellular automata-based region growing algorithm. In [GFL], given a number of user-specified seeds with label types corresponding to targeted objects, the algorithm uses an analytic method to determine, for every pixel in the image, the probability that a random walker starting its walk at that pixel first reaches a seed of a given label type. The pixel is then assigned the label with the greatest probability. Combinatorial Graph Cut methods [BF], which describe a segmentation in terms of optimal regional and boundary properties of objects, interpret object-region seeds as hard constraints for pixels that should be part of the corresponding object in the result, also using them to derive intensity statistics of objects in addition to allowing the user to edit and refine an initial segmentation efficiently by corrective seeds. Image region-based constraints entail interaction methods that are even less stringent than object-region based ones and often do not involve the marking of any particular pixel belonging to some object. For example, [RKB] only require that the user draws a bounding box around the targeted object in the image separating it from the background. The algorithm then refines background and foreground statistics iteratively as it converges to the optimal separating boundary. Figure . illustrates the three discussed interaction paradigms. Recent gesture-based multitouch interfaces offer new user interaction possibilities the use of which in simultaneous semiautomatic segmentation and visualization of anatomical data is undoubtedly an interesting research venue in its own right.



A number of successful variational and combinatorial segmentation methods have matured in the last decade, notably level-set based active contour [OF] and maximumflow based Graph Cut [BV] methods. While the first group of approaches formulate segmentation as a continuous optimization of contour or surface functionals and the latter as a discrete optimization of cost functions on the image graph, both groups of methods rely on implicit nonparametric representations of object boundaries. Level-set methods define the implicit representation of a contour or a surface as the zero-level set of a higher dimensional function and describe its evolution towards the desired object boundary by a partial differential equation. The contour movement is driven by external image-based and internal contour-based speed functions, which, respectively evolve the contour towards image features, such as intensity edges, and smooth its boundary. In Graph Cut segmentation, targeted objects are represented by labels and a discrete labeling cost function of variables corresponding to image pixels is formulated in terms of label likelihood and spatial regularization criteria. The cost function is then represented by a flow network the minimum cut of which, calculated by a maximum flow algorithm, finds the optimal segmentation boundary. As research in both fields has progressed, state-of-art Level-set and Graph Cut segmentation methods have come to incorporate different object appearance, boundary and shape based cues and prior information [CRD] (see Sections .., .. and . for corresponding Graph Cut literature) as well as topological constraints [HXP, BV], while Level-set methods too were shown to lend themselves to user-steered segmentation [YPH + ].

Simultaneous multiobject segmentation approaches are superior to those formulated as a sequence mono-object segmentations because they allow joint specifications of object properties and interrelations and entail simultaneous optimization with respect to these criteria. As such, they raise questions neither on the best segmentation sequence to follow nor on how to avoid the propagation of errors on individual segmentations [FAB]. Multiobject generalizations of the basic Level-set method usually present problem or application-dependent tradeoffs in either representing every object independently with its own internal and external speeds and relying on soft coupling terms to avoid overlaps and gaps at the expense of high memory requirements as the number of objects increases, or using so-called multiphase representations with a number of level-sets inferior to the number of objects which reduce storage requirements and ensure the absence of overlaps and gaps at the expense of reduced flexibility in the specification of internal and external speeds which often arise in practice [BPB]. By contrast, the basic Graph Cut segmentation method can be generalized to multiobject settings without modifying the form of the cost function and without representation overhead, although the complexity of the underlying graph structure is directly proportional to image size. In this case, the multilabel optimization of the cost function is expressed as a tightly coupled sequence of maximum-flow optimizations [BVZ].

Although the continuous level-set representation and the associated evolution strategies have a number of attractive properties, such as subpixel accuracy and natural handling of topological changes, the solution of underlying partial differential equations has to be addressed carefully in order to ensure numerical stability which often limits  contour evolution speeds making Level-set methods rather slow in practice, especially in multiobject formulations [TT]. By comparison, maximum flow algorithms underlying Graph Cut methods are deterministic polynomial-time algorithms which are numerically stable. For a particular class of cost functions covering a wide range of likelihood and prior formulations frequently arising in segmentation applications [KZ], Graph Cut methods produce provably good approximate solutions in multiobject and global optima in mono-object segmentation. The relative ease and the efficiency of multiobject formulation along with the speed and attractive numerical properties of associated optimization algorithms thus justify our choice to rely on Graph Cut methods as the foundation for our semiautomatic multiobject segmentation approach.

. The mission statement and the structure of dissertation

As we have seen, inherent limitations of quality and information content in medical images, their dependence of acquisition protocols and conditions coupled with tissue variability across patient populations make designing automatic general-purpose segmentation algorithms difficult, if not impossible. Therefore, our aim is to develop a general Graph Cut-based semiautomatic multiobject segmentation method, principally for medical images, that can be adapted to a particular application with minimal effort according to any or all of the following criteria and tradeoffs that exist between them: the number of objects to segment, the required anatomical level of detail, interrelations between objects to be respected, the accuracy of object and boundary retrieval, available computing resources of runtime and memory budget. The projected use is in routine applications involving low-resolution D to high-resolution D images on moderate to fairly hard multiobject segmentation tasks of up to few dozens of objects at a time. With the goal of minimizing user involvement and thus making the method viable in practical situations, we choose object-region style user interaction limited to initialization, allowing postsegmentation refinement if necessary. With a trained expert in mind at the human-computer interface, user input will also be relied upon as the principal source of information to collect intensity statistics of targeted objects as well as to constrain the assignment of user-marked pixels to some object.

Among segmentation criteria to adapt to a given application context, we address, in particular, the introduction of a fairly general graph-based prior model capturing spatial relationships between objects that, on one hand, is easy to specify for different relationships such as structure adjacency and distance, and on the other hand allows efficient optimization of resulting cost functions with multilabel Graph Cut algorithms. Targeting multiobject segmentation applications on high-resolution D images, we also address runtime and memory efficiency of Graph Cut segmentation with an imageadaptive clustering method paying particular attention to the robustness of segmentation and the invariance of its theoretical properties to the reduced data representation. The following list opens a window on the work's landscape summarizing its major contributions on a high level.



Chapter : Theoretical background, presents and discusses the theoretical foundations of our work, in particular the theory of Bayesian labeling and Markov random field modeling in image segmentation, concepts of flow networks, maximum flow techniques and Graph Cut algorithms in binary and multilabel optimization.

Chapter : A structural prior model for Graph Cut segmentation, investigates extensions of one of the most widely used prior models in Graph Cut segmentation, the Potts model, to account for multiple boundary interactions in multiobject segmentation allowing to minimize resulting multilabel cost functions by Graph Cut algorithms. Our main contribution is a structural prior model, called the vicinity prior, which captures the spatial configuration of objects. It is defined as shortestpath constraints on a graph model of interobject adjacency relations [KVDP].

Chapter : Reducing Graph Cut complexity, addresses runtime and storage efficiency of Graph Cut segmentation on high-resolution input and proposes to improve them both by clustering image voxels prior to segmentation by an imageadaptive centroidal Voronoi tessellation and reformulating Graph Cut segmentation according to the graph of clusters such that it is independent of the clustering resolution [KVDP]. The chapter also gives the results of a comparative evaluation of the clustering method with competing ones confirming its benefits in terms of runtime and quality of produced image partitions.

Chapter : Evaluations and applications, presents and discusses results of qualitative and quantitative evaluations of the proposed segmentation approach of multiobject Graph Cut segmentation on synthetic, simulated and real medical images confirming its advantages over its Potts prior-based counterpart, also showing that, compared to the application of segmentation on voxels, the clustering step improves both overall runtime and memory footprint of the segmentation process up to an order of magnitude without compromising the segmentation quality.

Chapter : Conclusions and perspectives, concludes the work outlining its major results and contributions and discusses future venues of research.

. Remarks on notation

A few remarks on general conventions we follow in spelling and notation are in order. We capitalize proper nouns when they refer to a particular method, an algorithm or a problem, designating them with small letters otherwise, for example, when speaking of the Maximum flow problem and maximum flow algorithms. In some contexts, we use the terms voxel and pixel interchangeably when such use does not affect the precision and the clarity of presentation. To simplify asymptotic notation of algorithm complexity, instead of O(|N |), we write O(N ), where N is some set the cardinality of which is proportional to the problem size. As it is customary in many treatments of discrete mathematics, we use braces to denote an unordered pair of indices or mathematical objects i, j as {i, j} using parentheses for ordered pairs (i, j).



Chapter 

Theoretical background

. Introduction

I

n this chapter, we present and discuss the theoretical foundations of our work for a fuller understanding and appreciation of our contributions. The presentation of methods and underlying theories is concise yet includes sufficient detail to enable unfamiliar readers to follow subsequent chapters without having to consult external references.

In Section . we present the theory of Bayesian labeling and Markov random field modeling in image segmentation discussing underlying probability distributions and how energies describing these models can be solved. Section . presents the theory of flow networks and discusses the algorithms that solve the Maximum flow problem. Maximum flow techniques constitute the foundation of our energy optimization method of choice, Graph Cut algorithms, which are discussed in Section ..

Throughout the presentation of methods, related literature is discussed in the appropriate context with state-of-art closely related to our contributions relegated to later chapters. The choice of a method over another for our problem is discussed and argued.

. Image segmentation as Bayesian labeling

We shall formulate segmentation as a labeling problem, defined as the assignment of a label from a set of m labels L = {l  , . . . , l m }, representing objects to be segmented, to each of the variables in a set of n variables, corresponding to image primitives indexed by S. In this chapter as well as in the following one, we shall assume that these primitives correspond to image voxels. Later in Chapter  we shall redefine them to correspond to clusters of image voxels. Assume that each variable i ∈ S is associated with a node in an undirected graph on image voxels G = S, E , where pairs of adjacent image voxels define the set of edges E = {i, j} | i, j ∈ S . Consequently, the neighborhood of a node i is defined as N i = j | j ∈ S, ∃ {i, j} ∈ E . Note that i N i because there is no edge from a node to itself. Although it can be arbitrarily complex, we assume that only immediate  neighbors of a node i constitute its neighborhood N i . Thus, the neighborhood of a variable is defined according to the directions parallel to the planes of a D image to get a -connected (-connected in D) system. When voxels on diagonals are also taken into account, the neighborhood becomes a -connected (-connected in D) system.

An assignment of labels to all variables is called a configuration, and is denoted by ∈ L, where L = L n is the space of all possible configurations. The assignment of a label to a single variable i ∈ S is denoted by i . The space of admissible configurations, which corresponds to possible segmentations, can be identical to L or, if user-specified constraints are taken into account, to a subset of L.

A label configuration ∈ L can also be expressed in probabilistic terms. Let X = {X  , . . . , X n } be a set of random variables on S taking its values in L. A configuration can be defined as the joint event X  =  , . . . , X n = n , abbreviated as X = . Defining a joint probability distribution Pr(X = ) in the space of configurations L, however, is rather difficult. Markov random field (MRF) theory [Li] is usually resorted to in order to define such joint probability distributions from interactions between neighboring variables which can be specified conveniently and meaningfully.

Given a joint probability distribution Pr(X), we recall that the set of random variables X is an MRF with respect to the graph G if it satisfies the following condition:

Pr(X i = i | X j = j , j ∈ S, j i) = Pr(X i = i | X j = j , j ∈ N i ). (.)
This is called the Markov local property. It states that the probability of a label assignment to a variable is conditionally dependent on assignments to its neighbors only; in other words, the state of a variable depends only upon the state of its neighbors. It still is not a trivial task to specify the set of conditional probabilities for an arbitrary joint probability distribution to satisfy the Markov local property (.). Fortunately, the Hammersley-Clifford theorem [Bes] provides a solution. According to this theorem, a positive joint probability distribution of label configurations Pr(X) >  satisfies the Markov local property with respect to G if and only if it is a Gibbs distribution with respect to G, that is, if it can be factorized over the cliques  of G.

A Gibbs distribution of configurations with respect to the graph G is given by:

Pr( ) = Z - × e - t U ( ) (.)
where Z is a normalizing constant, t is a hyperparameter called the temperature and U ( ) is an energy function specified as a sum of functions of label configurations, called clique potentials, defined on up to an arbitrary order of graph cliques

C = C  ∪ C  ∪ . . . : U ( ) = c∈C V c ( ) = {i}∈C  V  ( i ) + {i,j}∈C  V  ( i , j ) + . . . . (.)
Thus, the joint probability distribution Pr(X = ) can be conveniently defined by specifying appropriate clique potentials that encode prior domain knowledge about plausible configurations through local interactions between labels assigned to neighboring sites. Hence, Pr(X) is called a prior probability distribution.

 A clique in an undirected graph is a subset of vertices, each pair of which is connected by an edge.

The order of a clique is the number of its vertices.



Since we concern ourselves with finding the "best" segmentation with respect to some user-provided evidence, we are interested in the probability distribution of configurations Pr(X = | O) conditional on some, usually image-derived, observation O. This is called a posterior probability distribution. The observation constitutes a combination of spatially distinct records at each variable, for example, observed image intensity values. Evidence from observation is defined as the probability distribution Pr(O | X = ) and is called a likelihood distribution. Assuming that observations are local and spatially independent, it is reasonable to define likelihood probabilities at each variable Pr(O i | X i = i ) to be conditionally independent. Thus:

Pr(O | X = ) = i∈S Pr (O i | X i = i ) . (.)
Note that the posterior distribution Pr(X = | O) gives the probability of all possible configurations. The "best" configuration with respect to user-provided evidence can be considered the one with the highest probability given the observation, which corresponds to the highest mode of the posterior probability distribution:

* = arg max ∈L Pr (X = | O) . (.) Pr(X = | O) = Pr(O | X = ) Pr(X = )/ Pr(O).
Since Pr(O) does not depend on and therefore does not affect the MAP solution, the posterior probability is proportional to the joint likelihood and prior probabilities:

Pr(X = | O) ∝ Pr(O | X = ) Pr(X = ). (.)
Substituting the right-hand side of (.) from (.) and (.), taking the logarithm of both sides, ignoring the constant term ln z - and multiplying both sides by -t we get:

-t ln Pr(X = | O) ∝ t i∈S -ln Pr(O i | X i = i ) + c∈C V c ( ). (.)
Given likelihood and posterior probabilities, the negative logarithm of the posterior can be regarded as an energy function of label configurations, thus we write (.) as:

E( ) = t i∈S -ln Pr(O i | X i = i ) + c∈C V c ( ). (.)
Recall that the negative logarithm is a monotonically decreasing function. This means that the negative logarithm of the posteriorln Pr(X | O) achieves its minimum at the maximum posterior probability. Thus, the MAP solution can be equivalently found by: * = arg min ∈L E( ). (.)



Clearly, the maximum posterior probability configuration (.) is identical to the minimum posterior energy configuration (.).

Retaining only second-order clique potentials from the prior energy sum, we rewrite (.) as follows:

E( ) = t i∈S D i ( i ) +   i∈S j∈N i V i,j ( i , j ). (.)
This is the actual energy function we minimize to obtain the best segmentation with respect to user-supplied evidence and constraints, and some problem domain-dependent prior. Note that since either (i, j) or (j, i) is in the set of second-order binary cliques C  , either V i,j ( i , j ) or V j,i ( j , i ) should be in the prior energy sum (.). We consider that V i,j ( i , j ) = V j,i ( j , i ), and in order to use the double-sum notation over variables and their neighborhoods in (.), we multiply the double-sum by   .

.. Likelihood-based data cost and user constraints

The unary energy term  D i ( i ) in (.), which is also called a data term reflecting the fact that it is often derived from observed data, measures the cost of assigning a label to the variable i ∈ S so that unlikely label assignments incur higher penalties. This property stems from its definition as a negative log-likelihood function (.), since the higher a likelihood probability is, the lower its negative logarithm would be.

In many MAP-MRF formulations, the observation is taken to be the spatially distinct records of image intensity values at each variable, thus the likelihood probability for a variable i ∈ S given the label assignment i becomes Pr

(O i = I i | X i = i ).
A variety of approaches are followed to estimate such conditional probability distributions. For example, they can be represented as: object-specific intensity histograms [BF], kernel density estimates [LPPS], Gaussian distributions [GTC] and Gaussian mixture models (GMM) with hard [RKB] and soft membership probabilities [ZBS]. The latter approach sometimes follows an iterative Expectation-Maximization style two-step estimation [DLR] whereby, given an initial estimate of GMM parameters, the MAP solution is used in the expectation step to calculate membership probabilities, then GMM parameters are calculated in the maximization step in a standard fashion and the process is repeated until some convergence criterion is met. Such iterative schemes need to be carefully designed in order to ensure numerical stability and algorithm convergence [RRBCK]. In more recent conditional random field (CRF) image segmentation approaches [GRC + , FVS, LLB + ], observations often take the form of feature vectors over a set of different image features of intensity, texture and shape statistics [Low]. For each label class, feature vectors are constructed from ground-truth data and fed to classifiers, the output score of which, normalized as necessary, is used to estimate likelihood probabilities.

 Likelihood-based unary terms in (.) should not be confused with prior-based unary energy formulations frequently encountered in the MAP-MRF literature which stem from the first-order clique potentials of the prior energy sum (.). In the context of segmentation, such terms are derived, for example, from a probabilistic atlas of tissue location [STAG], a shape similarity measure [VM] or a distance transform of a shape template [KRBT].



In all results reported in our work, conditional probability distributions of intensity Pr(I | L = l) for every label are estimated as Gauss-smoothed and normalized intensity histograms from user-supplied "seeds", described in Section ., for every structure to be segmented. We could have employed a more elaborate approach, such as kernel density estimation [Par], to estimate likelihood densities. However such schemes require a choice of parameters, a good setting of which often entails a separate selection mechanism with no guarantees of yielding substantially better estimates in practice. We therefore have opted for a straightforward nonparametric approach.

Unary terms of the energy (.) can also encode user constraints on acceptable configurations. Recall that user-supplied seed labels indicate image voxels which are known to be part of some object. In addition to their use in deriving likelihood densities, in similar spirit to [BF], we take advantage of such seeds in reducing the solution space of the labeling problem by prohibitively increasing the cost of invalid configurations. In order to ensure the attribution of the respective label type to seed-marked voxels, the data term should create a maximum bias for the seed's label and, at the same time, incur the maximum penalty for the attribution of all other labels. Denoting the set of all seeds by H = {h l i | i ∈ S, l ∈ L}, we define the data term as follows:

D i ( i ) =            , ∃h k i ∈ H, k = i (.a) ∞, ∃h k i ∈ H, k i (.b) -ln Pr(I v | i ), h k i ∈ H, ∀k ∈ L (.c)
where h k i is a seed of label type k ∈ L attributed to voxel i ∈ S and I v is the gray-level of voxel v. Note that (.) is a continuous function where (.a) and (.b) correspond to negative log-likelihood values at maximum and zero probability respectively. In practice, we represent the infinite cost in (.b) by a large constant several orders of magnitude (  -  ) greater than the value of (.c) corresponding to the smallest estimated likelihood probability, so that for practical settings of temperature t > , (.a) and (.b) amount to imposing hard constraints for acceptable configurations. By comparison, in order to ensure the invariance of user constraints with respect to temperature settings, [BF] avoid multiplying (.b) by t in (.) and define a constant K representing its cost such that prior terms do not accidentally impose a labeling decision on corresponding variables. This, however, does not rule out the possibility of (.c) taking a value greater than K for a sufficiently small likelihood probability. In our formulation, as t greatly decreases, temperature-attenuated terms (.b) approach the order of magnitude of prior terms until the sum of data terms in (.) vanishes at t = , even though such a temperature setting is hardly interesting in practice.

In subsequent refinements of an initial segmentation, we allow the user to introduce additional "corrective" seeds to constrain the label preference of incorrectly segmented image regions. In order not to disrupt initial appearance statistics, we do not take such seeds into account for likelihood estimation.



.. Prior models for segmentation

Recall that the binary terms V i,j ( i , j ) in (.) encode prior information on pairwise interactions between labels assigned to pairs of neighboring sites. They encourage the spatial consistency of labeling with respect to a problem-dependent reference model, hence they are sometimes called spatial regularization terms.

Many classes of pairwise interaction models have been proposed in the MRF modeling literature for different image processing and computer vision problems [Li]. For image segmentation, the pairwise interaction has to encourage piecewise-constant labeling, which is the partitioning of the image into several regions where voxels within individual regions are identically labeled. Furthermore, it has to be discontinuity preserving, meaning that, on one hand, it should attribute identical labels only to image regions which do not contain intensity edges, and, on the other hand, it should not excessively penalize sharp label transitions. Convex functions, which have a number of advantageous properties for mathematical optimization [BV], unfortunately are not discontinuity-preserving because they do not have a bound on the largest possible penalty and potentially overpenalize large label differences. Optimizing an energy function with a convex pairwise potential would therefore lead to oversmoothing at object boundaries where abrupt transitions occur [BVZ], hence the use of convex priors in segmentation is limited to special contexts, such as the oversegmentation algorithm in [IGb]. Smooth convex priors however are desirable in some vision problems such as image restoration [BK] where labels correspond to image gray-levels and stereo matching [IGa] where labels represent disparities.

The earliest discontinuity-preserving piecewise-constant prior model is the Ising model in binary labeling where |L| = , usually L = {, }. It has first been introduced in statistical mechanics as a model of ferromagnetism and later used in image processing and computer vision by analogy between images and statistical mechanics systems [GG]. Given a pair of variables i, j ∈ S where j ∈ N i , the Ising model is defined as:

V i,j ( i , j ) =  -δ( i , j ) (.)
where δ is the Kronecker delta function. The identically-defined generalization of the Ising model for multilabeling problems where |L| >  is called the Potts model. We can readily see that both the basic model and its generalization incur a single level of penalization corresponding to the attribution of different labels to pairs of neighboring variables and incur no penalty for the identical labeling of such pairs. This penalization scheme, in fact, stems from the very nature of binary labeling where only a single type of boundary exists; the one between image regions labeled by either label.

Although the Potts prior is the most widely used spatial regularization model in MAP-MRF formulations of image segmentation and continues to enjoy popularity in recent CRF-based approaches [FVS, SWRC], many problem-specific pairwise models can be encountered in the literature, especially in the context of multiobject segmentation. For example, an ad-hoc extension of the Potts prior model to two penalization levels with respect to a label class cooccurrence model can be found in [CCB + ]. A shape template's distance transform has been used by [FZ] to bias the placement  of segmentation boundary near the shape's contour. Spatial configuration prior is also quite common in the context of multiobject segmentation, where inter-object geometric and spatial relationships are often defined as pairwise relations. Such use has appeared in the context of part-based object detection [WS] and scene layout estimation and segmentation where spatial relationships are either predefined on a limited number of classes [LVS] or learned from annotated training data [LLB + ].

In Chapter  we propose and investigate a principled approach to extend the basic discontinuity-preserving piecewise-constant prior as defined by the Potts model to account for multiple types of boundaries in multiobject segmentation, incurring multiple levels of penalization, furthermore allowing the application of a particular energy minimization algorithm which converges rapidly to a provably good local minimum.

.. Solving posterior energies

Having addressed common prior models and likelihood estimation approaches for MAP-MRF image segmentation, we now turn our attention to the minimization of the energy function (.). In the following, we shall discuss several widely used energy minimization methods without delving into the details of algorithms. Our method of choice is presented in detail in subsequent sections of the current chapter.

Since many early vision problems are naturally expressed as pixel-labeling problems in terms of a likelihood probability and a prior distribution, energy functions of the form (.) frequently arise in practice. However, the minimization of such, often nonconvex, energies in high-dimensional spaces is not a trivial task. In fact, [BVZ] show that the global minimization of (.) is NP-hard when |L| >  even with the simple discontinuity-preserving Potts prior model, since it reduces to the multiway cut problem which is known to be NP-complete [DJP + ]. This unfortunately implies that many similar energy formulations are intractable. Therefore, research has focused on efficient approximation algorithms which produce provably good local solutions. Having a bound on the local minimum is important, because when it is arbitrarily far from the global minimum, it would be difficult to determine the cause of the algorithm's unsatisfactory results, which could be due either to a bad energy formulation or a poorly performing minimizer.

Early energy minimization algorithms use an elementary search space exploration strategy where only a single pixel can change its label at a time. This is sometimes called a standard move. Starting with an initial estimate of the labeling, for each pixel, the algorithm chooses the label that gives the largest energy decrease either in a deterministic greedy fashion, like the Iterative Conditional Modes (ICM) algorithm [Bes], or according to probabilistic metaheuristics, like the Simulated Annealing (SA) algorithm [KGV], which both repeat the selection process for all pixels until convergence to a local minimum is achieved. The SA algorithm in particular has enjoyed popularity because of its implementation ease and ability to optimize an arbitrary energy. Theoretically, it would eventually find the global minimum if run for an infinite amount of time. In practice however, minimizing an arbitrary energy in finite time requires an  exponential cooling schedule to ensure convergence to a good local minimum [Cat], which means that the algorithm can be rather slow to terminate.

ICM and SA remained in mainstream use in vision until late 's when energy minimization approaches had a renaissance with the introduction of flow networks techniques from combinatorial optimization for the exact global minimization of some widely used classes of energy functions in early vision. In particular, maximum flow algorithms of low-order polynomial time were used to find the global minimum in binary image restoration with the Potts prior [GPS] and in stereo matching where |L| >  with linearly ordered labels and an absolute difference prior [IGa] which was later extended to handle arbitrary convex priors [Ish]. The binary exact optimization scheme in particular proved to be rather interesting in early vision and opened the door to a multitude of applications like interactive image and video segmentation, stereo, multiview reconstruction and texture synthesis (refer to the comprehensive list of references in [BF]). Subsequently, [BVZ] generalized the maximum flow-based binary optimization scheme to |L| >  for arbitrary labels and metric pairwise priors. In contrast to ICM and SA, their approximate labeling algorithm employed a search space exploration strategy, called an expansion move, whereby an arbitrarily large number of voxels can switch to a given label α simultaneously. The Expansion Moves algorithm, as we will refer to it henceforward, finds a local minimum such that no expansion move, for any label α, yields a labeling with lower energy. The local minimality criterion with respect to expansion moves is so strong that there are many fewer such minima in high-dimensional spaces compared to standard moves. [BVZ] prove that any labeling locally optimal with respect to these moves is, in fact, within a known factor of the global minimum. Usually, maximum flow-based energy minimization approaches in vision are collectively referred to as Graph Cut methods.

More recently, the revival of older belief propagation methods [Pea] has led to the development of the Sequential Max-Product Tree-Reweighted Message Passing algorithm (TRW-S) [Kol]. It is often seen as a rival to the Expansion Moves and the Loopy Belief Propagation algorithms (LBP) [WF] the latter of which has strong local minimality properties analogous to those of the former. [MYW] has shown that TRW-S is capable of finding the global minimum for certain instances of the stereo problem. This algorithm can be applied to arbitrary energies on arbitrary label sets, although, like LBP, sometimes it might oscillate and fail to converge.

In order to understand the tradeoffs between different energy minimization methods, [SZS + ] have devised a set of benchmarks on common vision problems and used them to compare the solution quality and runtime of ICM, LBP, Graph Cut and TRW-S algorithms. They conclude that while the latter three come very close to the global minimum on different multilabeling problems, the Expansion Moves variant of Graph Cut algorithms consistently outperforms others on nonconvex energies in terms of both speed and quality of obtained solutions. [BVZ] had earlier benchmarked SA and Expansion Moves algorithms and made similar conclusions. The comparative study of [SZS + ] also echoes the findings of [MYW] in that better minimization methods are unlikely to produce significantly more accurate labelings, since the global minimum,  when known, is no more accurate than the ground-truth labeling, the energy of which is furthermore frequently higher than that of the former. This suggests, as pointed out by [MYW], that often the problem is not in the optimization algorithm but in the energy formulation. Therefore, research should focus on finding better energy models.

Due to its optimality and runtime performance, we use Graph Cut algorithms for segmentation energy minimization. Closely following the method of [BF], we use a maximum flow-based binary optimization scheme in mono-object segmentation where |L| =  and apply the Expansion Moves algorithm of [BVZ] for multiobject segmentation with |L| > . In the remaining of this section, we shall give the necessary and sufficient conditions for the applicability of these methods to minimize the energy (.). Theoretical foundations and algorithmic details of maximum flow and Graph Cut methods are presented in Sections . and . respectively. Readers who are familiar with flow networks and Graph Cut methods can skip to Chapter  directly.

While the maximum flow-based binary energy minimization method is capable of finding the exact MAP solution [GPS], the Expansion Moves algorithm produces a locally-optimal labeling ˆ within a known factor of the global minimum * , that is:

E( ˆ ) ≤ cE( * ), c =  × max V i,j (a, b)/min V i,j (a, b), ∀i, j ∈ S, ∀a, b ∈ L. (.)
For example, c =  for the Potts model. As mentioned previously, in the original work of [BVZ], the Expansion Moves algorithm was shown to be applicable to energies of metric pairwise priors. [KZ] subsequently characterized the classes of energy functions that can be minimized via Graph Cut algorithms and relaxed the metricity requirement. In particular, they showed that the Expansion Moves algorithm can be applied to minimize (.) with |L| >  if and only if the pairwise term is submodular, that is, if the following condition holds ∀i, j ∈ S, ∀a, b, c ∈ L:

V i,j (a, a) + V i,j (b, c) ≤ V i,j (a, c) + V i,j (b, a). (.)
For L = {, }, the submodularity condition reduces to:

V i,j (, ) + V i,j (, ) ≤ V i,j (, ) + V i,j (, ). (.)
Submodular set functions and their properties and applications in combinatorial optimization are known to mathematicians for a long time [BKR]. Importantly, optimization problems of functions having such properties are polynomial-time solvable. In fact, submodularity can be regarded as a property of discrete functions analogous to the convexity property of continuous functions [Mur]. However, polynomial algorithms for minimizing arbitrary real-valued submodular functions have high complexity. For example, the algorithm proposed by [Orl] runs in O(n  γ +n  ) time where γ is the time to evaluate the submodular function. Therefore, Graph Cut-minimizable functions can be regarded as a subclass of submodular functions which can be minimized with much faster algorithms. The optimization of submodular functions remains an active research field.



. Flow networks

Flow networks are tools from combinatorial optimization used to solve problems that can be modeled as a steady-rate material flow through a directed capacitated network from a source node, where the material is produced, to a sink node where it is consumed. For example, we can think of liquid flows through pipes, and also electrical current and information through electrical and communication networks. Directed edges in such a network represent conduits each having a capacity indicating the maximum rate at which material can flow through it. Network nodes, other than the source and the sink, are conduit junctions through which material flows without collecting.

An important result in flow networks theory, due to [FF], is the duality of Maximum flow and Minimum cut problems. The maximum flow is the greatest rate at which material can be sent from the source to the sink without violating any capacity constraints. A cut of the flow network is a partition of its nodes into two disjoint sets such that the source is in one set and the sink is in the other. The minimum cut has the minimum capacity among all cuts, measured as the sum of the capacities of edges going from the source set to the sink set. Importantly, the minimum cut capacity is equal to the maximum flow value which can be computed by efficient algorithms as we shall see in this section. Let us first define flow networks and their properties formally.

A flow network is a directed weighted graph G = V , E with no self-loops where we distinguish two special nodes, called terminals, the source s and the sink t, sometimes also called the target. Each node in V lies on a path from the source to the sink, therefore the graph G is connected and |E| ≥ |V | -. In the special case of grid graphs with -connected and -connected neighborhood systems in D and D, |E| ≥ |V | and |E| ≥ |V |, respectively. A flow network G is equipped with a positive capacity function of edges; c :

E → R + so that each edge (u, v) ∈ E has a positive capacity c(u, v) ≥ . Note that c(u, v) =  when (u, v) ∈ E.
Given a flow network G = V , E and a capacity function c, a flow f in G, is a realvalued function f : V × V → R satisfying the following properties.

Capacity constraint: ∀u, v ∈ V , f (u, v) ≤ c(u, v) (.a) Skew symmetry: ∀u, v ∈ V , f (u, v) = -f (v, u) (.b) Flow conservation: ∀u, v ∈ V -{s, t}, v∈V f (u, v) =  (.c) Obviously f (u, v) =  when (u, v) E.
The value of flow in the network G is given by:

|f | = v∈V f (s, v). (.)
In the Maximum flow problem, given a flow network G we wish to find a flow of maximum value given by (.) respecting capacity constraints (.a) and the flow conservation property (.c). The skew symmetry property is merely a notational convenience indicating the direction of positive flow. Figure .a illustrates a flow network having a valid flow. Two general methods exist for solving the Maximum flow problem; the Ford-Fulkerson method, sometimes also called the Augmented-paths method, and the Push-relabel method. Several algorithms of different running times implementing either approach have been proposed in the literature. In the following two sections, we outline each method providing more detail on the first since we use an algorithm that falls in this category. It is worth mentioning that the Maximum flow and Minimum cut problems can also be stated as two primal-dual linear programs [Law]. However, such a formulation can be impractical for the large number of variables associated with vision and image processing applications.

 t s v  v  v  v    ,    ,   ,   ,     ,   ,  ,  ,  ,   t s v  v  v  v            (a) (b) (c)        t s v  v  v  v    ,     ,    ,     ,   ,  ,  ,  ,  

.. The Ford-Fulkerson method

The Ford-Fulkerson method is characterized by three principal ideas upon which algorithms that implement it depend. They are furthermore essential to the Max-flow Mincut theorem which ensures the correctness and the termination of such algorithms.

In the following, we present these concepts before outlining the general strategy of finding a maximum flow by the Ford-Fulkerson method.

Residual networks. Starting from |f | = , the Ford-Fulkerson method iteratively increases the value of flow until no more flow can be sent from the source to the sink due to capacity constraints. The method maintains a residual network G f whose topology is identical to that of G but the capacities of its edges represent the residual capacities of corresponding edges in G. The residual capacity of an edge (u, v) in G is the amount of additional flow we can push from u to v before exceeding c(u, v). It is given by:

c f (u, v) = c(u, v) -f (u, v). (.)
The residual network G f provides a roadmap for increasing flow in the original network with a simple rationale. Consider, for example, the edge (s

, v  ) in Figure .a with c(s, v  ) =  and f (s, v  ) = . We can increase flow on (s, v  ) by c f (s, v  ) =  before ex- ceeding the capacity constraint on the edge. Consider the negative flow f (v  , v  ) = -. In this case, the residual capacity c f (v  , v  ) =  is greater than the capacity c(v  , v  ), meaning that we may cancel the flow f (v  , v  ) =  by pushing a flow of  from v  to  v 
, and subsequently push a flow of  before violating the capacity constraint c Augmenting paths. Given a flow network G = V , E and a flow f , an augmenting path p is a simple path  from the source to the sink in the residual network G f . We can increase the flow on an edge (u, v) of an augmenting path by up to c f (u, v) without violating the capacity constraint on whichever of (u, v) or (v, u) is in the flow network G. The maximum amount by which we can increase the flow through the path p is called the residual capacity of p and is given by

(v  , v  ). Note that c f (u, v) = c f (v, u) =  when neither edge is in G, therefore an edge (u, v) ap- pears in G f if at least one of (u, v) and (v, u) are in G. Formally then, given a flow net- work G = V , E and a flow f , the residual network of G induced by f is G f = V , E f where E f = {(u, v) ∈ V × V | c f (u, v) > } and |E f | ≤ |E|. To illustrate,
c f (p) = min{c f (u, v) | (u, v) ∈ p}. Highlighted edges in Figure .b indicate an augmenting path with a residual capacity of c f (p) = .
Cuts of flow networks. The last concept we introduce before presenting the Max-flow Min-cut theorem and outlining the Ford-Fulkerson general procedure, is that of cuts of flow networks. As mentioned in the beginning of this section, a cut of a flow network G = V , E , denoted as (S, T ), is a partition of its set of nodes V into two disjoint sets S and T = V -S such that s ∈ S and t ∈ T . The capacity of a cut is given by:

c(S, T ) = u∈S v∈T c(u, v). (.)
Obviously, the minimum cut of the network is the cut that has the minimum capacity among all cuts of the network. • f is a maximum flow in G.

• There are no augmenting paths in the residual network G f .

• |f | = c(S, T ) where (S, T ) is the minimum cut of G.

We refer the reader to [CLRS] for the proof. Intuitively, the Max-flow Min-cut theorem states that the value of the flow in a network is bounded from above by the capacity of its minimum cut, and provides the termination condition for maximum flow algorithms following the Ford-Fulkerson method; the algorithm terminates and the flow is maximum if and only if its residual network contains no augmenting paths.

With all concepts in place, the general Ford-Fulkerson method of finding a maximum flow boils down to a simple procedure. Starting with In practice, the Maximum flow problem arises with integer capacities. If capacities are rational numbers, they can be converted to integers with an appropriate scaling transformation. With irrational capacities however, the Ford-Fulkerson method might iterate infinitely, finding augmenting paths of increasingly lower capacity. Worse, the infinite sequence of augmentations may not even converge to the maximum flow value. [Zwi] illustrates the smallest network on which this problem occurs. Therefore, irrational capacities have to be approximated by rational numbers or integers.

f (u, v) =  ∀(u, v) ∈ E
Implementations of the Ford-Fulkerson method address two questions; how augmenting paths are found, and in what order such paths are processed. Denoting the maximum flow by f * , a straightforward implementation executes |f * | iterations at most, assuming that the flow increases by one unit at each iteration. If augmenting paths in the residual network are found using either a breadth-first or a depth-first search with a complexity of O(V +E f ) = O(E f ) with no regard to path length, the total running time of the method would be O(|f * |E f ). The Edmonds-Karp algorithm [EK] improves the bound on the number of iterations to O(V E f ). Attributing unit weight to each edge in the residual network, in each iteration the algorithm choses the augmenting path corresponding to the shortest path from the source to the sink found by breath-first search. Therefore the algorithm runs in O(V E  f ) time. An algorithm due to [AO] uses a scaling technique and runs in O(E  f log  C) time where C = max ∀(u,v)∈E {c(u, v)}. Starting with K =  log  C , and while K ≥ , the algorithm finds augmenting paths of capacity at least K at each iteration halving K at the end of it.

In an attempt to improve the empirical performance of the Ford-Fulkerson method in vision applications involving grid graphs, [BK] observe that building search trees to find augmenting paths in every iteration, which is a O(E f ) operation, is prohibitively expensive especially in D where usually |E| > |V |. They propose reusing such trees. Their algorithm maintains two search trees rooted at the source and the sink respectively. In each iteration, each tree is grown through nonzero residual capacity edges acquiring new children from graph nodes until a leaf node from a tree comes in contact with a leaf node from the other tree resulting in an augmenting path. Flow is pushed through this path in the familiar manner, which breaks search trees into forests. Next, tree structures are restored by assigning nodes that have been disconnected by augmentation to new parents via nonzero capacity edges. If no parent can be found, disconnected nodes are made free nodes. The algorithm terminates when search trees cannot grow due to their leaves being separated by zero residual capacity edges. At this point, the flow is maximum and the nodes of either tree fall in the source and the sink set respectively. The obvious drawback of the algorithm proposed by [BK] is that augmenting paths found are not necessarily the shortest. As reported, the processing order of nodes at tree growth and restoration stages significantly impacts the algorithm's empirical runtime. Authors proceed in a First-In-First-Out (FIFO) fashion which ensures that the first path found is the shortest, with no guarantees on subsequently found paths since search trees may change unpredictably during the restoration stage. In order to make augmenting paths as short as possible, distance-to-root information is stored in each node so that disconnected nodes can be assigned to parents closer to a root. Although no asymptotic bound is given, the algorithm has supralinear empirical time in vision problems involving D and D grid graphs and outperforms the Edmonds-Karp algorithm and two algorithms among the best performing Push-relabel variants, namely those with queued and highest node selection rules discussed in the next section.

We use the algorithm of [BK] in our maximum flow-based image segmentation applications. As we shall see in Chapter , in order to reduce the redundancy of voxel data and improve runtime performance, we formulate our segmentation method in terms of clusters of voxels. As shown in Figure ., the neighborhood connectivity of such irregular cluster graphs is higher than that of image grid graphs. For example, we can see that with a clustered D image where the number of clusters equals % of the number of voxels, the number of neighbor edges of the cluster graph are as high as six times the number of voxels. This means that any algorithm that finds augmenting paths with breadth-first search would perform worse than the algorithm of [BK].



.. The Push-relabel method

For the completeness of the presentation of maximum flow approaches, we outline the general Push-relabel method in this section pointing at common implementations. The Push-relabel method differs from the Ford-Fulkerson method in two primary ways. Firstly, the Ford-Fulkerson method examines the entire residual network in search for augmenting paths, and each such path gives rise to an additional stream of flow from the source to the sink until no more flow can be pushed. The Push-relabel method works in a localized fashion, examining a single node at a time and its neighbors in the residual network, thus the method readily admits parallel implementations. Secondly, while the Ford-Fulkerson method maintains the flow conservation property (.c) throughout its execution, the Push-relabel method allows each node u ∈ V -{s, t} to have a positive flow excess given by:

e(u) = v∈V f (v, u) ≥ . (.)
In order to determine how flow is pushed from an overflowing node to its neighbors, the method maintains a height attribute for each node. Although flow from a lower node to a higher neighbor can be positive, flow can be pushed only from a higher node to a lower neighbor. In the beginning, the height of the source and the sink are set to |V | and  respectively, and those of all other nodes start at  and increase as the algorithm progresses. The algorithm first sends the largest possible flow from the source towards the sink, equal to the capacity of the cut (s, V -{s}), which collects in intermediate nodes. In order to rid an overflowing node of excess flow, the algorithm increases its height to one unit above the height of the lowest of its neighbors with which it shares a nonzero capacity residual edge. This operation is called relabeling, after which the node can push excess flow to at least one of its neighbors. All flow that can reach the sink eventually arrives there and in order to make flow in the network a valid flow, undeliverable flows are drained back to the source by continuing to relabel nodes above the fixed height |V | of the source. When all excess flow is drained, the flow in the network is a maximum flow.

The basic Push-relabel method, due to [Gol], thus boils down to an initialization followed by a sequence of push and relabel operations, whichever is applicable, executed in no particular order. It has an asymptotic bound of O(V  E), which is better than that of the Edmonds-Karp algorithm mentioned above. [GT] propose two improvements on this bound; O(V  ) by maintaining overflowing nodes in a FIFO queue and always selecting the front node for discharging, that is applying push and relabel operations successively until the said node has no positive excess, and adding newly overflowing nodes to the rear of the queue, and O(V E log  (V  /E)) by using an efficient dynamic tree data structure to reduce per push-operation time. [CG] address practical aspects of implementing the Push-relabel method underlining the importance of using heuristics. They conclude that the best Push-relabel algorithm is one that repeatedly discharges the highest overflowing node. It has an asymptotic bound of O(V  √ E).



It is worth mentioning as we conclude this section that the asymptotically fastest deterministic sequential algorithm for the Maximum flow problem of which we are aware is due to [GR] which is not a Push-relabel algorithm. It uses the concept of blocking flows and associates shorter lengths with higher-capacity edges in order for shortest st paths to have high capacity and the algorithm to perform fewer iterations. The algorithm runs in O(min(V / , E / )E log  (V  /E) log  C) time where C is the maximum edge capacity. Its performance on flow networks arising in vision is unknown to us.

. Graph Cut methods

Having discussed the Maximum flow problem and algorithms to solve it, we return our attention to minimizing the energy (.) with Graph Cut algorithms. In Section .. we outlined the general properties of this family of energy minimization algorithms mentioning that they apply to energy functions of the form (.) with submodular pairwise energy terms (.), and that in practice they outperform popular energy minimization algorithms, in terms of both runtime and quality of produced minima, on high-dimensional nonconvex energies frequently encountered in vision. In this section we examine the algorithmic details and the properties of Graph Cut algorithms.

Maximum flow methods lie at the heart of Graph Cut algorithms. In fact, every such algorithm entails a method of representing the energy (.) as a flow network so that a cut on the network and its capacity correspond to a label configuration and its energy respectively, furthermore ensuring that the minimum cut of the network minimizes the energy among all possible cuts.

Graph Cut methods give two interpretations to cuts of flow networks: ) a geometrical interpretation, typically arising in binary problems, which views such cuts as optimal separating hypersurfaces in an N-D manifold embedding the graph equipped with a cut metric approximating any continuous Riemannian metric [BK, KB], and ) as a classifier, finding the optimal attribution of two or more label classes to a set of variables [GPS, BVZ, Ish]. In both cases, maximum flow algorithms constitute the numerical optimization tool of energies describing such hypersurfaces and classes.

Maximum flow methods, as we have seen in Section ., are inherently binary techniques; they produce an optimal bipartition of the nodes of a flow network in terms of the capacity of the cut separating the two sets. In the following section, we shall show that binary problems in vision like mono-object segmentation, binary image restoration etc. can be most naturally and directly formulated and efficiently solved by maximum flow algorithms. Importantly, for energies of binary variables, the MAP solution (.), that is the global minimum of the energy (.), can be obtained via Graph Cut methods in polynomial time matching that of the underlying maximum flow algorithm. Moreover, multilabel energies, discussed in Section .., can be minimized with Graph Cut methods by transforming the multilabling problem to a sequence of binary subproblems each formulated as a Maximum flow problem on an appropriate network.

 s t i j tD j () tD i () tD i () tD j ()  n   Figure .:
A flow network representing the energy (.) with binary variables x i ∈ {, }. The configuration induced by the blue cut is x i = x j =  and that of the red cut is

x i = , x j = 

.. Graph Cut in binary optimization

In this section, we discuss one of the most straightforward Graph Cut applications of mono-object segmentation showing how the minimization of the energy (.) with binary variables can be formulated as a Maximum flow problem. The original formulation is due to [GPS] in the context of binary image restoration. Their approach remained relatively unknown to the vision community until it was reintroduced and popularized by [BJ] in mono-object interactive segmentation. Our discussion in this section closely follows the method presented in the latter work.

Let us begin by defining the flow network. As before, let G = V , E be a directed connected graph with a positive capacity function of edges. Let S be the set of integers indexing the variables of the energy (.) and the corresponding image primitives which are considered to be image voxels in this discussion. The set of nodes of the graph G is V = S ∪ {s, t} where s and t are the source and the sink terminals of the flow network respectively. The set of edges E consists of two types of edges; those connecting nodes in S to either terminal, called t-links (for terminal links), and those connecting pairs of nodes in S, called n-links (for neighbor links). The latter are dictated by the neighborhood system of the energy (.). The set of edges is therefore given by

E = i∈S {(s, i), (i, t)}, i,j∈S {(i, j) | j ∈ N i } . An illustration of the graph G is given in Figure . for a one-dimensional image.
Let us formulate mono-object segmentation according to the Bayesian labeling framework presented earlier in Section .. For this binary problem, the set of labels is L = {, }. We assume that "" represents the background and "" the sought object. We take the observation to be image intensity values at each variable, and assume that likelihood densities are defined for both labels. Disregarding user constraints, data terms of the energy (.) are given by (.c). For pairwise energy terms, the Potts prior model (.) suffices to ensure the spatial consistency of this instance of binary labeling. Note that data terms are positive by definition and that if a variable i ∈ S prefers a label, say the object label "", then D i () < D i ().

The one-to-one correspondence between cuts of the graph G and binary configurations of variables in S is fairly obvious. A cut of G, which partitions its nodes into  the source S and the sink T sets respectively, can be interpreted as a binary configuration which assigns labels {, } to variables in S corresponding to nodes such that i =  if i ∈ S and i =  if i ∈ T . Clearly, any cut (S, T ) on G corresponds to a configuration . There remains to show that the capacities of the edges of G can be chosen so that the capacity of a cut is equal to the energy of the configuration induced by it. Given a cut (S, T ), we denote the set of edges separating the source and the sink sets by

E (S,T ) = {(u, v) | u ∈ S, v ∈ T }. Consider the graph in Figure ..
Cuts on G can include t-links and n-links, however for each node, exactly one t-link is in E (S,T ) . Thus, if i ∈ S, (i, t) ∈ E (S,T ) , and if i ∈ T , (s, i) ∈ E (S,T ) . Choosing c(s, i) = tD i () and c(i, t) = tD i (), every node i ∈ S contributes the data penalty associated with the chosen label to the t-link part of the cut capacity. Note that an n-link (i, j) is in E (S,T ) if and only if i ∈ S and j ∈ T , thus by setting c(i, j) = , every n-link (i, j) ∈ E accounts for the Potts prior penalty in the cut capacity. The capacity of a cut on the network G thus defined is given by:

c(S, T ) = t i∈S D i ( i ) + (i,j)∈E i∈S,j∈T  (.)
which is clearly equal to the energy of the configuration induced by the cut (S, T ) as defined by (.) with the condition i j for pairwise terms. The minimum cut of G therefore yields the MAP configuration * and its capacity is equal to the energy E( * ).

... General-purpose constructions for binary energies

The graph construction we outlined in Section .. assumes that unary energy terms are positive and that pairwise terms follow the Potts model. The natural question that arises consequently is: could we construct such graphs representing arbitrary energies? The seminal work [KZ] answers this question giving a necessary graph representability condition for functions that can be written as a sum of terms of arbitrary order. Importantly, it shows that this condition is sufficient for functions which can be written as a sum of up to tertiary terms. In our work, we formulate mono and multiobject segmentation with binary and multilabel energies written as a sum of up to binary terms, therefore, in the following we shall present the main result of [KZ] for such energies.

Theorem .. (Graph Representability). Let F be a function of n binary variables indexed by S written as:

F(x  , . . . , x n ) = i∈S F i (x i ) + (i,j)∈N F i,j (x i , x j ), (.)
where N ⊂ S × S is a neighborhood system containing ordered pairs (i, j) for which i < j. F is graph representable if and only if each term F i,j satisfies the following inequality:

F i,j (, ) + F i,j (, ) ≤ F i,j (, ) + F i,j (, ). (.)
Note that (.) is the submodularity condition (.) we mentioned earlier in Section ... To make the distinction between functions of sets and functions of binary

 s t i s t i F i () -F i () F i () -F i () s t i j B + C -A -D D -C A -C (a) (b) Figure .:
Graph representations due to [KZ] for (a) unary and (b) pairwise terms of the function (.). variables, [KZ] call it the regularity condition, hence binary functions satisfying (.) are called regular functions. In our discussion, we use the more widely used terminology of submodularity and submodular functions.

F i,j = F i,j (, ) F i,j (, ) F i,j (, ) F i,j (, ) = A B C D = A +   C -A C -A +  D -C  D -C +  B + C -A -D  
In the following, we shall outline the constructive proof of theorem .. showing that submodularity is a sufficient condition for graph representability by describing the graph that represents an arbitrary submodular function of the form (.). We refer the reader to [KZ] for the proof of the reverse direction as well as for graph constructions representing submodular functions written as a sum of up to tertiary terms. Another important theoretical result reported in this work is that the minimization of nonsubmodular functions of the form (.) is an NP-hard problem.

Given a cut (S, T ) on the graph we defined earlier in Section .., which represents a binary energy function with positive data and Potts prior terms, the energy of the configuration induced by the cut is exactly equal to the latter's capacity E( ) = c(S, T ). Observing that the configuration which minimizes the binary function (.) does not change if we add a constant to it, for the constructive proof of the theorem .. we shall relax the last equality requiring that F(x  , . . . , x n ) = c(S, T ) + K, where K is some constant. Obviously, F is exactly represented by a graph if K = .

Let F be a submodular binary function of the form (.) and G = V , E a flow network. We shall construct a subgraph representing each term of F separately and merge all subgraphs into G afterwards. This is justified by the additivity theorem outlined in [KZ] which states that the sum of two graph-representable functions is graphrepresentable. Thus, we can illustrate the construction on the minimal example of a pair of variables indexed by i, j. As before, the set of nodes of the graph G is V = {s, t, i, j} where each nonterminal node encodes a binary variable. For each term in F, we shall add one or more edges to E so that edge capacities are positive as required.



Consider first the unary terms of F. As illustrated in Figure .a, the subgraph corresponding to a single variable i contains the nodes {s, t, i}. If F i () < F i (), we add the edge (s, i) to the subgraph with capacity c(s, i) = F i () -F i (), otherwise we add the edge (i, t) with capacity c(i, t) = F i () -F i (). We can see that the constructed subgraph represents the function F i in both cases so that F i (x i ) = c(S, T ) + K with K = F i () in the former and K = F i () in the latter case.

Consider now the pairwise terms of F. The subgraph corresponding to a pair of variables i, j naturally contains the nodes {s, t, i, j}. It is convenient to represent such a term in tabular form according to the four combinations of its variables' values as in Table .. As shown, we rewrite the pairwise term in four terms and consider each for edge addition separately. The first term A is a constant, therefore no edges need to be added for it. The second and third terms are unary functions which depend on x i and x j respectively, therefore we can use the previous construction for unary terms to represent these terms. To represent the fourth term, we add the edge (i, j) with capacity

c(i, j) = B + C -A -D. Note that B + C -A -D ≥  since F is submodular. Figure .b illustrates the construction corresponding to the case A-C >  and D -C > . Note that in this case F i,j (x i , x j ) = c(S, T ) + F i,j (, ).
As mentioned previously, we obtain the graph G that represents the function F by merging separately constructed subgraphs for unary and pairwise terms and adding edge capacities together. For referencing purposes, Table . summarizes graph constructions for submodular binary functions written as a sum of up to binary terms.

.. Generalized Graph Cut for multilabel energies

Many interesting vision tasks can be formulated as labeling problems with |L| >  and solved by multilabel optimization of energies of the form (.). To mention few instances, the set of labels L can represent disparities in stereo matching, motion vectors in motion estimation, intensity levels in grayscale image restoration and object identifiers in multiobject segmentation.

We have reported in Section .. that minimizing nonconvex multilabel energies of the form (.) is NP-hard, even with the simple discontinuity-preserving Potts prior model. We are aware of only two Graph Cut methods for exact multilabel minimization of the energy (.): [IGa] for stereo matching with an absolute difference prior and [Ish] which generalizes the first method to handle arbitrary convex priors. Both of these methods require that labels are linearly ordered integers L = {, . . . , k} with k = |L|, which limits their applicability to few problems. Unsurprisingly, therefore, for many interesting instances of multilabel minimization, researchers resort to approximate labeling algorithms which produce provably good local solutions. One Graph Cut method in particular, called the Expansion Moves algorithm [BVZ], produces local minima within a known factor of the global minimum. Its practical advantages have been demonstrated in a range of vision problems on which it was shown to outperform popular multilabel minimization algorithms in terms of runtime and quality of obtained minima [SZS + ]. In our work, we use the Expansion Moves algorithm to

 term V condition c K F i (x i ) {s, t, i} F i () -F i () >  c(s, i) = F i () -F i () F i () c(i, t) =  F i () -F i () >  c(s, i) =  F i () c(i, t) = F i () -F i () F i,j (x i , x j ) {s, t, i, j} F i,j (, )-F i,j (, ) >  F i,j (, )-F i,j (, ) >  c(s, i) = F i,j (, ) -F i,j (, ) F i,j (, ) c(s, j) = F i,j (, ) -F i,j (, ) c(i, t) = c(j, t) =  F i,j (, )-F i,j (, ) >  F i,j (, )-F i,j (, ) >  c(s, j) = F i,j (, ) -F i,j (, ) F i,j (, ) c(i, t) = F i,j (, ) -F i,j (, ) c(s, i) = c(j, t) =  F i,j (, )-F i,j (, ) >  F i,j (, )-F i,j (, ) >  c(s, i) = F i,j (, ) -F i,j (, ) F i,j (, ) + F i,j (, ) -F i,j (, ) c(j, t) = F i,j (, ) -F i,j (, ) c(s, j) = c(i, t) =  F i,j (, )-F i,j (, ) >  F i,j (, )-F i,j (, ) >  c(i, t) = F i,j (, ) -F i,j (, ) F i,j (, ) c(j, t) = F i,j (, ) -F i,j (, ) c(s, i) = c(s, j) =  F i,j (, ) + F i,j (, ) -F i,j (, )-F i,j (, ) ≥  c(i, j) = F i,j (, ) + F i,j (, ) -F i,j (, )-F i,j (, ) c(j, i) =  Table 
.: Summary of graph constructions due to [KZ] for submodular binary functions of the form (.). For every term, the subgraph definition corresponding to indicated positivity conditions is given along with the constant K so that F(x) = c(S, T ) + K. The graph representing F is obtained by merging subgraphs for unary and pairwise terms and adding edge capacities and constants K together respectively.

minimize multilabel segmentation energies, therefore in the following we shall discuss its properties and present its algorithmic details.

The Expansion Moves algorithm takes its name from its search space exploration strategy, the expansion move, whereby an arbitrarily large number of variables can switch to a label α simultaneously. The algorithm finds a local minimum such that no expansion move, for any label α, yields a labeling with lower energy. This strong local minimality criterion is actually what bounds solutions found by the algorithm. By contrast, traditional energy minimization algorithms, like the Iterative Conditional Modes (ICM) algorithm [Bes] and the Simulated Annealing (SA) algorithm [KGV], use a much simpler search space exploration strategy, called the standard move, where only a single variable can change its label at a time. In the following, we formally define expansion and standard moves and compare their local minimality criteria.

A configuration can alternatively be seen as a partition of the set of variables, indexed by S, on which the energy (.) is defined. It is given by P = {P l | l ∈ L}, where P l = {i ∈ S | i = l} is the subset of variables that are assigned the label l in the configuration . Given a label α ∈ L, a move from an initial configuration , mirrored by the partition P , to a new configuration attainable from , mirrored by the partition  Given a configuration , denote the set of possible standard or expansion moves from by M . We say that ˆ is a local minimum of the energy E with respect to M if:

∀ ∈ M , E( ˆ ) ≤ E( ). (.)
It is easy to see that for standard moves |M | is linear in the number of variables, which makes it easy to find a local minimum with respect to these moves. However, the local minimality condition (.) when standard moves are allowed is a very weak one, because when at a local minimum with respect to these moves, one cannot find a configuration with lower energy by changing the label of a single variable. This explains the low quality results frequently generated by untuned ICM and SA algorithms and their sensitivity to initialization especially in minimizing high-dimensional nonconvex energies which have many such local minima. By contrast, when expansion moves are allowed, |M | is exponentially large since if the move involves m variables of a total of n, the number of expansion moves is n m . This makes the condition (.) so demanding and the number of such local minima so few that [BVZ] are able to prove that any configuration locally minimal with respect to expansion moves is within a known factor of the global minimum. The strong local minimality property also makes expansion moves less sensible to the initial labeling. [BVZ] show that Given a configuration , there are, as we have seen, an exponential number of expansion moves from , therefore a direct search for the optimal expansion move is infeasible. This is where binary Graph Cut techniques discussed in Section .. come  into play. [BVZ] show that computing the optimal expansion on a label can be specified as a binary optimization problem, in effect, transforming the multilabling problem to a sequence of binary subproblems each formulated as a Maximum flow problem on an appropriate graph. In the following, using the results of [KZ] we show how to formulate an expansion move as a binary labeling problem and find the optimal move, that is the move that gives the largest decrease in energy, by computing a minimum cut. We also derive the multilabel submodularity condition mentioned earlier (.) that must be met by the pairwise terms of the energy (.) so that the method applies.

E( ˆ ) ≤ cE( * ) where c =  × max V i,j (a, b)/min V i,j (a, b), ∀i, j ∈ S,
Consider some initial configuration and a label we wish to expand on α ∈ L. Recall that after an expansion move on α, a variable i ∈ S either keeps its current label i , or switches to α and that variables that are already labeled α naturally keep their label. We can define such a label attribution as a function of binary variables

f i (x i ) = i , x i =  α, x i =  (.)
and write the multilabel energy function (.) as follows:

E(x  , . . . , x n ) = t i∈S D i (f i (x i )) + i,j∈S, j∈N i V i,j (f i (x i ), f j (x j )). (.)
As before, a graph with two terminal nodes can be constructed to represent the energy (.) with the following conventions. Given a cut (S, T ), if i ∈ S, then x i = , therefore f i () = i and the label keeps its current label, otherwise if i ∈ T , then x i = , therefore f i () = α and the label switches to α. The graph can then be constructed by the method of [KZ] presented in Section ... and summarized in Table ..

The original work which proposed the Expansion Moves algorithm [BVZ] relied on a graph construction requiring metric pairwise terms, that is ∀α, β, γ ∈ L:

V (α, β) =  ⇔ α = β, V (α, β) = V (β, α) ≥ , V (β, γ) ≤ V (β, α) + V (α, γ). (.)
In particular, symmetry is assumed by [BVZ] only to simplify graph constructions. [KZ] subsequently showed that metricity is an accidental property of elaborate graph constructions the original method employed. They conclude that the energy function (.) is graph representable if and only if each pairwise term satisfies the inequality:

V i,j (f i (), f j ()) + V i,j (f i (), f j ()) ≤ V i,j (f i (), f j ()) + V i,j (f i (), f j ()) (.) that is V i,j (β, γ) + V i,j (α, α) ≤ V i,j (β, α) + V i,j (α, γ) (.)
where β = i and γ = j . Setting V i,j (α, α) = , we obtain the triangle inequality, as required by [BVZ].



Algorithm ..: Expansion Moves energy minimization input : set of labels L, energy function of label configurations E( ) on a set of variables indexed by S, maximum cycles maxCycles output: label configuration

// Initialize configuration  l init ← Random (l ∈ L)  foreach i ∈ S  i ← l init // Perform energy minimization  cycles ←   repeat  success ← f alse  foreach α ∈ L  ˆ = arg min E( ), ∀ ∈ M  if E( ˆ ) < E( )  ← ˆ  success ← true  cycles ← cycles +   until success = f alse or cycles = maxCycles
Having formulated individual expansion moves as binary subproblems, we now return our attention to the Expansion Moves algorithm. Its listing is given in Algorithm ... We shall call the execution of lines - a cycle and an execution of lines - an iteration. Starting from a constant labeling by an arbitrarily chosen label (lines -), the algorithm cycles through labels in a fixed or a random order. In each iteration, the optimal expansion on a label is computed with the method outlined in the above (line ), and if the energy of the resulting configuration is strictly inferior (line ), the current configuration is updated (line ). A cycle is successful if a strictly better configuration is found in one of its iterations. The algorithm stops after the first unsuccessful cycle or if the required maximum number of cycles are executed.

The Expansion Moves algorithm .. is guaranteed to terminate since every iteration retains the configuration that decreases the value of the energy. [BVZ] mention that for an energy formulation which is independent of the number of variables, the algorithm performs O(V ) cycles, where V is the set of variables. Therefore, the running time of the algorithm is O(V L). However, the authors report that in practice the algorithm stops after a few cycles. This observation is confirmed by our experimental setup as well. A precise description of the algorithm's total running time requires taking into account the running time of the minimum cut computation (line ) for which a maximum flow algorithm is applied. For example, using a straightforward implementation of the Push-relabel method with the FIFO selection rule [GT] which runs in O(V  ) time, the total running time of the Expansion Moves algorithm becomes O(V  L).



. Conclusions

In this chapter, we have laid the theoretical foundation of our approach that we shall elaborate in subsequent chapters formulating image segmentation as a Bayesian labeling problem which is solved by minimizing a posterior probability energy, described in terms of likelihood and Gibbs distribution energies, to obtain the optimal segmentation with respect to user-provided evidence, constraints and a problem-specific prior.

We have addressed energy minimization methods discussing Graph Gut algorithms in exact binary and approximate multilabel optimization as well as the underlying maximum flow techniques. We choose the Expansion Moves algorithm for minimizing multilabel segmentation energies which finds a local solution within a known factor of the global minimum. Having a bound on local minima helps us to explain whether an inadequate solution is due to a bad energy definition or a poorly performing minimizer.

Even though good local minima correspond to good quality solutions, the quality of a labeling result in segmentation is usually assessed with respect to a ground-truth labeling. As reported, the global minimum, when known, is not more accurate than the ground-truth, and the energy of the ground-truth is frequently higher than that of the global minimum, therefore the focus is not on the optimization of energies, but on finding more expressive energy models. This is the objective of the following chapter.



Chapter 

A structural prior model for Graph Cut segmentation

. Introduction

M

ost successful segmentation methods usually incorporate some image or domain-specific prior information. Medical image segmentation methods in particular should ideally be able to incorporate a wide range of anatomical and physiological prior information to ensure the consistency of the result with respect to anatomical properties and variabilities thereof. From the standpoint of Graph Cut optimization of energies of label configurations, prior information can be regarded as a set of constraints on labels defined by a reference model. In Section ., we have encountered unary and pairwise constraints stemming from the energy potentials of a prior distribution expressed as a Markov random field (MRF) of label configurations with respect to the image graph. Such prior information is combined with data-derived likelihood and user constraints in an energy minimization framework to infer the optimal segmentation according to these criteria.

As mentioned in Section .., in image segmentation, pairwise prior interaction potentials have to encourage a piecewise-constant discontinuity-preserving labeling. The most widely used such prior, called the Potts model (.), presumes the presence of a single type of boundary in the image. With the assumption of equal violation of piecewise constance, it penalizes attributions of pairs of unequal labels to neighboring variables equally while incurring no penalty for the identical labeling of such pairs. This assumption obviously does not hold in multiobject segmentation where multiple boundaries exist between the objects.

In this chapter, we investigate extensions of the Potts prior model to account for multiple boundaries in multiobject segmentation which incur multiple levels of penalization and furthermore allow the minimization of resulting energies by the Expansion Moves Graph Cut algorithm discussed in Section ... We first consider a simple structure adjacency prior by adding a second level of penalization to the Potts model [JKD + ]. Due to the submodularity requirement, Expansion Moves optimization can  be applied only to a restricted form of this prior where the penalty of nonadjacency is at most twice that of adjacency. We outline a submodular extension of this prior for simple linear label layout models to penalization ratios greater than two by the introduction of auxiliary labels into the model. The main contribution of this chapter, motivated by a special case of the aforementioned submodular extension of the adjacency prior, is presented last. The vicinity prior, as we term it, is an intrinsically submodular prior model which incurs multiple levels of penalization, |L| at most, and captures the spatial configuration of structures in multiobject segmentation. It is defined as shortest-path pairwise constraints on a graph model of interobject adjacency relations [KVDP].

The following section reviews relevant literature and suggests a taxonomy of prior introduction approaches in Graph Cut segmentation.

. Prior information in Graph Cut segmentation

Many forms of prior information have been introduced in Graph Cut segmentation for different applications. We broadly classify the most relevant of such approaches in three categories according to the nature of introduced prior.

a) Location prior, formulated as spatially-varying prior probability maps of individual structure locations in the image space. Such prior information is usually encoded by the unary constraints of the prior energy (.) to indicate the label preference of individual variables [STAG, KRBT, LPPS] and sometimes via its pairwise constraints to favor the placement of the segmentation boundary near the object's contour [FZ]. In simplest cases, probabilistic maps of location are defined as a distance transform of a shape template [FZ, KRBT]. In medical applications, probabilistic atlases of tissue location are usually constructed by averaging label estimates for each pixel over a ground-truth dataset registered to some anatomical space [STAG, LPPS]. Naturally, such probabilistic maps have to be learned in advance and registered to an image before it is segmented. Registration cannot trivially handle topological and geometrical changes, due to pathologies for example, and is often computationally intensive.

b) Shape prior, where explicit or implicit shape constraints are used to match the segmented object with predefined or learned shapes. [VM] define an affine-invariant shape similarity measure and introduce it in a sequential multiobject Graph Cut segmentation framework. [MRT] first learn a statistical shape model from a training set via kernel principal component analysis, then use the projection of an initial segmentation in the learned shape space as the shape prior for subsequent segmentations in an iterative process. In their Graph Cut framework for the optimization of hypersurface functionals, [KB] show how to set unary constraints to favor cuts that maximize the flux of distance transform gradients for a shape's contour or its medial axis. Shape-prior segmentation approaches entail shape location initialization and alignment mechanisms, and often require model learning.

c) Spatial configuration prior, the use of which usually arises in the context of multiobject segmentation, where usually pairwise relations define interobject geometric and spatial relationships of the type "object a is above object b and next to object c". The  use of such prior information has appeared in the context of part-based object detection [WS] and scene layout estimation and segmentation where spatial relationships are either predefined on a limited number of classes [LVS] or learned from annotated training data [WS, LLB + ]. [DB] segment multiple or compound objects encountered in relative positions of containment and exclusion, also introducing graph-edge constraints on preferred distances between the boundaries of regions. Some useful configurations however cannot be represented by definitions therein, such as the inclusion of an adjacent pair of objects in a third.

There are a number of spatial configuration prior formulations closely related to our work which have appeared outside Graph Cut segmentation literature. In a theoretical framework similar to ours, although following a different numerical optimization scheme, [SSLS + , CCB + ] use a simple three-level piecewise-constant prior model for tissue classification in brain MR images. According to a brain tissue model, it encourages the identical classification of pairs of neighboring pixels, allowing different class adjacency patterns which are consistent with the model while penalizing adjacencies which are not allowed to occur. [GRC + ] proposes a relative location prior represented by probabilistic maps of interclass offset preferences for pairs of classes. Such probabilistic maps are learned in advance from training data and are used to define a local relative location feature which is introduced into the segmentation energy via unary constraints. In [SBZ + ], pairwise geometric relationships are defined by an anatomical network on easily detectable body landmarks and used in the organ detection phase of a sequential segmentation framework for thoracic and abdominal organs in whole-body CT scans. [OYH + ] addresses the same problem as [SBZ + ] with an automatic two-stage method whereby once relatively stable organs are identified and segmented, organ interrelations are used to proceed with the segmentation of variable organs by fitting previously learned statistical organ atlases and shape models. Organ interrelations are defined according to canonical correlation analysis and represented by a directed graph. Intuitively, these relations represent the degree of influence of an organ's surface upon another. A generic, graph-based scene model representing spatial information explicitly is also employed by [FAB] in the problem of segmenting multiple subcortical structures in brain MR images. Fuzzy representations of distance, adjacency and spatial directionality relations are introduced to account for model imprecisions and variabilities. The segmentation algorithm, however, is sequential and requires error detection and segmentation order rearrangement strategies.

. An adjacency prior

Let us first consider the image in Figure .a. It shows our synthetic phantom featuring  structures of different size, topology and adjacency patterns, labeled  to  as indicated. More details can be found in Section . where we present evaluation results on a dataset generated from the phantom by introducing variable Gaussian noise. In the remaining of this chapter, we shall use this phantom as a reference model to illustrate concepts as we develop definitions of prior models. In all our formulations, we shall rely on pairwise constraints, expressed by pairwise energy terms V i,j ( i , j ) in (.), to encode prior information on pairwise interactions between labels assigned to pairs of neighboring variables encouraging the spatial consistency of labeling with respect to a problem-specific reference model. Assuming a piecewise-constant prior model, the simplest manifestation of such consistency is that the spatial layout of image regions labeled differently should reflect the local spatial constraints defined by the prior model. In particular, given a pair of labels for which no close-range spatial interaction is described by the model, no pair of image regions labeled with this pair of labels should appear next to each other in the labeled image.



In a segmentation task of an image generated from the phantom in Figure .a, the use of the Potts prior model incurs equal penalty for all but pairs of identical labels. For example, given the adjacent and the nonadjacent pairs {, } and {, } respectively, the Potts model implies that attributions of the corresponding pairs of nonidentical labels to pairs of neighboring variables violate the spatial consistency of labeling equally, since ∀i, j ∈ S, V i,j (, ) = V i,j (, ) = .

In similar spirit to [SSLS + , CCB + ], we extend the Potts prior by an additional constraint level to penalize label assignments which violate layout consistency with respect to structure adjacency relations as described by the reference model. Such a formulation is able to distinguish pairings of labels corresponding to adjacent and nonadjacent structures in the model and penalizes them differently while preserving the piecewise constance property of the prior since V (a, a) = , ∀a ∈ L. For the example just illustrated, the effect would be ∀i, j ∈ S, V i,j (, ) > V i,j (, ).

Let R be the set of symmetric binary relations on pairs of distinct labels, R = {r | a r b, a, b ∈ L, a b}. Such binary relations can represent many types of relationships between pairs of objects, such as adjacency, distance, area of common surface, etc. At present, we let relations r ∈ R represent structure adjacency. A structure a is considered to be adjacent to another b, denoted as a r b, if they share a boundary (in D, a surface). We define the adjacency prior as follows:

 V i,j i , j =            , i = j (.a) γ sof t , ∃r ∈ R, i r j (.b) γ hard , r ∈ R, i r j (.c)
where γ sof t and γ hard correspond to adjacency and nonadjacency penalties respectively and γ hard > γ sof t . Ideally, γ hard → ∞, however, some energy terms would then violate the submodularity condition (.), namely those falling on the left-hand side of the inequality. One possible remedy to this problem is the "truncation" of nonsubmodular terms [RKKB, WS, DB], which is the modification of energy terms on either side of an unsatisfied submodularity condition so that the inequality holds. As pointed out by [RKKB], the Expansion Moves algorithm with truncated energy terms is a valid energy minimization technique for arbitrary functions and is guaranteed to converge. However, individual algorithm steps are no longer guaranteed to find the optimal labeling with respect to expansion moves and the final labeling does not necessarily have the strong local minimality property. They point out that the technique gives reasonable results as long as the number of nonsubmodular terms remains small. Another solution is using an efficient arbitrary energy minimization algorithm, like the TRW-S algorithm [Kol], which however does not have a bound on local minima.

The adjacency prior as defined by (.) is submodular when γ hard ≤  γ sof t , with the largest penalization ratio at equality. Clearly, this prior formulation is essentially limited to definitions where the penalty of nonadjacency is at most twice as much as the penalty of adjacency so that the prior is submodular and the energy (.) can be minimized via the Expansion Moves algorithm. In other words, it is impossible to penalize nonadjacencies by some large constant without penalizing adjacencies half as severely. In the following section, we outline a submodular extension of this prior for simple linear label layout models to penalization ratios greater than two by the introduction of auxiliary labels into the model.

.. Arbitrary penalization submodular formulations

Consider the simple three-label model illustrated in Figure .a with the label set L = {l, m, r}. The adjacency prior for this model is defined according to (.) as follows:

∀a ∈ L, V (a, a) = , V (l, r) = V (l, r) = γ hard , V (l, m) = V (m, l) = V (m, r) = V (r, m) = γ sof t . (.)
Let us suppose that γ hard >  γ sof t . Consider the set of inequalities of the submodularity condition for all permutations of labels in L. Because of symmetric pairwise definitions in (.), there are three such inequalities where each pair of distinct labels in {l, m, r} falls on the left-hand side of an inequality exactly once. The inequality which does not hold for γ hard >  γ sof t , is the following: To satisfy (.), we can either increase one of V (m, r) or V (l, m), or decrease V (m, m), as pointed out by [RKKB], or use an arbitrary energy minimizer, for example [Kol].

V (m, m) + V (l, r) ≤ V (m, r) + V (l, m). (.)  r l m r l m γ h γ s γ s r l m m γ h γ h - γ s γ s γ s γh -γs γ h -γ s (a) (b) (c) (d)
Without modifying the prior definition, in particular the ratio γ hard /γ sof t , we write the violated inequality (.) as two inequalities that can be satisfied on the set of labels L ∪ m . We call m the shadow label of m representing the same image structure. Consider the graphical representation of the triangle inequality in Figure .b, where triangle vertices correspond to labels l, m and r and lengths of its sides to prior penalties such that ∀a, b ∈ L, ab = V (a, b). Note that the illustration in .b corresponds to γ hard <  γ sof t and that the geometrical representation for γ hard =  γ sof t corresponds to collinear points. In the following, we shall denote γ hard and γ sof t shortly as γ h and γ s .

Our intuition is pulling the vertices l and m of the triangle lrm apart so that the side lr is longer, in particular lr > lm + rm . As illustrated in Figure .d, this splits the vertex m and results in two triangles, lrm and lrm , according to which we rewrite the violated inequality (.) as follows using line-segment notation:

lr ≤ lm + mr (.a) lr ≤ lm + m r (.b)
where lr = γ h and lm = m r = γ s and γ h >  γ s . We shall derive bounds on unknowns mr, lm and mm . From (.a), mr ≥ lrlm, that is, mr ≥ γ hγ s . Also from lrm, mr ≤ lr + lm, that is, mr ≤ γ h + γ s . Thus:

γ h -γ s ≤ mr ≤ γ h + γ s . (.)
Similarly, from (.b) and lrm we find:

γ h -γ s ≤ lm ≤ γ h + γ s . (.)
We now consider mm . From lmm we have, mm ≤ lm + lm , mm ≤ γ s + γ hγ s = γ h assuming lower bound on lm . We also have mm ≥ lmlm, that is, mm ≥ γ h - γ s under the same assumption. Hence: Since m and m represent the same structure, the upper bound γ h , equal to lm which is the highest penalty in the model, is too severe. Therefore, we shall assume lower bound on mm . This completes the construction of the extended adjacency prior that accounts for arbitrary ratios γ hard /γ sof t for the simple three-structure model in

γ h - γ s ≤ mm ≤ γ h . (.)  r l m x l r x l x (a) (b) (c) (d) r l m m γ h  γh - γs γ s γ s γh -γs γ h -γ s γ h - γ s  γ h -γ s γ h -γ s γh -γs r x γ h  γ h - γ s  γ h - γ s γ h - γ s
V (r, x) = V (x, r) = γ sof t , ∀a ∈ {l, m}, V (a, x) = V (x,



From lrr we have, lr ≤ lr + rr , that is, lr

≤ γ h + γ h - γ s =  γ h - γ s .
From the same triangle, we also have lr ≥ lrrr , that is, lr ≥ γ hγ h +  γ s =  γ s . Therefore:

 γ s ≤ lr ≤  γ h - γ s . (.)
Similarly, from xmm we find that:

 γ s ≤ xm ≤  γ h - γ s . (.)
Upper bounds of (.) and (.) can also be obtained from lr m and xrm respectively. Let us derive bounds on lx and mr . From lxr we have, lx

≤ lr + r x, that is, lx ≤  γ h - γ s + γ s =  γ h -γ a assuming upper bound on lr . We also have, lx ≥ lr -r x, which is, lx ≥  γ h - γ s -γ s =  γ h - γ s . Hence:  γ h - γ s ≤ lx ≤  γ h -γ a . (.)
From mrr we have, mr ≤ mr + rr , which is, mr

≤ γ h -γ s + γ h - γ s =  γ h - γ a . We also have, mr ≥ mr -rr , that is, mr ≥ γ h -γ s -γ h +  γ s = γ s . Hence: γ s ≤ mr ≤  γ h - γ a (.)
which completes the construction of the extended adjacency prior accounting for arbitrary ratios γ hard /γ sof t for the four-structure linear model in In order to understand the practical merit of extended versus basic adjacency prior models, we conduct an experiment of brain tissue segmentation on a single simulated D MR image from the second-generation BrainWeb simulated MRI dataset. The dataset is described in Section ... To test the four-label extended adjacency prior defined above, we formulate Graph Cut segmentation on the label set L = {"BKG", "CSF", "GM", "W M"} representing background, cerebrospinal fluid, gray and white matters respectively. The background comprises all nonbrain tissues in addition to air. The spatial layout of structures in L largely corresponds to that of the four-label model illustrated in Figure .a. We estimate likelihood distributions for the last three labels from user-supplied "seeds" as described in Section .., using distributions identical to those of the original label for shadow labels and a uniform distribution for the background label. We measure the quality of segmentation with the overall Dice similarly metric (DSM) (.) which indicates the agreement between the segmentation and the ground-truth as a mean weighted measure of individual-label DSM measures.

Figure . gives results of qualitative and quantitative evaluation using Potts, basic adjacency prior and its extension for three ratios γ hard /γ sof t equal to ,  and . Presented cross-sections and DSM figures suggest that there are no real advantages of using the extended adjacency prior compared to the basic adjacency prior on this instance of segmentation described by the four-label model. Furthermore, segmentation quality deteriorates with increasing gap between the two penalty levels and false positives of the label class "CSF" begin to appear, which however obey adjacency constraints. It is interesting to note that for γ hard /γ sof t = , the extended adjacency prior gives a definition that accounts for relative structure distance for the linearly ordered fourlabel model:



V (l, m) = V (m, r) = V (r, x) = γ sof t V (l, r) = V (m, x) =  γ sof t , V (l, x) =  γ sof t .
(.)

Given the complexity of extended adjacency prior constructions, we seek a simpler specification of the idea of stepwise penalization which would furthermore be easier to define for more complex nonlinear models, counting several tens of labels.

. A shortest-path formulation: the vicinity prior

The adjacency prior discussed in Section . ensures the spatial consistency of labeling by penalizing pairings of labels corresponding to adjacent and nonadjacent pairs of structures differently (.). The formulation however makes the assumption that pairings corresponding to nonadjacent structures are equally unfit, independently of the mutual distance between structures. Higher penalization for farther away pairs of structures while maintaining the submodularity property of the prior is desired.



We have observed that a special case of the extended adjacency prior, corresponding to γ hard =  γ sof t , accounts for relative structure distance for the simple linearly ordered four-label model .a incurring multiple levels of penalization (.). In this section, we present a different formulation of the idea of variable penalization according to stepwise relative distance which can be furthermore specified easily for much more complex nonlinear models defined on an arbitrary number of labels. Consider the phantom in Figure .a again, and observe the layout of structures "", "", "" and "". Notice how the structure labeled "" is nested in the structure "" which, in turn, is adjacent to the structure "". Due to this adjacency pattern, it makes sense to encourage the assignment of the pair of labels {, } while penalizing the attribution of the pair {, } to neighboring variables in a segmentation task of an image generated from the phantom. Hence, ∀i, j ∈ S, V i,j (, ) > V i,j (, ). Furthermore, since the structure "" is farther still, V i,j (, ) > V i,j (, ) > V i,j (, ) should hold. By comparison, the adjacency prior penalizes nonadjacent pairs {, } and {, } equally implying that they are equally unfit.

Recall that R is the set of symmetric binary relations on pairs of distinct labels representing adjacency relations of corresponding structures. Assume furthermore that each relation r ∈ R has a corresponding characteristic function r which gives a quantitative measure for the relationship r represents, r :

L × L \ {(a, a)} → R + .
The set of relations R can be equivalently represented by a weighted undirected graph A = L, W without self-loops  on labels L, where the set of edges is given by W = {{a, b} | ∃r ∈ R, a r b, a b}. Edge weights are defined by w ({a, b}) = r (a, b), such that

w ({a, b}) = ∞, r ∈ R, a r b.
Define a path of length  k from a vertex l  to another l k in A by a sequence π = (l  , l  , . . . , l k ). We say that the vertex l k is reachable from the vertex l  via the path π and denote it as l  π ; l k . Define the weight of a path w(π) as the sum of the weights of its constituent edges w(π) = k- i= w ({l i , l i+ }). The shortest-path weight from a vertex a to another b is defined by:

(a, b) = min w(π), ∀π, a π ; b (.a) ∞, otherwise (.b)
where (.b) corresponds to the case where there is no path π through which the vertex b can be reached from the vertex a. Note that the shortest path from a to b is defined as any path π with weight w(π) = (a, b) and that the shortest-path weight from a vertex to itself is (a, a) = .

Given the graph A, we define the pairwise prior term in (.) as:

V i,j i , j = (a, b) , i = a, j = b. (.)
There are several algorithms to calculate shortest-path weights the choice of which depends on the underlying graph and the application [CLRS]. For weighted graphs with positive edge weights, shortest paths from a given vertex to all others can be found  A self-loop is an edge from a vertex to itself.  The length of a path is equal to the number of edges on it.

 via Dijkstra's algorithm with a running time of O(L  ). If the graph is sparse and every vertex can be reached from the starting vertex, the running time becomes O(W log L) when the algorithm's priority queue is implemented with a binary heap. For graphs in which each edge has unit weight, shortest paths from a given vertex to all others can be found via the breadth-first search (BFS) algorithm with a running time of O(L + W ).

In our formulation, relations r ∈ R represent structure adjacency. For two such structures a and b, the adjacency relation's characteristic function is given by:

r (a, b) = , ∃r ∈ R, a r b (.a) ∞, otherwise. (.b)
This results in a connected undirected graph A with unit-weight edges between adjacent vertices, where the degree of each vertex l ∈ L is equal to the number of structures adjacent to the structure associated with it. Shortest paths from a vertex to all others, conveniently calculated by an application of BFS to the graph A, define the vicinity of the associated structure with respect to all others, which is introduced in the energy function of segmentation (.) via pairwise terms (.). Clearly, the graph A constitutes a prior model of explicit spatial configuration of objects. To illustrate, the graph prior model of our synthetic phantom in Figure .a is given in Figure .b. Note that the graph A is not limited to unit-weight edges. It can encode any quantity that accumulates linearly along a path that we would want to minimize.

The vicinity prior as defined by (.) is intrinsically submodular. Therefore, the Expansion Moves algorithm can by applied to minimize segmentation energies using this prior without resorting to model truncation, unlike [RKKB, WS, DB]. The following theorem establishes the submodularity property of the vicinity prior.

Theorem .. (Submodularity of vicinity prior). Given the undirected weighted graph A = L, W and the prior definition V i,j (a, b) = (a, b), the inequality V i,j (a, a) + V i,j (b, c) ≤ V i,j (a, c) + V i,j (b, a) holds ∀i, j ∈ S for any a, b, c ∈ L. (.)

First assume that the graph A = L, W is connected, that is ∀a, b ∈ L, ∃π, a π ; b. Let π ba be the shortest path with weight w(π ba ) = (b, a), similarly, let π ac be the shortest path with weight w(π ac ) = (a, c). The shortest path π bc either passes by a, that is, π bac = (b, . . . , a, . . . , c), in which case its weight is w(π bac ) = w(π ba ) + w(π ac ), hence (b, c) = (b, a) + (a, c), or it passes by a vertex q ∈ L not on π bac through a different path π bqc = (b, . . . , q, . . . , c) with weight w(π bqc ) < w(π ba ) + w(π ac ), therefore (b, c) < (b, a) + (a, c).

If the graph A is unconnected, infinite shortest-path weight might appear on the right-hand side of (.) when either or both b and c are unreachable from a, and the inequality holds. However, infinite shortest-path weight cannot appear only on the left-hand side of (.), because if b ; c with w(π bac ) = w(π ba ) + w(π ac ) thus completing the proof.



. Discussion and conclusions

In this chapter, we have investigated discontinuity-preserving piecewise-constant submodular prior formulations for multiobject segmentation which, unlike the Potts prior model, incur multiple levels of penalization accounting for multiple types of boundaries. Starting from the Potts model, we have devised a simple structure adjacency prior by introducing a second level of penalization for nonadjacent structures. The submodularity requirement restricts this model to definitions where the penalty of nonadjacency is at most twice that of adjacency. We have outlined submodular extensions of this prior for linear spatial layout models to penalization ratios greater than two by the introduction of auxiliary labels into the model. A special case of the extended formulation gives rise to a variable penalization scheme accounting for relative structure distance on the linear spatial layout models. We proposed a graph-based formulation of this scheme which can be specified easily for complex nonlinear models defined on an arbitrary number of labels. The vicinity prior, as we name it, is defined as shortest-path pairwise constraints on a graph model of interobject adjacency relations. It is intrinsically submodular and incurs |L| levels of penalization at most, furthermore capturing the spatial configuration of structures in multiobject segmentation. In Chapter , we give the results of quantitative and qualitative evaluation of Potts, adjacency and vicinity prior-based segmentation on synthetic, simulated and real medical images.

While theoretically interesting, we did not observe any practical advantages of the extended compared to the basic adjacency prior model in the problem of brain tissue segmentation. One explanation is that the Expansion Moves algorithm has a worse bound on local minima (.) for the extended model compared to that of the basic model, in addition to the labeling problem being defined in a higher-dimensional space in this case. At any rate, further empirical evidence is necessary to reach a conclusion on the practical merit of the extended adjacency prior model.

Compared to shape prior-based approaches, structural prior models capturing the spatial configuration of objects are quite robust to shape deformations since relative positions of objects remain largely stable. Graph models describing such priors are fairly easy to define at an arbitrary level of detail from simple specifications of spatial relationships on pairs of objects, and can be generated from existing anatomical models for medical applications in particular. Interobject adjacency patterns can be learned from training data or several model variants thus capturing model variabilities in which case continuous graph-edge weighs may be used, for example, to represent degrees of adjacency. Pairwise spatial configuration prior formulations learned from data are frequently nonsubmodular, therefore authors resort to model truncation to apply Graph Cut minimization [WS] or use arbitrary-energy minimizers [GRC + , LLB + ] thus sacrificing the accuracy of the model in the former and the quality of local minima in the latter case. Compared to such learning-based approaches, the vicinity prior formulation has the advantage that of all pairings on L × L only those corresponding to structures that exhibit adjacency in the dataset need to be explicitly specified or learned; on the remaining pairs, the prior is defined as shortest-path lengths on the prior graph.



Chapter 

Reducing Graph Cut complexity . Introduction W e have seen in Section . that Graph Cut algorithms usually employ rectangular grid graphs as their data structure. In segmentation applications on high-resolution D medical images where several tens of millions of voxels are common, such graphs have to be defined on a number of nodes equal to the number of image voxels and, when -connected voxel neighborhoods are used, a number of edges at least thrice the number of nodes for neighborhood links plus twice the number of nodes for terminal links. With data structure size touching on the limit or exceeding the storage capacity of computer memory along with supralinear time complexity of best performing maximum flow algorithms [BK], numerical optimization becomes impractical.

Several approaches have been proposed to address the complexity of Graph Cut and other graph-based image segmentation approaches on high-resolution input. They can be broadly classified into two categories; algorithm parallelization and image simplification. Parallel versions of some maximum flow algorithms of the Push-relabel variant discussed in Section .. have been devised for shared-memory multiprocessor architectures [DB] and graphics processing units (GPU) [VN]. While these methods attain good speedups, they offer no solutions for the memory footprint problem since they operate on the entire image domain. Image simplification approaches address this problem directly observing that the full-resolution pixel representation is often redundant because objects usually comprise many similar pixels that could be grouped. In [LSGX], Graph Cut segmentation is first applied to a low-resolution version of the image and the solution is propagated to the full resolution where Graph Cut segmentation is applied in a narrow band surrounding the projected foreground-background boundary. In order to accelerate feedback in their interactive Graph Cut segmentation scheme, [LSTS] exploit the fact that Graph Cut algorithms lend themselves equally well to irregular nongrid graphs and oversegment the image into small regions via the Watershed algorithm [VS]. Graph Cut segmentation is then formulated according to the watershed regions graph. Starting from a grid partition, [CA] cluster image pixels  by an iterative k-means algorithm augmented by color similarity and shape compactness constraints. The image is then represented by a cluster graph which is partitioned via the Normalized Cuts approach [SM] based on color information only yielding a coarse oversegmentation of the image. Recent conditional random field (CRF) image segmentation approaches such as [GRC + , FVS, LLB + ] rely on image clustering not only to reduce memory footprint and accelerate runtime, but also to collect image features from clusters and their neighborhood.

In order to reduce memory footprint and accelerate Graph Cut optimization, we cluster image voxels prior to segmentation by an intensity and gradient-adaptive centroidal Voronoi tessellation (CVT). Image partitions produced by this clustering method strike a good balance between the conflicting criteria of cluster compactness and structure boundary adherence furthermore allowing to control the tradeoff between them and helping towards the fine-grained placement of segmentation boundaries. We formulate Graph Cut segmentation according to the graph of CVT clusters so that the resulting labeling energy, defined after (.), is independent of the clustering resolution [KVDP]. Compared to a direct application of segmentation on voxels, the clustering step improves overall runtime and memory footprint of the segmentation process up to an order of magnitude virtually without compromising the segmentation quality.

. Related work

Several image tessellation algorithms have been developed by vision and pattern analysis communities for image classification and segmentation applications. Of special interest is a class of image partition algorithms which produce a dense oversegmentation of compact segments of relatively uniform size and shape which adapt to local intensity edges in the image. These are often called superpixels, a term introduced by [RM] who proposed a Normalized Cuts-based algorithm that recursively bipartitions an image minimizing a cut cost function. Superpixels algorithms with better runtime performance and comparable oversegmentation quality have been introduced since, such as the TurboPixels algorithm of [LSK + ] and the Graph Cut superpixels algorithm of [VBM]. The oversegmentation produced by these algorithms achieves a better balance between the conflicting goals of cluster compactness and structure boundary adherence than some well-known image partition algorithms which produce segments of highly variable shape and size, like Watershed [VS] and Mean-shift [CM] algorithms and the oversegmentation algorithm of [FH]. Unsurprisingly, a large segment of irregular shape is more likely to span more than one object, especially in the absence of boundary cues in the image due to insufficient contrast or the presence of shadows. A compact, regular shape, on the other hand, is less likely to cross object boundaries unless they are too wiggly. Moreover, if such a segment does encompass more than one object, the overall error rate can be controlled if its size is not allowed to grow arbitrarily. Cluster compactness, however, comes at the expense of boundary adherence, therefore it is desirable to be able to control the tradeoff between the two criteria by an appropriate parameterization of the algorithm. Later in this chapter, we shall compare  the CVT clustering algorithm we employ with the aforementioned TurboPixels and Graph Cut superpixels algorithms, therefore we review them here before presenting our clustering method in the following section.

The TurboPixels algorithm [LSK + ] is a geometric flow algorithm cast in a banded level-set curve evolution that iteratively dilates an initial set of uniformly distributed contour seeds, each corresponding to one superpixel. Superpixel compactness is maximized by a term in the curve evolution equation which produces constant motion in the outward normal direction in image regions of uniform intensity maximizing the rate of area growth while maintaining the minimum possible isoperimetric ratio. The smoothness of superpixel boundaries and their alignment with local image edges are ensured, respectively, by curvature and gradient-based curve attraction terms in curve evolution equation velocity coefficients. To prevent expanding seeds from overlapping, background skeletonization is used to stop boundary evolution when two distinct dilating seeds are about to collide. We note that upon termination, the algorithm produces more superpixels than the initial number of seeds, since any remaining unassigned large connected regions are treated as superpixels.

Similar to our image segmentation theoretical framework presented in Section ., [VBM] formulate superpixel partition as a pixel labeling problem and optimize with the multilabel Expansion Moves Graph Cut algorithm [BVZ]. The image is initially covered with overlapping square patches of the maximum allowed superpixel size, each corresponding to a label, and the optimal "stitching" of these patches is sought, in other words, the optimal assignment of each pixel to only one such patch. In order to align patch boundaries with intensity edges, the stitching cost at image intensity edges is made cheaper by using a weighted Potts pairwise model, which, for the assignment of pairs of neighboring pixels to different patches, incurs a penalty inversely proportional to their intensity difference. The explicit energy formulation allows for different unary constraints to be specified in order to encourage different properties of superpixels. A formulation which was experimentally found to perform best specifies an absolute intensity difference-based unary term that encourages intensity homogeneity inside superpixels at the expense of slightly less regular superpixel sizes and shapes. We note that, instead of the number of superpixels, the maximum superpixel size is considered to be a more natural parameter for the algorithm which could be chosen by the user as appropriate for an application. The number of patches is deduced subsequently by spreading square patches on the image at certain intervals.

. Image clustering by centroidal Voronoi tessellation

Given a grayscale image I and n sites c i ∈ I such that  ≤ i ≤ n, a Voronoi tessellation partitions the image I into n disjoint clusters C i associated with each site c i , where every cluster C i consists of the set of all voxels in the image which, according to a distance measure d, are closer to the associated site than to all other sites:

C i = {v ∈ I | d(v, c i ) < d(v, c j );  ≤ j ≤ n, j i}. (.)
Unless otherwise mentioned, d is usually the Euclidean distance.



A Centroidal Voronoi Tessellation (CVT) [DFG] has the additional property that each site c i corresponds to the mass centroid of the associated cluster and is defined by:

c i = v∈C i ρ(v) v v∈C i ρ(v) (.)
where ρ(v) is a density function defined by (.). Moreover, a CVT minimizes the following energy function:

F(v; c i ) = n i=         v∈C i ρ(v) v -c i          . (.)
We extend the cluster geometry compactness property expressed by (.) by adding an intensity-space norm as well:

E(v; c i ) = n i=         v∈C i ρ(v) ( v -c i  + α I v -I i  )         (.)
where I v is the gray-level of the voxel v, I i is the gray-level of the cluster C i defined as the mean intensity of its voxels, and α is a positive scalar allowing to control the intensity adaptivity of the CVT partition. Intuitively, minimizing (.) now corresponds to maximizing cluster compactness in terms of both geometry and gray-level. In (.) and (.), we define ρ(v) as a linear function of intensity-gradient magnitude at voxel v: A number of algorithms exist to construct a CVT minimizing (.), including Lloyd's relaxation method for k-means clustering [Llo]. We use a variant of the clustering algorithm first proposed by [VC] which approximates a CVT in a computationallyefficient manner involving only local queries on voxels located on boundaries of pairs of clusters. The original algorithm was applied in the context of the uniform coarsening of triangle meshes, and was later extended to discrete image volumes by [DVS + ]. Results of quantitative analysis and qualitative comparison, reported in Section .., with the TurboPixels and the Graph Cut superpixels algorithms on the Berkeley segmentation dataset [MFTM] confirm that the CVT clustering algorithm is substan- tially faster than the other two algorithms while generating tessellations of comparable quality to theirs. Runtime evaluations on medical images, given in Section ., show that the total runtime of CVT clustering and Graph Cut segmentation is - times better than voxel segmentation alone with no compromises in segmentation quality for most practical purposes. The following section outlines the clustering algorithm.

 v u C i C j v u C i C j v u C i C j (a) (b) (c)

.. The clustering algorithm

Here we outline the CVT clustering algorithm of [VC] which we use to partition an image before proceeding with Graph Cut segmentation. Algorithm .. gives its listing.

The algorithm starts by randomly selecting n voxels and assigning each to one of the n projected clusters C i (lines -). It assigns all of the remaining voxels to a special null cluster representing the fact that these voxels do not yet belong to any cluster (lines -). The algorithm then proceeds to minimize the energy (.) iteratively, examining all pairs of adjacent clusters and updating voxel membership in clusters according to boundary tests for local minimum-energy configurations (lines -). If one of the clusters is the null cluster (lines -), the algorithm assigns all of the voxels in its boundary set, denoted by ∂C, which are neighbors with voxels in the nonnull cluster to the latter without performing any energy calculations, thus growing the nonnull cluster. If both clusters are nonnull, the algorithm checks whichever of the three possible voxel assignments, illustrated in Figure ., gives the largest decrease in energy and updates the current clustering configuration accordingly. Note that boundary tests rely on voxel neighborhood systems as defined in Section . for voxels in the boundary sets of clusters. The algorithm's main loop terminates either when no clusters are modified during its last iteration or when the maximum number of iterations are executed.

We note, without giving details, that [VC] formulates the energy (.) so that it depends only on the assignment of voxels to clusters, therefore cluster energy computations on line  do not involve any explicit calculation of cluster centroids. Moreover, in order to ensure implementation efficiency, cluster energy terms are stored in computer memory and updated incrementally as the configuration changes so that the energy computation for a cluster is a constant-time operation. Hence, the asymptotic time complexity of the clustering algorithm is that of its inner loop: O(∆ n e), where ∆ = max ∀i,j∈S |∂C i ∩ ∂C j | is the greatest cluster common-boundary length, and e is the



 for i ←  to n  C i ← Random (v ∈ I )  foreach v ∈ I | C i , v ∈ C i  C null ← v // Perform energy minimization  S = S ∪ {null}  iterations ←   repeat  modif ications ←   foreach {i, j} ∈ S × S | |∂C i ∩ ∂C j |   if i = null  foreach u ∈ ∂C i | ∃v ∈ ∂C j , v ∈ N u  C j ← u  else if j = null  foreach u ∈ ∂C j | ∃v ∈ ∂C i , v ∈ N u  C i ← u  else  foreach {u, v} ∈ ∂C i × ∂C j | v ∈ N u  conf init ←-{C i ← u, C j ← v}  conf i ←-{C i ← u, v}  conf j ←-{C j ← u, v}  if Energy(conf i ) < Energy(conf init ) or Energy(conf j ) < Energy(conf init )  conf new ←-MinEnergy conf i , conf j  UpdateCVT (conf new )  modif ications ← modif ications +   iterations ← iterations +   until modif ications =  or iterations = maxIterations
ratio of the number of edges to the number of vertices in the dual graph of the CVT giving a mean estimate of the number of neighbors to consider per cluster during boundary tests. The convergence of the Algorithm .. is guaranteed since the energy is positive and each modification decreases its value.



.. Evaluating the clustering algorithm

In this section, we present the results of quantitative evaluation and qualitative validation of the CVT clustering algorithm comparing its performance with those of the TurboPixels algorithm [LSK + ] and the Graph Cut supepixels algorithm [VBM].

We shall first define the two performance metrics used in quantitative analysis; undersegmentation error and boundary recall. Given a CVT-clustering of an image C = {C  , . . . , C n } indexed by S and a ground-truth partition of it T = {g  , . . . , g k }, the undersegmentation error for segment g j is defined as the total amount of surface bleeding induced by clusters which overlap the segment g j , normalized by the segment's area:

seg err j (C, T ) =          i:C i ∩g j ∅ |C i |          -|g j | |g j | (.)
It is obvious that the above fraction equals zero when all clusters overlapping the segment g j fall inside it, that is, no cluster crosses the segment's boundary. When measuring the undersegmentation error of an algorithm for a given clustering, we average the above fraction over all segments of the corresponding ground-truth partition and give the mean error. Boundary recall is defined straightforwardly as the fraction of ground-truth boundary pixels which fall within a given small distance from at least one cluster boundary. Since we seek high precision in segmentation boundary placement for medical applications, we use a distance of one pixel.

In order to stay in line with quantitative results reported in [LSK + ] and [VBM], we also use the Berkeley segmentation database [MFTM] for quantitative evaluation. This dataset comprises  generic photos of size  ×  or  ×  pixels which have been presented to human subjects both in color and grayscale modes for manual segmentation. For evaluation purposes, we select  of the  grayscale images from the test subset, and choose one ground-truth segmentation for each such image so that it does not contain small segments of few pixels. Such small ground-truth segments can create large outliers and disrupt mean undersegmentation error measurements.

Quantitative evaluations are performed with respect to increasing number of clusters in the range [, ]. In all evaluations, the temperature hyperparameter of the Graph Cut superpixels algorithm is set to , a default recommended by the author's implementation, and two iterations of the algorithm are applied. The maximum number of iterations for the CVT clustering algorithm is set to . In practice, however, the algorithm terminates much earlier. The intensity-norm weighting coefficient α in (.) is set to one, and the algorithm's performance is evaluated with respect to two definitions of the density function (.); m = , yielding a gradient and intensity adaptive clustering, and m = , yielding an intensity adaptive clustering.

Results of quantitative evaluation for TurboPixels, Graph Cut superpixels and the CVT clustering algorithms are presented in sures, averaged over all images for a given clustering resolution in Figures .d-.f.

From these graphs, we can readily see that the CVT clustering algorithm outperforms the other two in terms of execution time by a wide margin and that its boundary recall is mostly higher than those of the other algorithms. Its mean undersegmentation error is comparable to that of the TurboPixels algorithm and is better than that of the Graph Cut superpixels algorithm with lower clustering resolutions.

We can see in Figures .a-.c that both TurboPixels and Graph Cut superpixels algorithms produce a variable number of output clusters. As previously mentioned, the TurboPixels algorithm assigns remaining connected regions to new superpixels upon termination and the Graph Cut superpixels algorithm deduces the number of superpixels from the specified maximum superpixel size, although, in theory, the number of superpixels could be specified directly as the number of labels to optimize. Mean performance figures for these algorithms in Figures .d-.f are therefore reported with respect to the average number of output clusters for each input resolution. Note that the CVT clustering algorithm produces the exact number of clusters specified, which is important in applications where the required cluster budget has to be respected.

We present the results of quantitative evaluation of the CVT clustering algorithm with respect to different density function definitions separately in Figure .. Here we can see that the average performance of the gradient-based CVT clustering is slightly lower, more significantly in terms of mean boundary recall, it also takes longer to execute. These results are not surprising, as the dataset on which evaluations are performed features highly textured images with nonuniform illumination and substantial presence of shadows. The gradient-based density function, while decreasing the size of clusters, increases their density in textured areas of the image. This reduces the cluster budget for other image areas, which, in the presence of shadows or with insufficient contrast, provide poor gradient information for precise cluster placement. This point is illustrated in



. Segmentation energy redefined

In Section ., we formulated segmentation as a labeling problem, defined as the assignment of a label from a set of m labels L = {l  , . . . , l m } representing objects to be segmented to each of the variables in a set of n variables corresponding to image voxels. Henceforward, we shall assume that these variables correspond to the clusters of a CVT-clustered image C = {C  , . . . , C n } indexed by S. We shall assume furthermore that each variable i ∈ S is associated with the corresponding node in the dual graph of the CVT, G, illustrated in Now the prior probability distribution of label configurations can be expressed as an MRF of label configurations with respect to the graph G, and if image-derived likelihood densities can be defined for each label l ∈ L, then the MAP configuration yielding the best segmentation with respect to user-supplied evidence and constraints can be found by minimizing the energy function (.) with the obvious replacement of the voxel neighborhood system N i by the cluster neighborhood system Ñi . We reproduce the energy function of segmentation here for convenience.

E( ) = t i∈S D i ( i ) +   i∈S j∈ Ñi V i,j ( i , j ). (.)
Recall that the unary term D i ( i ) in (.), also called the data term, is a negative loglikelihood function derived from observed data and measures the cost of assigning a label to the variable i ∈ S, and that pairwise terms V i,j ( i , j ) encode prior information on pairwise interactions between labels assigned to pairs of neighboring sites.

In the following sections we redefine unary and pairwise energy terms of (.) to suit the Bayesian labeling of CVT-clustered images. Importantly, this redefinition maintains the clustering-resolution independence of segmentation energy, which we establish formally in Section ...

.. Likelihood-based data cost redefined

Recall that the likelihood probability for a variable i ∈ S given the label assignment i is given by the conditional probability Pr(O i = I i | X i = i ), where the observation O constitutes the spatially distinct records of image intensity levels at each variable. The conditional probability distribution of intensity Pr(I | L = l) for every label l ∈ L is estimated as Gauss-smoothed and normalized intensity histograms from user-supplied "seeds", as described in Section ... The definition of the term D i ( i ) is slightly different however. To facilitate the understanding of different cases, we define it as follows:

 D i ( i ) =                , ∃h k ∈ C i , k = i (.a) ∞, ∃h k ∈ C i , k i (.b) -ln v∈C i Pr(I v | i ), h k ∈ C i , ∀k ∈ L (.c)
where h k represents a seed of label type k ∈ L attributed to a voxel v ∈ C i and I v is the gray-level of voxel v. (.a) (respectively (.b)) corresponds to the negative loglikelihood value at maximum (zero) probability and indicates no cost (infinite cost) for assigning the label k ∈ L (any label but k) to a variable i when there is at least one voxel in the associated cluster C i marked by a seed of label type k. Recall that conditions (.a) and (.b) amount, in practice, to imposing hard constraints on configurations in the space L, and thus encourage plausible segmentations. (.c) accounts for all other variables for which there is no seed to impose a labeling decision. This is the case for the vast majority of variables, since the user is expected to introduce only a small number of seeds with respect to the number of image voxels. Note that in (.c) the likelihood probability for a variable i ∈ S, which corresponds to a cluster C i , is defined on intensity observations at image voxels v

∈ C i as Pr(O i = I i | X i = i ) = v∈C i Pr(I v | i )
under the assumption of probabilistic independence of such observations. An alternative specification of likelihood probability Pr(O i = I i | X i = i ) = Pr(I i | i ) defined on intensity observations at clusters, the intensities of which could be calculated as the mean intensity of its voxels, would not take intensity variance within clusters into account and furthermore the resulting energy would not maintain resolution independence.

.. Prior redefined

We weight the pairwise prior term V i,j ( i , j ) in (.) by the length (area, in D) of the common boundary (surface) of adjacent clusters |∂C i ∩ ∂C j | so that the sum of pairwise terms in the energy of some configuration E( ) for any pair of labels would be equal to the length (area) of the common boundary (surface) between the corresponding pair of objects multiplied by a constant, equal to the length of the shortest path between corresponding nodes in the prior graph of adjacency relations as defined in Section .. The definition of the pairwise prior thus becomes:

V i,j i , j = |∂C i ∩ ∂C j | (a, b) , i = a, j = b. (.)
Note that clusters sharing longer boundaries (larger surfaces, in D) would prefer similar labels, or label assignments incurring smaller penalties.

In addition to its theoretical soundness, as shown in the following section, results of quantitative evaluation presented in Section . establish the practical advantage of this weighting scheme compared to the empirical inverse-distance prior weighting method frequently encountered in MRF-based segmentation literature, for example in [BF] and [CCB + ].



.. Energy invariance

We have mentioned previously that the definitions of data and prior terms of (.) given in Sections .. and .. respectively maintain the clustering-resolution independence of segmentation energy. This means that the energy of identical label configurations on CVT clusterings of different resolution is invariant. In particular, given CVT clusterings of different resolution, the energies of segmentations corresponding to identical MAP solutions * are identical.

The following theorem establishes the energy invariance property with respect to data energy formulations disregarding hard constraints as defined by (.a) and (.b). In the presence of hard constraints, energy invariance no longer holds. The corollary of the energy invariance theorem shows that when hard constraints are considered, the energy of a MAP solution * on a higher resolution CVT clustering is greater than the energy of an identical MAP solution obtained on a clustering of lower resolution. Let us write the sums of unary energies of * for T and T as follows:

E  ( * ) = t i∈S D i ( i ) = t i∈S M D i ( i ) + t i∈S M D i ( i ) (.a) E  ( * ) = t i∈S D i ( i ) = t i∈S M D i ( i ) + t i∈S M D i ( i ). (.b)
We shall show that the respective sums on the right hand side of (.a) and (.b) are equivalent.



If T = T , then M = M therefore |M| = |M |.
Expanding both sides on voxels

i∈S M v∈C i  = i∈S M v∈C i 
and summing over the negative log-likelihood for every v ∈ C i , we get

i∈S M v∈C i -ln Pr(I v | i ) = i∈S M v∈C i -ln Pr(I v | i ).
Rewriting the sum of logarithms on both sides as the logarithm of the product we get

i∈S M -ln v∈C i Pr(I v | i ) = i∈S M -ln v∈C i Pr(I v | i )
which, by definition (.c), gives

i∈S M D i ( i ) = i∈S M D i ( i ).
Similarly, we can show that

i∈S M D i ( i ) = i∈S M D i ( i ). Hence E  ( * ) = E  ( * ). (.)
Now we show that the sums of pairwise energies of * for T and T are equivalent. As illustrated in Figure ., denote the boundary between the pair of segments in T and T by F and F respectively, and let us define the corresponding sets of cluster pairs falling on either side of each boundary as

B = {(i, j) | C i ∈ M, C j ∈ M} and B = {(i, j) | C i ∈ M , C j ∈ M }. If T = T , then F = F therefore |F| = |F |
which can also be expressed as the length of the common boundary between the pair of segments in either segmentation

|∂M ∩ ∂ M| = |∂M ∩ ∂ M |.



Expanding both sides on common boundaries between pairs of clusters on either side of each segmentation boundary, we get i,j∈B

|∂C i ∩ ∂C j | = i,j∈B |∂C i ∩ ∂C j |.
Let l ∈ L denote the label of segments M and M and l ∈ L that of segments M and M .

Multiplying both sides by a constant k = V i,j

( i = l, j = l)/   i,j∈B |∂C i ∩ ∂C j |V i,j ( i , j ) =   i,j∈B |∂C i ∩ ∂C j |V i,j ( i , j )
and expanding both sums on cluster neighborhood, we get

  i∈S M j∈ Ñi j∈S M |∂C i ∩ ∂C j |V i,j ( i , j ) =   i∈S M j∈ Ñi j∈S M |∂C i ∩ ∂C j |V i,j ( i , j ). Hence E  ( * ) = E  ( * ). (.)
From (.) and (.) we get E( * ) = E ( * ), thus completing the proof.

Corollary ... Given two CVT clusterings C and C of the image I indexed by S and S , and two identical segmentations on them T and T resulting from minimizing (.) in some labeling space with respect to hard constraints and label likelihood-based data expressed by (.) and pairwise prior defined by (.).

If |C| < |C |, then E( * ) < E ( * ).
Proof. Denote the energy of the optimal labeling * with respect to label likelihoodbased data (.c) and pairwise prior(.) by E( * ). Let S H index the set of variables subject to hard constraints

S H = {i ∈ S | ∃h k ∈ C i , k ∈ L}.
The energy of * with respect to label likelihood and hard constraints-based data (.) and pairwise prior(.) is

E c ( * ) = E( * ) - i∈S H -t ln v∈C i Pr(I v | i ) (.) because ∀i ∈ S H , D i ( i ) = . Obviously, E c ( * ) < E( * ). Let |C| = m and |C | = n. If m < n, then the average cluster size of C is greater than that of C  m i∈S |C i | >  n i∈S |C i |.
Therefore, for S H and S H i∈S H

|C i | > i∈S H |C i |.



Expanding both sides on voxels

i∈S H v∈C i  > i∈S H v∈C i 
and summing over the negative log-likelihood for every v ∈ C i , we get

i∈S H v∈C i -ln Pr(I v | i ) > i∈S H v∈C i -ln Pr(I v | i ).
Rewriting the sum of logarithms on both sides as the logarithm of the product we get

i∈S H -ln v∈C i Pr(I v | i ) > i∈S H -ln v∈C i Pr(I v | i ).
From the Energy Invariance theorem we have E( * ) = E ( * ). Multiplying both sides of the above inequality by -t and adding E( * ) to both sides we get

E( * ) - i∈S H -t ln v∈C i Pr(I v | i ) < E( * ) - i∈S H -t ln v∈C i Pr(I v | i ) hence E c ( * ) < E c ( * )
, thus completing the proof.

. Conclusions

In this chapter, we have addressed the complexity of Graph Cut image segmentation on high-resolution input, in particular on D medical images often comprising several millions of voxels. We presented an image simplification approach to partition such an image prior to segmentation by an efficient intensity and gradient-adaptive centroidal Voronoi tessellation (CVT) method. It allows to control the tradeoff between the conflicting criteria of cluster compactness and structure boundary adherence, and furthermore produces dense partitions of relatively small clusters near intensity edges helping towards the fine-grained placement of subsequent segmentation boundaries.

We have redefined the Bayesian labeling problem and the Graph Cut minimization of its energy according to the graph of CVT clusters. As we have established formally, our redefinition maintains the clustering-resolution independence of the segmentation energy, which means that identical image segmentations on CVT partitions of different resolution have identical energies.

Comparative evaluation results of our CVT image clustering method with closelyrelated TurboPixels and Graph Cut superpixels algorithms on the Berkeley segmentation dataset confirmed the advantages of our method in terms of both runtime and quality of produced partitions. In the next chapter, quantitative and qualitative evaluation of our CVT-clustered segmentation approach on synthetic, simulated and real medical images confirm that, compared to a direct voxel segmentation approach, the clustering step improves the total runtime and memory footprint of the segmentation process up to an order of magnitude without compromising the segmentation quality for most practical purposes.



Chapter 

Evaluations and applications

. Introduction T his chapter summarizes the results of quantitative analysis and qualitative validation of the vicinity prior-based Graph Cut segmentation approach proposed in Section . on synthetic, simulated and real medical CT and MR images. In all evaluations, comparisons are drawn with the standard Potts prior-based segmentation and in some cases also with the adjacency prior-based approach introduced in Section .. Let us first outline the conditions under which evaluations presented in this chapter are performed.

For comparison purposes with clustered image segmentation as formulated in Section ., we perform some of the evaluations on image pixels (voxels). We apply the basic framework of Bayesian labeling outlined earlier in Section . by letting variables i ∈ S correspond to image pixels (voxels) v ∈ I and defining an -connected (connected) neighbourhood N i for each variable (respectively in D).

In performance measurements with respect to varying number of clusters, we use  clustering resolutions, reported in different figures and tables in this chapter as percentages of the number of image voxels. In all evaluations, we set the coefficients m and b of the clustering density function (.) to m = b = , although we have seen in Section .. that in some cases, especially with images having high noise levels, careful setting of these parameters can improve the clustering aspect and yield higher boundary recall and lower undersegmentation error.

We minimize the energy (.) using the aforementioned three prior models with the Expansion Moves algorithm discussed in Section ... In all evaluations, we set the algorithm's maximum number of cycles to . In practice, however, the algorithm converges to a local minimum rapidly in - cycles. We note that all our software implementations use the multilabel energy optimization library developed by the Computer Vision Research Group at the University of Western Ontario [BVZ, KZ, BK].

Before presenting the results in detail, we shall define performance metrics used in quantitative analysis in the following section.



. Performance measures

We measure the quality of segmentation with respect to ground-truth via the widelyused Dice similarity metric (DSM) [Dic]. Denote the set of voxels labeled with l ∈ L in the segmented image I and the ground-truth T by S l and T l respectively, the DSM for label class l ∈ L is defined as:

dsm l (I , T ) =  |S l ∩ T l | |S l | + |T l | . (.)
It is obvious that the DSM measures the overlap between the segmentation and the ground-truth. It takes its values in [, ] where  indicates a perfect match and  no match at all. In addition to the DSM metric for individual label classes, we employ an "overall" segmentation quality metric calculated from mean weighted DSM measures for all labels in L:

dsm L (I , T ) = l∈L dsm l (I , T ) |T l | |T | . (.)
The two other classification performance metrics we use in quantitative analysis are sensitivity and specificity, defined respectively as the true-positive and the truenegative rates for individual label classes [Faw]:

sens l (I , T ) = |S l ∩ T l | |T l | , (.) spec l (I , T ) = k∈L |S k ∩ T k | -|S l ∩ T l | |I | -|T l | . (.)
Similarly to (.), overall sensitivity and specificity metrics are defined from mean weighted sensitivity and specificity measures for all labels in L.

In order to justify the significance of observed improvements brought by the vicinity prior on the aforementioned performance measures, we complement comparison results on medical case studies with statistical testing according to Welch's t test [Wel], which is an adaptation of Student's t-test intended for use with samples having possibly unequal variances. For two such samples, the test statistic

t = µ  -µ  σ   n  + σ   n  (.)
approximately follows Student's t-distribution with degrees of freedom

df = σ   n  + σ   n   σ   n   (n  -) + σ   n   (n  -) (.)
where µ i , σ  i and n i are respectively the i th sample mean, variance and size. Given t and df , we calculate a p-value using the t-distribution in a standard fashion and perform  We evaluate the precision of segmentation boundary placement with an error metric measuring the Euclidean distance of each point on the surface of a segmented structure from the surface of the corresponding ground-truth structure. Let M l S and M l T be the triangle meshes representing isosurfaces of structure volumes labeled by l ∈ L in the segmented image I and the ground-truth T . We obtain such meshes by applying the Marching Cubes algorithm to labeled volumes [LC]. The error at x ∈ M l S is given by:

err(x) = inf {d(x, y) | ∀y ∈ M l T }. (.)
Later in this chapter, we give the cumulative histogram of errors calculated for all points on the surface of a given structure. The histogram representation gives a clear indication of the percentage of segmented surface lying within a given distance from the true surface, and consequently allows us to understand how precise the segmentation is in retrieving true structure boundaries.



. Synthetic phantom and dataset

In this subsection, we present the results of quantitative and qualitative evaluation of segmentation on a dataset generated from a synthetic phantom we have created for the very purpose inspired by the Shepp-Logan phantom [SL] which has been used in tomographic reconstruction. We have already encountered this phantom in Chapter . For convenience, its structure is reproduced in Figure .a. Due to the small number of structures, their simplistic spatial arrangement and intersection patterns, the Shepp-Logan phantom, illustrated in Figure .c, is not well-suited for evaluating a segmentation approach using prior information on spatial relationships of structures.

The layout of the  structures our phantom features mimics the spatial configuration of some thoracic and abdominal organs in a coronal cross-section of a CT scan. This can be seen in Figure .b. Roughly speaking, identically-colored phantom structures  and  correspond to lungs, right and left,  and  to pulmonary arteries in each lung,  and  to the liver and the spleen and  to the cavity of the digestive tract, which has zero intensity in areas filled by air.  Visceral fat, , fills the space between these structures and the whole is covered by a thin envelope of skin, , and is placed against a background of air, , which is distinct from the digestive tract cavity, although both have identical intensities. The intensity levels of the  structures are set from mean intensity values of corresponding structures in a CT image for which ground-truth estimates for these structures are available.

Let us outline the conditions under which evaluations discussed in this section are performed. Graph Cut segmentation with Potts, adjacency and vicinity priors is applied to the pixels as well as to clusterings of  ×  images generated from the synthetic phantom by the introduction of additive zero-mean Gaussian noise of variable standard deviation σ sweeping the range [, .] with a step of .. For each noise level, we generate  images clamping intensity values to [, ]. Of ranges swept by the temperature parameter t of the energy (.), we report results in the representative window [, ] sampled at .. The reason behind the choice is that in preliminary experiments, where we probed .-sampled temperature ranges with an upper bound much higher than , we did not find any subranges where Potts prior segmentation, as measured by mean overall DSM, outperformed the best result of vicinity prior segmentation over all subranges. Moreover, segmentation quality with both priors usually deteriorated with temperature settings t > . The set of "seeds" used to estimate likelihood densities as well as to impose hard constraints on segmentation in all evaluations on the synthetic dataset is given in Figure .. The number of seed-marked pixels for each label is in the range [, ] and roughly proportional to the size of the corresponding structure, otherwise hard constraints of equal image support for all labels would create uneven bias of segmentation ease for smaller structures, like  and .

In the remaining of this section as well as in Figures ., ., ., ., . and . any references to DSM are actually for overall DSM measures as defined by (.).

 The digestive tract cavity is not visible in the same plane as the other structures in the coronal cross- Comparisons are drawn between best segmentations obtained with either prior on all images generated at a given noise level for a fixed temperature t = . and over all temperature settings. The vicinity prior was defined according to the graph model of structure adjacency given in Figure .b on page . Improvements and consistency in segmentation results brought by the vicinity prior, which are furthermore confirmed by DSM measures, are clear on both pixel and cluster segmentations. We can see that the pair of structures  and  having identical intensity levels and the pair  and  having close intensity levels, with the latter structure having a rather close level to structure , are always segmented correctly as distinct structures when the vicinity prior approach is able to retrieve them. The approach also shows better ability to retrieve thinner structures ,  and the envelope  in particular.

                                                         (d) (e) (f) (g) (h) (i) (j) (k) (l) (a) (b) (c)
In Figures .a-.c and Figures .a-.c, we give mean DSM maps for segmentation with Potts, adjacency and vicinity priors applied to clusters and pixels respectively. In the latter group, we calculate mean DSM values for segmentations corresponding to a particular temperature setting on all images generated for a given noise level. In the former group, such mean DSM values corresponding to a particular clustering resolution are calculated at temperature settings producing the best segmentation on all images generated for a given noise level. provements brought on mean DSM by the vicinity prior are remarkable, particularly for noise levels up σ = .. Note that the adjacency prior segmentation displays overall performance patterns similar to its vicinity prior counterpart albeit being less robust to temperature settings and more sensitive to noise levels higher than σ = . as it can be observed in Figures .c and .c respectively. Note that the low mean DSM values corresponding to low temperature settings in Figures .a, .b and .c are due to high regularization by the prior since data terms do not contribute as much at low temperatures. In particular, the minimum energy configuration corresponding to the setting t =  is a constant labeling, since V (a, a) =  for any a. Recall that the vicinity prior incurs higher regularization than the Potts prior since its highest penalty in this case study is four times greater than its Potts counterpart. To observe the difference between the two priors on a finer scale, we give a magnified view of DSM maps on the low temperature range [., .].

In order to understand the effect of additive Gaussian noise parameter σ alone, we give mean DSM comparisons of Potts, adjacency and vicinity prior segmentations on pixels and  clustering resolutions with respect to increasing noise levels at a fixed  temperature t = . (Figures .d-.f, .d) and at the temperature setting yielding the best segmentation for all images generated at a given noise level (Figures .g-.i, .e). As for reliability curves, they are given in Figure .f for pixel segmentation, and Figures .j-.l for segmentation corresponding to the  clustering resolutions. Quantitative measures in both groups of evaluations confirm that the vicinity prior segmentation approach outperforms the Potts prior segmentation, and the adjacency prior approach to a lesser extent.

We also evaluate the robustness of Potts, adjacency and vicinity prior segmentations to clustering resolution. Results of this evaluation for  noise levels are presented in Figures .a-.c. We can readily see that with vicinity and adjacency priors, segmentation is quite robust to clustering resolution, even at a noise level as high as σ = ., although it attains lower mean DSM values with the latter prior. This means that similar segmentation results can be obtained using lower clustering resolutions thus improving both runtime and memory footprint. For this evaluation, we calculate mean DSM measures for segmentations corresponding to the best temperature setting on all images generated at given noise level.

Lastly, in Figure . we present the results of evaluating different weighting schemes of pairwise prior terms in the segmentation energy (.) for  clustering resolutions. Let C i and C j be a pair of adjacent clusters sharing a common boundary (a surface in D), and denote their centroids by c i and c j respectively. We compare the following four weighting schemes. ) Cluster common-boundary length weighting (.) where the pairwise prior is weighted by |∂C i ∩ ∂C j |. As shown by the Energy Invariance theorem .., this weighting scheme along with the likelihood definition (.) maintain the clustering-resolution independence of segmentation energy. ) Inverse centroiddistance weighting frequently encountered in MRF-based segmentation literature, for example in [BF] and [CCB + ], where the pairwise term is weighted by /d(c i , c j ). ) Both weighting schemes applied together. ) No weighting. As it can be seen from graphs in Figure ., the common-boundary length weighting scheme is a clear winner independently of the clustering resolution thus confirming the theoretical result.

. Simulated MRI: BrainWeb

In this section we summarize results of quantitative and qualitative validation of Graph Cut segmentation with Potts and vicinity priors on simulated MRI images. We first give quantitative evaluation results on the original BrainWeb simulated MRI dataset [CZK + ] and compare them to  tissue-class voxel segmentation results reported by [GTC] who use an iterative variant of the Expansion Moves algorithm with the Potts prior model. Then we present a segmentation case study for a head image from the second generation BrainWeb simulated MRI dataset [AEC].



.. Original BrainWeb simulated MRI dataset

As mentioned, we perform a comprehensive quantitative evaluation on the original BrainWeb simulated MRI dataset [CZK + ]. In the creation of this dataset, a semirealistic anatomical head and brain phantom was used by an MRI simulator of T-, Tand PD-weighted modalities to generate images according to MR acquisition physics with  noise levels (%, %, %, %, % and %) and  intensity inhomogeneity levels (%, % and %). For evaluation purposes, we use  images from the T-weighted modality covering the entire noise and inhomogeneity ranges. Image volumes are of size  ×  ×  voxels with isotropic  mm voxel size. The phantom itself used to generate images comprises  structure and tissue classes.

It should be noted that the phantom exhibits certain anatomical inconsistencies and structural artifacts. As shown in Figure ., some structures corresponding to tissue classes "muscle" and "skin" are located inside the skull within the volume of tissue class "cerebrospinal fluid" and next to "gray matter" tissue. We therefore define our graph-based anatomical model of structure adjacency, given in Figure .a, as a compromise between phantom imprecisions and true anatomical properties. We note that we use  out of the  tissue classes the phantom features, as we merge the tissue class "glial matter" into the class "gray matter" without affecting anatomical consistency.

In Tables .a, .b and .c, we give means and standard deviations for overall and individual-class DSM, sensitivity and specificity measures calculated on voxel and cluster segmentations corresponding to best temperature settings for all images in the dataset. We define the vicinity prior according to the graph model presented in We compare our results on the original BrainWeb simulated MRI dataset with those reported by Gorthi et al. in [GTC] which is the only recent work we are aware of that evaluates the Expansion Moves algorithm for brain tissue classification using the Potts prior model. There are, however, some differences between their approach and ours. Firstly, they use Gaussian likelihood distributions the parameters of which are learned during an iterative application of the algorithm. Secondly, and more importantly, using the phantom as a volume mask they retain only brain structures i.e. cerebrospinal fluid, gray and white matters and discard other head tissues [Gor]. Thus they formulate segmentation as a multilabeling problem with  labels plus  for the background. This not only reduces computational complexity, but also makes the segmentation task easier since other head structures, some of which have intensity levels similar to those of brain structures such as muscle and gray matter tissues, are eliminated beforehand. The aforementioned  brain structures, on the other hand, exhibit good contrast in images in general, and have reasonably distinct mean intensity values which further facilitate the segmentation task. In our evaluation, we use the same subset of the original BrainWeb dataset as [GTC] corresponding to the  T-weighted images generated with inhomogeneity levels % and % and all noise levels. In Figure ., we compare means and standard deviations for DSM, sensitivity and specificity measures for the  brain structures from full-image vicinity prior segmentations with corresponding measures in [GTC]. Unsurprisingly, results we report are not superior. Specialized approaches are expected to, and indeed should, produce better results. Therefore, the present comparison should be viewed as an indication of how close a generic approach can come to a specialized one on this problem instance. Nevertheless, a practical advantage of our approach is that no manual tissue removal is required to segment the brain volume.



.. Second-generation BrainWeb simulated MR image

In this section, we present a segmentation case study of a single image from the second generation BrainWeb simulated MRI dataset [AEC]. This dataset features  simulated T-weighted head MR images generated by an MRI simulator at a fixed noise level of % according to subject-specific phantoms. These phantoms are superior to that of the original BrainWeb dataset in terms of anatomical precision and consistency. They comprise  tissue classes of which we use  merging tissue classes "dura mater" and "around fat" respectively into classes "cerebrospinal fluid" and "muscle" without affecting anatomical consistency. The image we choose for evaluation purposes is that of subject no. . Its volume is of size ×× voxels with isotropic  mm voxel size. We apply Potts and vicinity prior-based Graph Cut segmentation to a % clustering and to the voxels of the image at a temperature setting of t = .. Tables .a and .b give overall as well as individual-class DSM, sensitivity and specificity measures for cluster and voxel segmentations with either prior. In these tables, vicinity prior mean performance measures indicating statistically significant improvements at the level α = . are highlighted. Qualitative validation and comparison of cluster and voxel segmentations is given in Figures . These evaluations confirm that the vicinity prior outperforms the Potts prior in Graph Cut segmentation on the simulated head MR image, and that with -fold reduction in image size, we are capable of producing a high-quality result which, for most practical purposes, is virtually identical to that obtained on the full-resolution image with two-fold improvements on runtime and more than three-fold savings of memory. In this subsection, we present a segmentation case study of in vivo acquired thoracicabdominal D CT scan realized during the arterial phase in inhaled position. The image comes from the D-IRCADb dataset [IRC] and is accompanied by ground-truth estimates of  anatomical structures, which are basically manual segmentations performed by clinical experts. Of these  estimates, we use  in our evaluations for the following two reasons. ) Some of these estimates do not correspond to any particular structure but to many, for example the "skin" tissue class (refer to Figures .d-.f) which comprises skin, fat, muscle and visceral fat among other tissues. ) Some estimates are partial, for example those of intestines for which ground-truth is available only for the colon part. We think that these  structures of different size, shape and topology are sufficient for evaluation purposes.



(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
The image volume we use is of size ×× voxels with .×.×. mm voxel size. We segment all visible structures in this image, including those for which no ground-truth estimate is available, applying Graph Cut segmentation with Potts and vicinity priors to % and % clusterings of the image at a temperature setting of t = .. Seeds, illustrated in Tables .a and .b give overall as well as individual-class DSM, sensitivity and specificity measures for both clustered segmentations with either prior. In these tables, vicinity prior mean performance measures indicating statistically significant improvements at levels α = . and . are highlighted with blue and red respectively. Runtime and memory footprint figures for segmentations on both clustering resolutions are reported in Table .. Qualitative validation and comparison is given in Figure . on the  principal cross-sections of ground-truth, vicinity prior and Potts prior-segmented images for the clustering resolution of %. Disagreements between ground-truth and Potts prior segmentations are encircled. Notice in particular how right and left lungs are correctly segmented as distinct organs with vicinity prior. Note that, for the reasons mentioned above, ground-truth estimates can be relied on for qualitative inspection only for structures listed in Table .. Figure . gives D views of simplified surface meshes generated via the mesh simplification method of [VC] from different structure volumes of the vicinity prior-segmented %-clustered image.

Lastly, in Figure . we give the surface error measurements of the aforementioned Potts and vicinity prior-based % and % clustered image segmentations for  structures, excluding adrenal glands and the hyperplastic lesion since they are too small for their surface error measurements to be interesting. Cumulative error histograms of structures for each segmentation are broken into  graphs to improve readability. Presented graphs indicate the fraction of segmented surface lying within a given distance from the true surface. We can see that for reasonably low error levels of - cm,  vicinity prior segmentation retrieves a considerably larger fraction of structure surfaces than its Potts prior counterpart suffering much lower maximum error. Also take note of the similarity between error measurements of vicinity prior segmentations applied to the %-clustered and the %-clustered image. Nevertheless, vicinity prior segmentation produces some noticeable outliers, namely for the volumes of left lung, heart and spleen. In fact, the CT scanner table in the image is erroneously segmented as part of the body. Such outliers constitute a small percentage of the labeled volume, and can be corrected via "corrective" seeds, as mentioned in Section .., or by a separate modeling of the scanner table with its own set of labels and pairwise relations with other structures in the image. Note that surface error measures alone do not constitute an indicator of segmentation quality. Consider, for example, the case of erroneous segmentation of lungs with the Potts prior in .a indicates little agreement with the corresponding ground-truth estimate. The higher surface error for the left lung in Figure .a is due to most of the right lung volume being segmented as left lung, and consequently its surface having larger per-point minimum distance (.) with respect to the left lung ground-truth surface.

Like many Graph Cut approaches, our generic method is prone to shrinking bias affecting thin elongated structures. Take note of fragmentary abdominal and hepatic arteries and portal veins in Figure .b. This is due to the sum of pairwise costs along the boundary of a vessel being higher than the cost of a shortcut across it when likelihood terms do not incur a strong preference for the vessel's class. A number of solutions have been proposed for this problem in Graph Cut segmentation. [KB] segment vessels using surface functionals combining image-gradient flux and image-based Riemannian area, [VKR] improve a result with user-supplied connectivity constraints indicating object regions which need to be connected, [PWL] propose a Hessianbased "vesselness" likelihood for segmenting hepatic arteries and portal veins in CT.

We conclude this section recalling that quantitative and qualitative evaluation confirm the significant improvements brought by the vicinity prior to this instance of multiobject segmentation even when it is applied to a coarsely clustered image, and that doubling the clustering resolution does not improve the quality but comes at the cost of doubling the runtime and increasing memory footprint by one and a half times.



(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

.. Visible Human female CT

In this section, we present a case study of complete skeleton segmentation in a highresolution CT image reconstructed from the Visible Human female dataset. Due to good contrast, CT is frequently the imaging modality of choice for bone structures.

The Visible Human dataset is the result of an initiative by the U.S. National Library of Medicine (NLM) to create a detailed digital image library of the normal adult male and female anatomy [SASW]. The dataset comprises radiographs, axial and coronal MR and axial CT images and high-resolution cryosectional photographs all acquired on donated bodies of two deceased male and female individuals, considered representative of an evaluated population. All images in both male and female datasets are provided in the convenient PNG format as -bit RGB, for cryosections, and as -bit grayscale, for radiological images. A representative view is given in Figure ..

Since its creation, the Visible Human dataset has generated significant interest in medical, image processing and computer vision research communities, in particular for the creation of anthropomorphic phantoms for applications in health physics and nuclear medicine, such as whole-body and individual-organ radiation dosimetry for radiation protection, radiotherapy and medical imaging. Computational procedures of such studies usually involve Monte Carlo simulations to track radiation interactions and energy deposition in the body. The accuracy of such simulations largely depend on the precision and the realism of the anatomical representation, especially for internal organs. Among digital phantoms created from the Visible Human male dataset, we mention the NURBS-based cardiac-torso (NCAT) phantom [Seg] and voxel-based VIP-Man [XCB] and Voxel-Man [PHP + ] phantoms. Me mention that all of these phantoms were created from cryosectional photographs mostly by manual contouring and volume editing complemented by color-space classification and thresholding. We choose to segment only the complete skeleton from Visible Human female CT dataset scanned from the fresh cadaver for two reasons. ) We have not attempted to segment any other structures because of the lack of ground-truth estimates to validate any results on soft tissue and internal organs, as the aforementioned phantoms were created from the Visible Human male dataset. The skeleton, on the other hand, especially a D surface representation of it, is relatively easy to validate by qualitative inspection. ) Slices from the Visible Human male CT dataset have variable thickness and are missing  slices from just superior to the knees to superior to the ankles as well as an one from the abdominal area. The Visible Human female CT dataset, which comes in  slices, is complete and has uniform  mm thickness throughout. Pixel size is variable though, ranging from . mm in the head and neck area to . mm in the torso. We reconstruct a  ×  ×  D volume from D slices by rescaling and padding images and setting voxel size to . × . ×  mm. Given computer memory resources at our disposal, segmentation applied to the voxels of this image is an impossible task, therefore, we would also like to demonstrate the possibility of accomplishing it by clustering the image first.



Naturally, we formulate skeleton segmentation as a binary Bayesian labeling problem and obtain the best segmentation with respect to user-supplied evidence and constraints by minimizing the energy (.). The set of labels is L = {, }, where "" represents the sought object, that is, all bone structures, and "" the background, which corresponds to everything else in the image. For pairwise energy terms, the Potts prior suffices to ensure the spatial consistency of this instance of binary labeling. We weight pairwise terms, in addition to the cluster common boundary length mentioned in Section .., by a smooth contrast function the value of which is inversely proportional to the mean intensity difference between neighbouring clusters. It is defined as:

f (I i -I j ) = e -  I i -I j σ  (.)
where I i and I j are the intensities of clusters i, j ∈ S calculated from mean intensities of their voxels respectively and σ can be estimated as the sample standard deviation from the image so that the higher it is, the more the segmentation result is spatially  regularized and its boundaries are smooth. Its effect is demonstrated in Figure ..

The data terms of the segmentation energy (.) follow the definition in (.) on the set of labels L = {, }. Thus, the segmentation is cast as a binary optimization problem which can be solved exactly, as shown in Section .., by a single application of a maximum flow algorithm to find the minimum cut of a directed weighted graph representing the energy (.). In this case, the graph corresponds to the CVT dual graph on image clusters G = S ∪ {s, t}, Ẽ with the addition of two terminal nodes. We follow the graph construction method presented in Section ... For completeness, the edge weights of this construction are summarized in Table .. To obtain the minimum cut of the graph thus defined, we use the maximum flow algorithm proposed by [BK], which we discussed in Section ... Even though its authors do not give an asymptotic bound on the running time, the algorithm has been empirically shown to outperform standard Push-relabel algorithms, generally regarded as the fastest by the combinatorial optimization community, on vision tasks formulated as labeling problems.

|∂C i ∩ ∂C j | e -  I i -I j σ  {i, j} ∈ Ẽ ∨ {j, i} ∈ Ẽ (s, i)  ∃h k ∈ C i , k = "" ∞ ∃h k ∈ C i , k = "" -t ln v∈C i Pr(I v | i = "") h k ∈ C i , ∀k ∈ L (i, t)  ∃h k ∈ C i , k = "" ∞ ∃h k ∈ C i , k = "" -t ln v∈C i Pr(I v | i = "") h k ∈ C i , ∀k ∈ L



. Runtime evaluations

We conclude this chapter with evaluations of image clustering, graph construction and Expansion Moves numerical optimization runtimes giving also total runtimes and global program memory footprint figures in terms of the maximum resident size of the process. The clustering step would obviously be less interesting should the average time to cluster an image and segment it subsequently be longer than applying segmentation to the image's voxels directly, regardless of the higher memory footprint. In Figures ., . and . respectively, we give evaluation results for three images of increasing size; a  ×  pixel D image from the synthetic phantom dataset (Section .), the  ×  ×  voxel simulated MR image from the second generation BrainWeb dataset (Section ..) and the  ×  ×  voxel thoracic-abdominal CT image from the D-IRCADb dataset (Section ..). For Expansion Moves optimization step, the temperature hyperparameter was set to t = . for all executions. We can see from the bar graphs in Figures ., . and . that the clustering runtime is inversely proportional to the clustering resolution, as previously encountered in Section ... To see how this follows from the asymptotic bound on the clustering algorithm runtime we gave in Section .., recall that the connectivity of the dual graph of the CVT is inversely proportional to the clustering resolution, as shown in Figure ., and that higher cluster densities yield shorter mean cluster commonboundary length. Graph Cut data structure construction and optimization runtime and program memory footprint figures are, unsurprisingly, directly proportional to the clustering resolution. We note that graph construction on relatively high clustering resolutions (%-) takes as long as or longer than graph construction for voxel segmentation because of the longer calculation time involved in evaluating segmentation energy data terms for image clusters from its constituent voxels (.).

Total runtime bar graphs in Figures ., . and . confirm that the clustering step does indeed accelerate the entire segmentation process compared to an application on the image's voxels directly. Thus, for practical clustering resolutions in the range [%-%], two to ten-fold improvements in total runtime can be observed, with better speedups on higher-resolution images naturally. On such images, up to ten-fold memory savings can be achieved.



. Discussion and conclusions

Qualitative and quantitative evaluation on synthetic, simulated and real medical images in this chapter confirmed the advantages of our vicinity prior model over the standard Potts model for Graph Cut multiobject segmentation counting up to a few dozen of objects. Improvements brought by our model were remarkable particularly in the correct segmentation of distinct objects having identical intensities, the accurate placement of structure boundaries and the robustness of segmentation with respect to clustering resolution. These results are also important from a theoretical perspective.

Recall that the Expansion Moves algorithm's bound on local minima is equal to twice the ratio of the largest nonzero pairwise potential to the smallest (.). Even though the Potts model yields the best bound, equal to , it was outperformed by the vicinity prior in all evaluations, which has a worse bound in all cases. For example, the vicinity prior defined after the graph model in Figure . gives a bound of .

All graph models of adjacency that were used to derive vicinity priors were defined manually on per-dataset basis according to visible anatomy in the image and the chosen level of detail. As suggested in Section ., such graph models could also be learned from training data or several model variants. Increasing the detail of a model simply amounts to representing new structures, in particular those on finer scales, by introducing new labels and associated graph vertices, and properly defining their neighbourhood with respect to other vertices. In this way, it would have been possible, for example, to represent right and left ventricles in brain MRI segmentation case studies.

Evaluations of runtime and memory footprint confirmed that the clustering step is a worthwhile undertaking, particularly for higher-resolution D images. Compared to the voxel-based approach, it accelerates the total segmentation process and reduces memory overhead remarkably. The robustness of vicinity prior segmentation with respect to clustering resolution means that in general there are no advantages of applying segmentation on image voxels or fine-grained partitions, especially that in this case the labeling space has a higher dimension and potentially many more local minima that are hard to escape. In practice, the clustering resolution can be set according to a predefined budget based on the storage and processing capacities of computing platforms.



Chapter 

Perspectives of future work

. Mission restatement I n this work, we have developed a general Graph Cut-based semiautomatic multiobject image segmentation method principally for use in routine medical applications ranging from moderately difficult tasks involving few objects in D images to fairly complex near whole-body D image segmentation of few dozens of objects for the creation of subject-specific models. We have shown that the flexible formulation of the method allows its straightforward adaption to a given application context, which does not necessarily have to be medical. In particular, spatial relationships between objects in a given problem instance can be accounted for by a simple redefinition of the general graph-based vicinity prior model we propose, and the segmentation algorithm can be tailored to the runtime requirements of the application and the online storage capacities of the computing platform by an efficient and a controllable clustering of the input image, which can be carried out in offline mode prior to presenting the image to the user for initializing the segmentation algorithm.

In order to save the reader the tedium of reiterating contributions and closely related conclusions already addressed in respective sections of chapters where they were presented and discussed, the current chapter shall point at venues of future research discussing ideas which have not been touched upon. Given the method described in this volume and our experimental setup, future research directions encompass near and long-term objectives.

. Near-term objectives

Such objectives include alternative energy definitions and algorithm implementations which do not require addressing the theoretical properties of the segmentation method.

Combination of energy models. Due to the flexibility of the energy-based theoretical framework of Graph Cut segmentation, different likelihood and prior energy models can be combined. In particular, the vicinity prior can replace the Potts model in some segmentation methods we discussed while presenting related work, such as [VM, GTC] which, respectively, use a shape similarity measure for prior and Gaussian intensity distributions for likelihood-based data energies. In our work, we have used intensity histograms in such data energies since our focus was on the prior formulation. Nevertheless, it is interesting to investigate the influence of more expressive intensity distribution models, such as Gaussian mixtures, even though they would require an appropriate learning support.



Algorithm optimization. The clustered-image Graph Cut segmentation process can be further accelerated by using a parallel implementation of a maximum flow algorithm, typically of the Push-relabel variant [BS], which would take advantage of symmetric multiprocessor architectures of modern computers. The parallelization of the clustering algorithm is more involved however. We have experimented with a rudimentary parallelization of the algorithm on image subvolumes, but the lack of congruence between clusters at the boundaries of adjacent subvolumes created blocky partition artifacts which were too degrading of the subsequent segmentation quality for the resulting speedup to be worthwhile.

Interactive segmentation. All results reported in this work relied on an implementation that executes in batch mode; user-provided seeds are collected, passed to the segmentation program which performs clustering and Graph Cut optimization writing the segmentation result onto a disk to be read and displayed by a visualization front-end. If the segmentation result does not satisfy the user, he supplies corrective seeds and the program is rerun. This basic mechanism is sufficient for developing and evaluating methods, but not for an interactive segmentation application integrating the user gesture and the segmentation algorithm. Maximum flow algorithms of the Augmenting-paths variant allow a minimum cut to be recomputed efficiently after modifying relatively few edge capacities [KT]. This property can be taken advantage of in designing an interactive application which allows efficient editing of a current segmentation with the addition or the removal of seeds for correction and refinement purposes. Naturally, parallel maximum flow algorithms would help attain interactive execution rates in such applications.

. Long-term objectives

Such objectives require further research and investigations which entail redesigning algorithms, extending prior models and revisiting the method's theoretical foundations.

Learning prior models. As mentioned previously, adjacency graphs used to derive the vicinity prior were defined manually, which is acceptable for the purposes of method evaluation. In practical applications however, automatic definition of graph models is interesting. For example, such models can be readily generated from existing anatomical phantoms by structure common-surface analysis. In general, automatic definition relies on a learning support of adjacency patterns from presegmented training data, possibly also capturing model variabilities and relation uncertainty, especially with partial data. While graph models of continuous edge weight can represent "degrees" of adjacency, care should be taken for the prior ultimately defined from learned relations  to satisfy submodularity in order for Graph Cut optimization to apply. Priors defined as shortest-path constraints, such as the vicinity prior, do respect this condition, however there remains to be seen whether the formulation is meaningful for variable or different types of relations.

Introducing other relations. It is no doubt interesting to define graph models representing object relationships other than adjacency, such as distance and spatial directionality relationships of the type "object a is on the left of object b". Distance relationships in particular would require a more precise description of the locations and the surfaces of sought objects, therefore we think that prior formulations taking account of object distance would require a learning support as well. Spatial directionality relations rely on orientation information and would call for directed graph structures. The shortest-path formulation on a directed graph however produces asymmetric pairwise constraints which would require revisiting some of the fundamental properties of the underlying theory of Bayesian labeling, in particular, to express the prior distribution according to a directed Markov model [Lau].

Are vicinity prior-based energies tractable? It is known that the shortest-path length between two vertices defines a natural distance metric on graphs bringing notions of convex functions and convex sets to graphs [BC]. Therefore, it is alluring to see whether tools of convex analysis could be applied to energies defined on shortest-path constraints for designing optimizers having better runtime, capable of finding better solutions than the Expansion Moves algorithm. We develop a generic Graph Cut-based semiautomatic multiobject image segmentation method principally for use in routine medical applications ranging from tasks involving few objects in D images to fairly complex near whole-body D image segmentation. The flexible formulation of the method allows its straightforward adaption to a given application.



In particular, the graph-based vicinity prior model we propose, defined as shortest-path pairwise constraints on the object adjacency graph, can be easily reformulated to account for the spatial relationships between objects in a given problem instance. The segmentation algorithm can be tailored to the runtime requirements of the application and the online storage capacities of the computing platform by an efficient and controllable Voronoi tessellation clustering of the input image which achieves a good balance between cluster compactness and boundary adherence criteria. Qualitative and quantitative comprehensive evaluation and comparison with the standard Potts model confirm that the vicinity prior model brings significant improvements in the correct segmentation of distinct objects of identical intensity, the accurate placement of object boundaries and the robustness of segmentation with respect to clustering resolution. Comparative evaluation of the clustering method with competing ones confirms its benefits in terms of runtime and quality of produced partitions. Importantly, compared to voxel segmentation, the clustering step improves both overall runtime and memory footprint of the segmentation process up to an order of magnitude virtually without compromising the segmentation quality.

  Figure .: Anatomical representation between the past and the present. (a) A dissectional drawing of internal organs and genitals from [Ves], (b) a D dissectional view of the Voxel-Man inner organs virtual model [PHP + ] constructed from the Visible Human male image dataset [SASW].

Figure

  Figure .: Commonly encountered artifacts in CT and MR images. (a) (left) Partial volume effect on the clavicle and (right) motion artifacts on lung and esophagus boundaries in a thoracic-abdominal CT scan. (b) Intensity inhomogeneity field affecting muscular and adipose tissues in an MR image of the shoulder. Image courtesy of David Sarrut of Léon Berard cancer research center and Jean-Baptiste Pialat of Édouard Herriot hospital, Lyon, France.

  Figure .: Three commonly encountered user constraints applied to the semiautomatic segmentation of femoral head and hip bones. (a) Object-boundary constraints with part of the boundary traced by the algorithm, (b) object-region constraints for object (blue) and background (red), (c) image-region constraints.

Figure

  Figure .: Illustration of concepts of flow networks after [CLRS]. (a) A flow network with |f | =  indicating the flow and the capacity of each edge, negative flows are not shown, (b) the residual network induced by the previous flow with an augmenting path of capacity c f (p) =  highlighted, (c) the flow network after augmentation through the highlighted path producing the indicated minimum cut with capacity c(S, T ) = .

  the complete residual graph of the flow network in Figure .a is given in Figure .b.

  Figure .c shows the minimum cut of the earlier flow network example having a capacity of c(S, T ) = . Theorem .. (Max-flow Min-cut theorem). Given a flow f in a flow network G = V , E , the following conditions are equivalent:

  giving an initial flow of |f | = , at each iteration, we find an augmenting path p in the residual network G f induced by the current flow f in G, and for each edge (u, v) in p we  A simple path in a directed graph is one the nodes of which are distinct.  increase the flow f (u, v) in G by c f (p) and set the flow on the reverse edge f (v, u) to the negative of the increased flow -f (u, v). After each iteration, the value of flow in G equals |f | + c f (p). When no augmenting paths exist in G f , this flow is a maximum flow. Figure . illustrates the last steps of the method's application. Given the network in Figure .a with the flow |f | = , the path p = (s, v  , v  , t) with c f (p) =  is found in the residual network G f given in Figure .b and the flow in G is updated along the edges of p giving a new flow of |f | =  in Figure .c, which is the maximum flow.

  Figure .: Number of flow network edges to number of nodes ratio for cluster graphs of  images, calculated at  increasing resolutions with % corresponding to the voxel grid graph. Note that edges connecting nodes to flow network terminals are not taken into account as their number is always |V |.

Figure

  Figure .: A labeling instance with L = {"black", "gray", "white"}. (a) The initial configuration, (b) a standard move on the white label, (c) an expansion move on the white label.

  a, b ∈ L and * is the global minimum of E. Note that the algorithm has the best bound on local minima for the Potts model c = .

Figure

  Figure .: (a) Our Synthetic phantom and (b) the graph model of adjacency relations on its structures.

Figure

  Figure .: (a) -label structure adjacency model, (b) geometrical representation of the adjacency prior defined on the model, (c) the model after shadowing the label m, (d) geometrical representation of the extended adjacency prior on the new model. Note that trapezoid side lengths are set for illustration clarity and do not precisely correspond to indicated quantities.

Figure

  Figure .: (a) -label structure adjacency model, (b) the model after shadowing the label m and (c) after shadowing the label r , (d) geometrical representation of the extended adjacency prior on the new model. Note that line segment lengths are set for illustration clarity and do not precisely correspond to indicated quantities. Segments lr and mx are drawn curved for the same purpose.

  Figure .a. Its definition is summarized in Figure .d. Note that for a ratio of two, γ hard =  γ sof t , the extended definition gives the original (.). In particular, vertices m and m coincide and mm = , and the diagonals of the trapezoid lrm m in Figure .d coincide with the sides lm and rm of the triangle lrm in Figure .b with lm = rm = γ sof t . Now we shall consider a more elaborate four-structure model, illustrated in Figure .a, defined on L = {l, m, r, x} by the addition of a label x adjacent to r to the previous three-label model in Figure .a. The adjacency prior definition for this model follows (.), with the following addition:

  a) = γ hard . (.) Once again, we assume that γ hard >  γ sof t . Let us consider the triple {l, m, r} first. Shadowing the label m gives the five-label construction shown in Figure .b. The definition of the extended adjacency prior for the set {l, m, m , r} follows the derivation discussed above which is shown in Figure .d. Furthermore, triples involving labels {l, m, r, x} satisfy submodularity inequalities. However, the condition does not hold for the triple m rx, as expected. As before, we shadow the label r producing the six-label model illustrated in Figure .c. The geometrical representation of the extended adjacency prior for this model is given in Figure .d. Definitions of trapezoids lrm m and xm rr , each corresponding to an extended adjacency prior on a three-structure model, are identical to that in Figure .d. There remains to derive bounds on lr , xm, lx and mr .

  Figure .a. Figure .d summarizes the definition.

  Figure .: Qualitative and quantitative comparison of best brain tissue segmentations, as measured by overall DSM, in a second-generation BrainWeb simulated MR image. Sagittal cross-sections show segmentations using  different prior formulations: (a) Potts, (b) adjacency prior (.) with γ hard =  γ sof t and arbitrarily-penalization submodular adjacency prior for the -structure model with (c) γ h =  γ s , (e) γ h =  γ s and (f) γ h =  γ s . (d) Shows the ground-truth. Color legend; black: background, red: cerebrospinal fluid, green: gray matter, blue: white matter.

Proof.

  By definition, V i,j (a, a) = (a, a) = . Substituting (.) in (.), we try to prove (b, c) ≤ (b, a) + (a, c), ∀a, b, c ∈ L.

  ∃π bac = (b, . . . , a, . . . , c) hence b π bac

Figure

  Figure .: A binary image of a solid circle (a) and  different CVT clusterings corresponding to different parameter settings. (b) Uniform geometric clustering with no regard to intensity or gradient, (c) intensityadaptive clustering with constant density, (d) gradient-adaptive clustering, and (e) intensity and gradientadaptive clustering. For each clustering, the dual graph of the CVT is overlaid on the cluster map.

Figure

  Figure .: The three voxel assignment possibilities to a pair of clusters, in gray and blue, during energy minimization. (a) Initial configuration, (b) C i grows; C i ← u, v , and (c) C j grows; C j ← u, v.

  Algorithm ..: Image clustering via CVT approximation input : image I , number of clusters n, maximum iterations maxIterations output: image partition C = {C  , . . . , C n } indexed by S // Calculate initial clustering

  Figure .: Qualitative and quantitative comparison of -cluster tessellations of a Berkeley dataset grayscale image via the CVT clustering algorithm with (a) m =  and (b) m = , (c) the TurboPixels and, (d) the Graph Cut superpixels algorithms.

  Figure .: Quantitative evaluation of CVT clustering, TurboPixels and Graph Cut superpixels algorithms on  grayscale images from the Berkeley segmentation dataset.

  Figure .: Quantitative evaluation of the CVT clustering algorithm on  grayscale images from the Berkeley segmentation dataset with respect to gradient-based and constant density functions.

  Figure .: Qualitative and quantitative evaluation of the CVT clustering algorithm on a textured image from the Berkeley segmentation dataset with respect to (b) gradient-based and (c) constant density functions. Ground-truth segmentation boundary is given in (a).

  Figure ., which we define in the following. Denote the boundary (surface, in D) of a cluster C i by ∂C i , and the length (area) of the common boundary (surface) of a pair of clusters i, j ∈ S by |∂C i ∩ ∂C j |. Therefore, for a pair of adjacent clusters |∂C i ∩ ∂C j | . Define the dual graph of the CVT, G = S, Ẽ , as an undirected graph on cluster centroids where pairs of adjacent clusters define the set of edges Ẽ = {i, j} | i, j ∈ S, |∂C i ∩ ∂C j |  . Consequently, the neighborhood of a node i ∈ S is defined as Ñi = j | j ∈ S, ∃ {i, j} ∈ Ẽ where i Ñi .

F

  Figure .: Identical binary segmentations on two CVT clusterings of different resolution.

Figure

  Figure .: (a) Our synthetic phantom with structures numbered. (b) A coronal cross-section of a CT image showing corresponding structures. (c) The Shepp-Logan phantom.

Figure

  Figure .: Qualitative inspection and comparison of Potts and vicinity prior pixel segmentations on the synthetic phantom dataset. (top row) The phantom along with "seeds" used in all evaluations on the dataset and three images generated from the phantom with additive zero-mean Gaussian noise of σ = ., . and .. Columns left and right give best Potts and vicinity prior segmentations at a fixed temperature and over all temperature settings on all images generated at a given noise level. Refer to color code in Figure .b.

Figure

  Figure .: Quantitative evaluation and comparison of Potts, adjacency and vicinity prior-based clustered image segmentation on the synthetic phantom dataset.

  Figure .: Quantitative evaluation and comparison of pairwise prior energy weighting schemes in vicinity prior segmentation for  clustering resolutions.

  Figure .: A sagittal cross-section of the head and brain phantom of the original BrainWeb simulated MRI dataset showing artifacts and anatomical inconsistencies.

  Figure .a. Reported clustered image segmentation results correspond to a clustering resolution of %. Vicinity prior performance measures indicating statistically significant improvements at levels α = . and . are highlighted with blue and red respectively. It is easy to observe in these tables that Graph Cut segmentation attains better average performance levels with the vicinity prior on the original BrainWeb database. Moreover, it performs better on clusters, compared to its performance on image voxels. et al. on DSM

Figure

  Figure .: Graph-based anatomical models of structure adjacency relations for (a) the original and (b) the second generation BrainWeb simulated MR dataset phantoms. Structures corresponding to graph vertices; AIR: background, CNT: connective tissue, CSF: cerebrospinal fluid, FAT: fat, GM: gray matter, MRW: bone marrow, MSC: muscle, SKL: skull, SKN: skin, VSL: vessels and WM: white matter.

  Figure .: Surface error measurements for (a) Potts and (b) vicinity prior segmentation applied to the %-clustered second generation BrainWeb simulated MR image, and (c) vicinity prior segmentation applied to the voxels of the image.

  and . respectively on the  principal cross-sections of ground-truth, Potts prior and vicinity prior-segmented images. Disagreements between ground-truth and Potts prior segmentations are encircled. Figure . gives D views of simplified surface meshes generated from different brain and head tissue volumes of the vicinity prior-segmented image via the mesh simplification method outlined in[VC]. Runtime and memory footprint figures are reported in Table.. Lastly, we give the surface error measurements of the aforementioned Potts and vicinity prior-based voxel and clustered image segmentations in Figure . for the  brain tissues. Cumulative error histograms presented in these graphs indicate the fraction of segmented surface lying within a given distance from the true surface. For example, we can see in Figure .b that the vicinity prior-based segmentation applied to the % clustered image retrieves more than % of brain surfaces with an error of  mm, whereas with Potts prior, as shown in Figure .a, the error for the same percentage is well above  mm. Moreover, we notice that surface error measurements of vicinity prior segmentation applied to image voxels, given in Figure .c, are very close to the aforementioned measurements of clustered image segmentation with the same prior.

FigureFigureFigure

  Figure .: Qualitative comparison of Potts and vicinity prior segmentations of the %-clustered second generation BrainWeb simulated MR image. Axial (left column), sagittal (middle column) and coronal (right column) cross-sections correspond to (a)-(c) simulated MR, (d)-(f) ground-truth, (g)-(i) Potts prior and (j)-(l) vicinity prior segmented images. Ellipses on (g)-(i) indicate erroneous segmentations with Potts prior. Figure (b) shows seeds introduced by the user in the corresponding sagittal cross-section.

  Figure .a, are introduced by a user in the axial view of a graphical interface by marking few dozens of pixels from any visible target structure in approximately every tenth volume slice. The vicinity prior is defined according to the graph model of  structures illustrated in Figure ..

  Figure .h, where both lungs are segmented as the left lung. Corresponding error graphs in Figure .a indicate lower surface error for the right lung whereas its DSM entry in Table

Figure

  Figure .: Qualitative comparison of Potts and vicinity prior segmentations of the %-clustered D-IRCADb thoracic-abdominal CT image. Axial (left column), coronal (middle column) and sagittal (right column) cross-sections correspond to (a)-(c) CT, (d)-(f) ground-truth, (g)-(i) Potts prior and (j)-(l) vicinity prior segmented images. Refer to the color code in Figure . or the color legend in Figure .. Ellipses on (g)-(i) indicate erroneous segmentations with Potts prior. Notice how right and left lungs were correctly segmented as distinct organs with vicinity prior. Figure (a) shows seeds introduced by the user in the corresponding axial cross-section.

Figure

  Figure .: Images from the Visible Human female dataset. (a) a cryosectional photo at the heart level and (b) the corresponding CT slice.

Figure

  Figure .: (a) An axial cross-section of femoral head and hip bones from the Visible Human male dataset, bones segmented (b) with and (c) without contrast function weighted Potts prior.



  Figure .: Lateral ventral and dorsal D views of surface mesh generated from Potts prior skeleton segmentation of a %-clustered reconstructed CT image from the Visible Human female dataset.

  Figure .: Performance evaluation for segmentation program total and individual step runtimes measured on a  ×  pixel image from the synthetic phantom dataset generated at σ = .. Figures are recorded on a desktop computer powered by a quad-core processor running at . GHz with  GB of RAM.

  de Structure pour la Segmentation Multi-objet d'Images Médicales D par Partition d'Images et Coupure de Graphes NATURE: Doctorat Numéro d'ordre: XXXXXX École doctorale: EEA Spécialité: Traitement d'images RÉSUMÉ: Nous développons une méthode générique semi-automatique multi-objet de segmentation d'image par coupure de graphe visant les usages médicaux de routine, allant des tâches impliquant quelques objets dans des images D, à quelques dizaines dans celles D quasi corps entier. La formulation souple de la méthode permet son adaptation simple à une application donnée. En particulier, le modèle d'a priori de proximité que nous proposons, défini à partir des contraintes de paires du plus court chemin sur le graphe d'adjacence des objets, peut facilement être adapté pour tenir compte des relations spatiales entre les objets ciblés dans un problème donné. L'algorithme de segmentation peut être adapté aux besoins de l'application en termes de temps d'exécution et de capacité de stockage à l'aide d'une partition de l'image à segmenter par une tesselation de Voronoï efficace et contrôlable, établissant un bon équilibre entre la compacité des régions et le respect des frontières des objets. Des évaluations et comparaisons qualitatives et quantitatives avec le modèle de Potts standard confirment que notre modèle d'a priori apporte des améliorations significatives dans la segmentation d'objets distincts d'intensités similaires, dans le positionnement précis des frontières des objets ainsi que dans la robustesse de segmentation par rapport à la résolution de partition. L'évaluation comparative de la méthode de partition avec ses concurrentes confirme ses avantages en termes de temps d'exécution et de qualité des partitions produites. Par comparaison avec l'approche appliquée directement sur les voxels de l'image, l'étape de partition améliore à la fois le temps d'exécution global et l'empreinte mémoire du processus de segmentation jusqu'à un ordre de grandeur, sans compromettre la qualité de la segmentation en pratique. MOTS-CLÉS: segmentation d'images, partition d'images, modèle de Markov, a priori spatial, coupure de graphe. Laboratoire de recherche: CREATIS; INSA de Lyon, GIPSA-lab; Grenoble INP. Directeurs de thèse: Rémy PROST, Michel DESVIGNES.

  Keywords: image segmentation, image clustering, Markov model, spatial prior, Graph Cut A Priori de Structure pour la Segmentation Multi-objet d'Images Médicales D par Partition d'Images et Coupure de Graphes Résumé:Nous développons une méthode générique semi-automatique multi-objet de segmentation d'image par coupure de graphe visant les usages médicaux de routine, allant des tâches impliquant quelques objets dans des images D, à quelques dizaines dans celles D quasi corps entier. La formulation souple de la méthode permet son adaptation simple à une application donnée. En particulier, le modèle d'a priori de proximité que nous proposons, défini à partir des contraintes de paires du plus court chemin sur le graphe d'adjacence des objets, peut facilement être adapté pour tenir compte des relations spatiales entre les objets ciblés dans un problème donné. L'algorithme de segmentation peut être adapté aux besoins de l'application en termes de temps d'exécution et de capacité de stockage à l'aide d'une partition de l'image à segmenter par une tesselation de Voronoï efficace et contrôlable, établissant un bon équilibre entre la compacité des régions et le respect des frontières des objets. Des évaluations et comparaisons qualitatives et quantitatives avec le modèle de Potts standard confirment que notre modèle d'a priori apporte des améliorations significatives dans la segmentation d'objets distincts d'intensités similaires, dans le positionnement précis des frontières des objets ainsi que dans la robustesse de segmentation par rapport à la résolution de partition. L'évaluation comparative de la méthode de partition avec ses concurrentes confirme ses avantages en termes de temps d'exécution et de qualité des partitions produites. Par comparaison avec l'approche appliquée directement sur les voxels de l'image, l'étape de partition améliore à la fois le temps d'exécution global et l'empreinte mémoire du processus de segmentation jusqu'à un ordre de grandeur, sans compromettre la qualité de la segmentation en pratique. Mots-Clés: segmentation d'images, partition d'images, modèle de Markov, a priori spatial, coupure de graphe.CREATIS: Centre de Recherche en Acquisition et Traitement de l'Image pour la Santé Université de Lyon; CNRS UMR ; Inserm U; INSA-Lyon; Université Lyon   avenue Jean Capelle,  Villeurbanne, France

Table 

 

.: Factorization of a pairwise function of binary variables into a constant, two unary functions and a pairwise function respectively.

  Compared to the Potts prior approach, im-Quantitative evaluation and comparison of the robustness of Potts, adjacency and vicinity prior-based image segmentation to clustering resolution on the synthetic phantom dataset.

			Robustness to clustering resolution at noise level σ = .					Robustness to clustering resolution at noise level σ = .					Robustness to clustering resolution at noise level σ = .
																																																						vicinity
																																																						adjacency
	DSM	. .																			DSM	. .																				DSM	. .												Potts
	Mean	.																			Mean	.																				Mean	.
		.																				.																					.
					vicinity																		vicinity
		.			adjacency														.				adjacency														.
						Potts																			Potts
		.																				.																					.
																																																														
			Number of clusters as % of voxel count				Number of clusters as % of voxel count				Number of clusters as % of voxel count
										(a)																			(b)																			(c)
	. . . . .  Figure .: . Mean DSM Prior weighting schemes at % clustering common surface inv. distance both schemes no weighting	Mean DSM	. . . . .  .	Prior weighting schemes at % clustering common surface inv. distance both schemes no weighting	Mean DSM	. . . . .  .	Prior weighting schemes at % clustering common surface inv. distance both schemes no weighting
		 . . . . . . . . . .			 . . . . . . . . . .			 . . . . . . . . . .
						Noise standard deviation										Noise standard deviation									Noise standard deviation
																																			

Table  .

  . . ± . . ± . . ± . Cerebrospinal fluid . ± . . ± . . ± . . ± . Gray matter . ± . . ± . . ± . . ± . White matter. ± . . ± . . ± . . ± . Fat . ± . . ± . . ± . . ± . Muscle . ± . . ± . . ± . . ± . Skin . ± . . ± . . ± . . ± . Skull . ± . . ± . . ± . . ± . Connective tissue . ± . . ± . . ± . . ± . Overall . ± . . ± . . ± . . ± . . . ± . . ± . . ± . Cerebrospinal fluid . ± . . ± . . ± . . ± . Gray matter . ± . . ± . . ± . . ± . White matter . ± . . ± . . ± . . ± . Fat . ± . . ± . . ± . . ± . Muscle . ± . . ± . . ± . . ± . Skin . ± . . ± . . ± . . ± . Skull . ± . . ± . . ± . . ± . Connective tissue . ± . . ± . . ± . . ± . . . ± . . ± . . ± . Cerebrospinal fluid . ± . . ± . . ± . . ± . Gray matter . ± . . ± . . ± . . ± . White matter . ± . . ± . . ± . . ± . Fat . ± . . ± . . ± . . ± . Muscle . ± . . ± . . ± . . ± . Skin . ± . . ± . . ± . . ± . Skull . ± . . ± . . ± . . ± . Connective tissue. ± . . ± . . ± . . ± . : Means and standard deviations of overall and individual-class DSM, sensitivity and specificity measures for voxel and %-clustered image segmentations of the original BrainWeb dataset with Potts and vicinity priors. Significant improvements at α = . and . are given in blue and red respectively.

			DSM
			voxels	clusters %
	Structures	Potts	vicinity	Potts	vicinity
	Air	. ± (a)
			Sensitivity
			voxels	clusters %
	Structures	Potts	vicinity	Potts	vicinity
	Air . ± Overall . ± . . ± . . ± . . ± .
			(b)
			Specificity
			voxels	clusters %
	Structures	Potts	Vicinity	Potts	Vicinity
	Air . ± Overall . ± . . ± . . ± . . ± .
			(c)

  The motivation behind keeping the temperature at a "factory default" value of . in the present as well as in other single-image case studies in the following sections comes from observing that in Figures .a and .b Graph Cut segmentation with both priors, and especially with vicinity prior, is reasonably robust to temperature settings in the range [., .] for low to moderate noise levels [, .]. Seeds, illustrated in Figure .b, are introduced by a user in the sagittal view of a graphical interface by marking about  pixels for every structure via brush strokes of size  ×  pixels at most distributed evenly in the volume. The vicinity prior is defined according to the graph model given in Figure .b.  mm new BrainWeb % clustered seg. with Potts prior new BrainWeb voxel seg. with vicinity prior new BrainWeb % clustered seg. with vicinity prior

			Cumulative error histogram for			Cumulative error histogram for			Cumulative error histogram for
	fraction w.r.t. ground-truth	. . . . 		fraction w.r.t. ground-truth	. . . . 		fraction w.r.t. ground-truth	. . . . 
	Surface	. .	CSF, max. err. . mm GM, max. err. . mm	Surface	. .	CSF, max. err. . mm GM, max. err. . mm	Surface	. .	CSF, max. err.  mm GM, max. err. . mm
		.	WM, max. err. . mm		.	WM, max. err. . mm		.	WM, max. err.
			               			               			               
			Error (in millimeters)			Error (in millimeters)			Error (in millimeters)
						

Table  .

  : Overall and individual-class performance measures for (a) %-clustered and (b) voxel segmentations of the new BrainWeb simulated MR image with Potts and vicinity priors. Significant improvements at α = . are given in red.

		DSM	Sensitivity	Specificity
	Structures	Potts vicinity Potts vicinity Potts vicinity
	Air	.	. .	. .	.
	Cerebrospinal fluid .	. .	. .	.
	Gray matter	.	. .	. .	.
	White matter	.	. .	. .	.
	Fat	.	. .	. .	.
	Muscle	.	. .	. .	.
	Skin	.	. .	. .	.
	Skull	.	. .	. .	.
	Vessels	.	. .	. .	.
	Marrow	.	. .	. .	.
	mean	.	. .	. .	.
	Overall	.	. .	. .	.
			(a)
		DSM	Sensitivity	Specificity
	Structures	Potts vicinity Potts vicinity Potts vicinity
	Air	.	. .	. .	.
	Cerebrospinal fluid .	. .	. .	.
	Gray matter	.	. .	. .	.
	White matter	.	. .	. .	.
	Fat	.	. .	. .	.
	Muscle	.	. .	. .	.
	Skin	.	. .	. .	.
	Skull	.	. .	. .	.
	Vessels	.	. .	. .	.
	Marrow	.	. .	. .	.
	mean	.	. .	. .	.
	Overall	.	. .	. .	.
			(b)
	Image + size (voxels)	Segmentation Prior Clust. (m) Constr. (m) Optim. (m) Mem. (Mb)
	BrainWeb subject no. 	voxels	Potts vicinity	-	. .		. .	. .
	 ×  × 	clusters %	Potts vicinity	.	. .		. .	. .
		

Table  .

  : Performance figures for Potts and vicinity prior segmentation applied to the voxels and a % clustering of the second generation BrainWeb simulated MR image presented in Figures . and . respectively. Last four columns indicate execution time in minutes for the clustering algorithm, graph construction, Expansion Moves optimization and the program memory footprint in megabytes. Figures are recorded ona desktop computer powered by a quad-core processor running at . GHz with  GB of RAM.

 . Real CT images .. Thoracic-abdominal CT image

Table  .

  : Performance figures for Potts and vicinity prior segmentation applied to the % and %clustered D-IRCADb thoracic-abdominal CT image. Last four columns indicate execution time in minutes for the clustering algorithm, graph construction, Expansion Moves optimization and the program memory footprint in megabytes. Figures are recorded ona desktop computer powered by a quad-core processor running at . GHz with  GB of RAM. . adrenal gland (L) . hyperplasia . aorta . vena cava . portal vein . skin . fat . visceral fat . muscle . marrow . renal pyramids (R) . renal pyramids (L) . digestive tract tissue . digestive tract cavity . pulmonary arteries (R) . pulmonary arteries (L) . bronchi (R) . bronchi (L) . intervertebral discs . hepatic arteries . hepatic vena cava . hepatic portal veins.

				1
			35	22
	26	3	23
			15	25	16
	27	13		14	28
				24
					29	30
	19	20 21		10
	36 37 38
	9			4
		7			8
	18		31	32
		5			6
			33	34
	Image + size (voxels)	Segmentation Prior Clust. (m) Constr. (m) Optim. (m) Mem. (Mb)
	D-IRCADb CT	clusters %	Potts vicinity	.	. .	. .	. .
	 ×  × 	clusters %	Potts vicinity	.	. .	. .	. .
				

Figure .: Graph-based anatomical model of structure adjacency relations for the D-IRCADb thoracicabdominal CT image. Structures corresponding to graph vertices: . air . bone . heart . lung (R) . lung (L) . liver . spleen . gallbladder . pancreas . kidney (R) . kidney (L) . adrenal gland (R)

Table  .

  : Overall and individual-class performance measures for (a) %-clustered and (b) %-clustered segmentations of D-IRCADb thoracic-abdominal CT image with Potts and vicinity priors. Significant improvements at α = . and . are given in blue and red respectively.

		DSM	Sensitivity	Specificity
	Structures	Potts vicinity Potts vicinity Potts vicinity
	Bone	.	. .	. .	.
	Heart	.	. .	. .	.
	Lung (R)	.	. .	. .	.
	Lung (L)	.	. .	. .	.
	Liver	.	. .	. .	.
	Spleen	.	. .	. .	.
	Gallbladder	.	. .	. .	.
	Pancreas	.	. .	. .	.
	Kidney (R)	.	. .	. .	.
	Kidney (L)	.	. .	. .	.
	Adrenal gland (R) .	. .	. .	.
	Adrenal gland (L) .	. .	. .	.
	Hyperplasia	.	. .	. .	.
	Aorta	.	. .	. .	.
	Vena Cava	.	. .	. .	.
	Portal vein	.	. .	. .	.
	mean	.	. .	. .	.
	Overall	.	. .	. .	.
			(a)
		DSM	Sensitivity	Specificity
	Structures	Potts vicinity Potts vicinity Potts vicinity
	Bone	.	. .	. .	.
	Heart	.	. .	. .	.
	Lung (R)	.	. .	. .	.
	Lung (L)	.	. .	. .	.
	Liver	.	. .	. .	.
	Spleen	.	. .	. .	.
	Gallbladder	.	. .	. .	.
	Pancreas	.	. .	. .	.
	Kidney (R)	.	. .	. .	.
	Kidney (L)	.	. .	. .	.
	Adrenal gland (R) .	. .	. .	.
	Adrenal gland (L) .	. .	. .	.
	Hyperplasia	.	. .	. .	.
	Aorta	.	. .	. .	.
	Vena Cava	.	. .	. .	.
	Portal vein	.	. .	. .	.
	mean	.	. .	. .	.
	Overall	.	. .	. .	.
			(b)
				

  Two views of surface meshes generated via the mesh simplification method of[VC] from a skeleton segmentation of the %-clustered reconstructed D CT volume are given in Figure .. Runtime and memory footprint measures are given in Table.. Due to insufficient memory resources, we are unable to run our segmentation program on image voxels to draw comparisons with the clustered variant.

	edge	weight
	(i, j)

Table  .

  : Definitions of maximum flow graph edge weights for Visible Human female CT skeleton segmentation.

	Image + size (voxels) Segmentation Clust. (m) Constr. (m) Optim. (m) Mem. (Mb)
	VH CT female	voxels	-	N/A	N/A	N/A
	 ×  × 	clusters %	.	.	.	.

Table  .

  : Performance figures for Potts prior segmentation applied to %-clustered reconstructed CT image of the Visible Human female dataset presented in Figure .. Last four columns indicate execution time in minutes for the clustering algorithm, graph construction, Maximum flow optimization and the program memory footprint in megabytes. Figures are recorded on a dedicated cluster computer powered by a -core processor running at . GHz with  GB of RAM. Due to insufficient memory resources we are unable to apply segmentation to the voxels of the image.

  Structural Priors for Multiobject Semiautomatic Segmentation of Three-dimensional Medical Images via Clustering and Graph Cut Algorithms Abstract:

	Président de jury:
	Composition du jury:
	Philippe BOLON	Professeur, Polytech Annecy-Chambéry	Examinateur
	Fabrice HEITZ	Professeur, Université de Strasbourg	Rapporteur
	Jean-Philippe THIRAN	Professeur, EPFL de Lausanne	Rapporteur
	Sébastien VALETTE	Chargé de recherche, CNRS	Co-encadrant
	Michel DESVIGNES	Professeur, Grenoble INP	Co-directeur
	Rémy PROST	Professeur, INSA de Lyon	Co-directeur

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0058/these.pdf © [R. Kéchichian],[2013], INSA de Lyon, tous droits réservés

Cette thèse est accessible l'adresse : http://theses.insa-lyon.fr/publication/2013ISAL0058/these.pdf © [R. Kéchichian],[2013], INSA de Lyon, tous droits réservés

Acknowledgements

Acknowledgements vii

Table of contents ix



Error (in centimeters) 