This PDF file contains the front matter associated with SPIE Proceedings volume 7623, including the Title Page, Copyright information, Table of Contents, Introduction, and the Conference Committee listing.

Cochlear implant surgery is a procedure performed to treat profound hearing loss. Since the cochlea is not visible in surgery, the physician uses anatomical landmarks to estimate the pose of the cochlea. Research has indicated that implanting the electrode in a particular cavity of the cochlea, the scala tympani, results in better hearing restoration. The success of the scala tympani implantation is largely dependent on the point of entry and angle of electrode insertion. Errors can occur due to the imprecise nature of landmark-based, manual navigation as well as inter-patient variations between scala tympani and the anatomical landmarks. In this work, we use point distribution models of the intra-cochlear anatomy to study the inter-patient variations between the cochlea and the typical anatomic landmarks, and we implement an active shape model technique to automatically localize intra-cochlear anatomy in conventional CT images, where intra-cochlear structures are not visible. This fully automatic segmentation could aid the surgeon to choose the point of entry and angle of approach to maximize the likelihood of scala tympani insertion, resulting in more substantial hearing restoration.

Different imaging modalities or protocols of a single patient may convey different types of information regarding a disease for the same anatomical organ/tissue. On the other hand, multi-modal/multi-protocol medical images from several different patients can also provide spatial statistics of the disease occurrence, which in turn can greatly aid in disease diagnosis and aid in improved, accurate biopsy and targeted treatment. It is therefore important to not only integrate medical images from multiple patients into a common coordinate frame (in the form of a population-based atlas), but also find the correlation between these multi-modal/multi-protocol data features and the disease spatial distribution in order to identify different quantitative structural and functional disease signatures. Most previous work on construction of anatomical atlases has focused on deriving a population-based atlas for the purpose of deriving the spatial statistics. Moreover, these models are typically derived from normal or healthy subjects, either explicitly or implicitly, where it is assumed that the inter-patient pathological variation is not large. These methods are not suitable for constructing a disease atlas, where significant differences between patients on account of disease related variations can be expected. In this paper, we present a novel framework for the construction of a multi-parametric MRI-based data-driven disease atlas consisting of multi-modal and multi-protocol data from across multiple patient studies. Our disease atlas contains 3 Tesla structural (T2) and functional (dynamic contrast enhanced (DCE)) prostate in vivo MRI with corresponding whole mount histology specimens obtained via radical prostatectomy. Our atlas construction framework comprises 3 distinct modules: (a) determination of disease spatial extent on the multi-protocol MR imagery for each patient, (b) construction of a multi-protocol MR imaging spatial atlas which captures the geographical proclivity of the disease, and (c) feature extraction and the construction of the data-driven multi-protocol MRI based prostate cancer atlas. The marriage of data driven and spatial atlases could serve as a useful tool for clinicians to identifying structural and functional imaging disease signatures so as to make better, more informed diagnoses. Each spatial location in this atlas can be associated with a high dimensional multi-attribute quantitative feature vector. Additionally, since the feature vectors are extracted from across multiple patient studies, each spatial location in the data-driven atlas can be characterized by a feature distribution (in turn characterized by a mean and standard deviation). Preliminary investigation in quantitatively correlating the disease signatures from across the spatial and data driven atlases suggests that our quantitative atlas framework could emerge as a powerful tool for discovering prostate cancer imaging signatures.

Distortion of Left Ventricle (LV) external anatomy is related to some dysfunctions, such as hypertrophy. The architecture of myocardial fibers determines LV electromechanical activation patterns as well as mechanics. Thus, their joined modelling would allow the design of specific interventions (such as peacemaker implantation and LV remodelling) and therapies (such as resynchronization). On one hand, accurate modelling of external anatomy requires either a dense sampling or a continuous infinite dimensional approach, which requires non-Euclidean statistics. On the other hand, computation of fiber models requires statistics on Riemannian spaces. Most approaches compute separate statistical models for external anatomy and fibers architecture. In this work we propose a general mathematical framework based on differential geometry concepts for computing a statistical model including, both, external and fiber anatomy. Our framework provides a continuous approach to external anatomy supporting standard statistics. We also provide a straightforward formula for the computation of the Riemannian fiber statistics. We have applied our methodology to the computation of complete anatomical atlas of canine hearts from diffusion tensor studies. The orientation of fibers over the average external geometry agrees with the segmental description of orientations reported in the literature.

Computed tomography angiography (CTA), a non-invasive imaging technique, is becoming increasingly popular for cardiac examination, mainly due to its superior spatial resolution compared to MRI. This imaging modality is currently widely used for the diagnosis of coronary artery disease (CAD) but it is not commonly used for the diagnosis of ventricular and atrial function. In this paper, we present a fully automatic method for segmenting the whole heart (i.e. the outer surface of the myocardium) and cardiac chambers from CTA datasets. Cardiac chamber segmentation is particularly valuable for the extraction of ventricular and atrial functional information, such as stroke volume and ejection fraction. With our approach, we aim to improve the diagnosis of CAD by providing functional information extracted from the same CTA data, thus not requiring additional scanning. In addition, the whole heart segmentation method we propose can be used for visualization of the coronary arteries and for obtaining a region of interest for subsequent segmentation of the coronaries, ventricles and atria. Our approach is based on multi-atlas segmentation, and performed within a non-rigid registration framework. A leave-one-out quantitative validation was carried out on 8 images. The method showed a high accuracy, which is reflected in both a mean segmentation error of 1.05±1.30 mm and an average Dice coefficient of 0.93. The robustness of the method is demonstrated by successfully applying the method to 243 additional datasets, without any significant failure.

Using an atlas, an image can be segmented by mapping its coordinate space to that of the atlas in an anatomically correct way. In order to find the correct mapping between the two different coordinate spaces e.g. diffeomorphic demons registration can be applied. The demons algorithm is a popular choice for deformable image registration and offers the possibility to perform computationally efficient non-rigid (diffeomorphic) registration. However, this registration method is prone to image artifacts and image noise. Therefore it has been the main objective of the presented work to combine the efficiency of diffeomorphic demons and the stability of statistical models. In the presented approach a statistical deformation model that describes "anatomically correct" displacements vector fields for a specific registration problem is used to guide the demons registration algorithm. By projecting the current displacement vector field, which is calculated during any iteration of the registration process, into the model space a regularized version of the vector field can be computed. Using this regularized vector field for the update of the deformation field in the subsequent iteration of the registration process the demons registration algorithm can be guided by the deformation model. The proposed method was evaluated on 21 CT datasets of the right hip. Measuring the average and maximum segmentation error for all 21 datasets and all 120 test configurations it could be demonstrated that the newly proposed algorithm leads to a reduction of the segmentation error of up to 13% compared to using the conventional diffeomorphic demons algorithm.

We present a novel variational framework for deformable multi-modal image registration. Our approach is based on Renyi's statistical dependence measure of two random variables with the use of reproducing kernel Hilbert spaces associated with Gaussian kernels to simplify the computation. The popularly used method of maximizing mutual information based optimization algorithms are complex and sensitive to the quantization of the intensities, because it requires the estimation of continuous joint probability density function (pdf). The proposed model does not deal with joint pdf but instead observed independent samples. Experimental results are provided to show the effectiveness of the model.

An image, being a continuous function, is commonly discretely represented as a set of sample values, namely the intensities, associated with the spatial grids. After that, all types of the operations are then carried out there in. We denote such representation as the discrete function representation (DFR). In this paper we provide another discrete representation for images using the point sets, called the point set representation (PSR). Essentially, the image is normalized to have the unit integral and is treated as a probability density function of some random variable. Then the PSR is formed by drawing samples of such random variable. In contrast with the DFR, here the image is purely represented as points and no values are associated. Besides being an equivalent discrete representation for images, we show that certain image operations benefit from such representation in the numerical stability, performance, or both. Examples are given in the Perona-Malik type diffusion where in the PSR there is no such problem as the numerical instability. Furthermore, PSR naturally bridges the fields of image registration with the point set registration. This helps handle some otherwise difficult problems in image registration such as partial image registration, with much faster convergence speed.

In registration-based analyses of lung biomechanics and function, high quality registrations are essential to obtain meaningful results. Various criteria have been suggested to find the correspondence mappings between two lung images acquired at different levels of inflation. In this paper, we describe a new metric, the sum of squared vesselness measure difference (SSVMD), that utilizes the rich information of blood vessel locations and matches similar vesselness patterns in two images. Preserving both the lung tissue volume and the vesselness measure, a registration algorithm is developed to minimize the sum of squared tissue volume difference (SSTVD) and SSVMD together. We compare the registration accuracy using SSTVD + SSVMD with that using SSTVD alone by registering lung CT images of three normal human subjects. After adding the new SSVMD metric, the improvement of registration accuracy is observed by landmark error and fissure positioning error analyses. The average values of landmark error and fissure positioning error are reduced by about 30% and 25%, respectively. The mean landmark error is on the order of 1 mm. Statistical testing of landmark errors shows that there is a statistically significant difference between two methods with p values < 0.05 in all three subjects. Visual inspection shows there are obvious accuracy improvements in the lung regions near the thoracic cage after adding SSVMD.

In this work we evaluate a novel method for multi-modal image registration of MR images. The key feature of our approach is a new distance measure that allows for comparing modalities that are related by an arbitrary gray-value mapping. The novel measure is formulated as least square problem for minimizing the sum of squared differences of two images with respect to changing gray-values of one of the images. It turns out that the novel measure can be computed explicitly and allows for very simple and efficient implementation. We compare our new approach to rigid registration with cross-correlation, mutual information, and normalized gradient fields as distance measure.

Quadratic regularizers are used in nonparametric registration to ensure that the registration problem is well posed and to yield solutions that exhibit certain types of smoothness. Examples of popular quadratic regularizers include the diffusion, elastic, fluid, and curvature regularizers.¹ Two important features of these regularizers are whether they account for coupling of the spatial components of the deformation (elastic/fluid do; diffusion/curvature do not) and whether they are robust to initial affine misregistrations (curvature is; diffusion/ elastic/fluid are not). In this article, we show how to extend this list of quadratic regularizers to include a second-order regularizer that exhibits the best of both features: it accounts for coupling of the spatial components of the deformation and contains affine transformations in its kernel. We then show how this extended taxonomy of quadratic regularizers is related to other families of regularizers, including Cachier and Ayache's differential quadratic forms2 and Arigovindan's family of rotationally invariant regularizers.^{3, 4} Next, we describe two computationally efficient paradigms for performing nonparametric registration with the proposed regularizer, based on Fourier methods5 and on successive Gaussian convolution.^{6, 7} Finally, we illustrate the performance of the quadratic regularizers on the task of registering serial 3-D CT exams of patients with lung nodules.

A novel approach for coupling brain tumor mass effect with a continuous model of cancer progression is proposed. The purpose of the present work is to devise an efficient approximate model for the mechanical interaction of the tumor with its surroundings in order to aid registration of brain tumor images with statistical atlases as well as the generation of atlases of brain tumor disease. To model tumor progression a deterministic reaction-diffusion formalism, which describes the spatio-temporal dynamics of a coarse-grained population density of cancerous cells, is discretized on a regular grid. Tensor information obtained from a probabilistic atlas is used to model the anisotropy of the diffusion of malignant cells within white matter. To account for the expansive nature of the tumor a parametric deformation model is linked to the computed net cell density of cancerous cells. To this end, we formulate a constrained optimization problem using an inhomogeneous regularization that in turn allows for approximating physical properties of brain tissue. The described coupling model can in general be applied to estimate mass effect of non-convex, diffusive as well as multifocal tumors so that no simplification of the growth model has to be stipulated. The present work has to be considered as a proof-of-concept. Visual assessment of the computed results demonstrates the potential of the described method. We conclude that the analogy to the problem formulation in image registration potentially allows for a sensible integration of the described approach into a unified framework of image registration and tumor modeling.

In this paper we address the problem of 3D shape reconstruction from sparse X-ray projections. We present a correspondence free method to fit a statistical shape model to two X-ray projections, and illustrate its performance in 3D shape reconstruction of the femur. The method alternates between 2D segmentation and 3D shaoe reconstruction, where 2D segmentation is guided by dynamic programming along the model projection on the X-ray plane. 3D reconstruction is based on the iterative minimization of the 3D distance between a set of support points and the back-projected silhouette with respect to the pose and model parameters. We show robustness of the reconstruction on simulated X-ray projection data of the femur, varying the field of view; and in a pilot study on cadaveric femora.

In this work, we present a four-dimensional reconstruction technique for cardiac gated SPECT images using a content-adaptive deformable mesh model. Cardiac gated SPECT images are affected by a high level of noise. Noise reduction methods usually do not account for cardiac motion and therefore introduce motion blur-an artifact that can decrease diagnostic accuracy. Additionally, image reconstruction methods typically rely on uniform sampling and Cartesian griding for image representation. The proposed method utilizes a mesh representation of the images in order to utilize the benefits of content-adaptive nonuniform sampling. The mesh model allows for accurate representation of important regions while significantly compressing the data. The content-adaptive deformable mesh model is generated by combining nodes generated on the full torso using pre-reconstructed emission and attenuation images with nodes accurately sampled on the left ventricle. Ventricular nodes are further displaced according to cardiac motion using our previously introduced motion estimation technique. The resulting mesh structure is then used to perform iterative image reconstruction using a mesh-based maximum-likelihood expectation-maximization algorithm. Finally, motion-compensated post-reconstruction temporal filtering is applied in the mesh domain using the deformable mesh model. Reconstructed images as well as quantitative evaluation show that the proposed method offers improved image quality while reducing the data size.

Several studies based on biplanar radiography technologies are foreseen as great systems for 3D-reconstruction applications for medical diagnoses. This paper proposes a non-rigid registration method to estimate a 3D personalized shape of bone models from two planar x-ray images using an as-rigid-as-possible deformation approach based on a moving least-squares optimization method. Based on interactive deformation methods, the proposed technique has the ability to let a user improve readily and with simplicity a 3D reconstruction which is an important step in clinical applications. Experimental evaluations of six anatomical femur specimens demonstrate good performances of the proposed approach in terms of accuracy and robustness when compared to CT-scan.

The automatic recognition of vascular trees is a challenging task, required for roadmapping or advanced visualization. For instance, during an endovascular aneurysm repair (EVAR), the recognition of abdominal arteries in angiograms can be used to select the appropriate stent graft. This choice is based on a reduced set of arteries (aorta, renal arteries, iliac arteries) whose relative positions are quite stable. We propose in this article a recognition process based on a structural model. The centerlines of the target vessels are represented by a set of control points whose relative positions are constrained. To find their position in an angiogram, we enhance the target vessels and extract a set of possible positions for each control point. Then, a constraint propagation algorithm based on the model prunes those sets of candidates, removing inconsistent ones. We present preliminary results on 5 cases, illustrating the potential of this approach and especially its ability to handle the high variability of the target vessels.

Accurate definition of the aorta and pulmonary artery from three-dimensional (3D) multi-detector CT (MDCT) images is important for pulmonary applications. This work presents robust methods for defining the aorta and pulmonary artery in the central chest. The methods work on both contrast enhanced and no-contrast 3D MDCT image data. The automatic methods use a common approach employing model fitting and selection and adaptive refinement. During the occasional event that more precise vascular extraction is desired or the method fails, we also have an alternate semi-automatic fail-safe method. The semi-automatic method extracts the vasculature by extending the medial axes into a user-guided direction. A ground-truth study over a series of 40 human 3D MDCT images demonstrates the efficacy, accuracy, robustness, and efficiency of the methods.

In cardiac MR images the slice thickness is normally greater than the pixel size within the slices. In general, better segmentation and analysis results can be expected for isotropic high-resolution (HR) data sets. If two orthogonal data sets, e. g. short-axis (SA) and long-axis (LA) volumes are combined, an increase in resolution can be obtained. In this work we employ a super-resolution reconstruction (SRR) algorithm for computing high-resolution data sets from two orthogonal SA and LA volumes. In contrast to a simple averaging of both data in the overlapping region, we apply a maximum a posteriori approach. There, an observation model is employed for estimating an HR image that best reproduces the two low-resolution input data sets. For testing the SRR approach, we use clinical MRI data with an in-plane resolution of 1.5 mm×1.5 mm and a slice thickness of 8 mm. We show that the results obtained with our approach are superior to currently used averaging techniques. Due to the fact that the heart deforms over the cardiac cycle, we investigate further, how the replacement of a rigid registration by a deformable registration as preprocessing step improves the quality of the final HR image data. We conclude that image quality is dramatically enhanced by applying an SRR technique especially for cardiac MR images where the resolution in slice-selection direction is about five times lower than within the slices.

Tissue contrast and resolution of magnetic resonance neuroimaging data have strong impacts on the utility of the data in clinical and neuroscience tasks such as registration and segmentation. Lengthy acquisition times typically prevent routine acquisition of multiple MR tissue contrast images at high resolution, and the opportunity for detailed analysis using these data would seem to be irrevocably lost. This paper describes an example based approach using patch matching from a multiple resolution multiple contrast atlas in order to change an image's resolution as well as its MR tissue contrast from one pulse-sequence to that of another. The use of this approach to generate different tissue contrasts (T2/PD/FLAIR) from a single T1-weighted image is demonstrated on both phantom and real images.

A wide range of medical applications in clinic and research exploit images acquired by fast magnetic resonance imaging (MRI) sequences such as echo-planar imaging (EPI), e.g. functional MRI (fMRI) and diffusion tensor MRI (DT-MRI). Since the underlying assumption of homogeneous static fields fails to hold in practical applications, images acquired by those sequences suffer from distortions in both geometry and intensity. In the present paper we propose a new variational image registration approach to correct those EPI distortions. To this end we acquire two reference EPI images without diffusion sensitizing and with inverted phase encoding gradients in order to calculate a rectified image. The idea is to apply a specialized registration scheme which compensates for the characteristical direction dependent image distortions. In addition the proposed scheme automatically corrects for intensity distortions. This is done by evoking a problem dependent distance measure incorporated into a variational setting. We adjust not only the image volumes but also the phase encoding direction after correcting for patients head-movements between the acquisitions. Finally, we present first successful results of the new algorithm for the registration of DT-MRI datasets.

Arterial spin labeling (ASL) is a noninvasive MRI method that uses magnetically labeled blood to measure cerebral perfusion. Spatial resolution of ASL is relatively small and as a consequence perfusion from different tissue types is mixed in each pixel. An average ratio of gray matter (GM) to white matter (WM) blood flow is 3.2 to 1. Disregarding the partial volume effects (PVE) can thus cause serious errors of perfusion quantification. PVE also complicates spatial filtering of ASL images as apart from noise there is a spatial signal variation due to tissue partial volume. Recently, an algorithm for correcting PVE has been published by Asllani et al. It represents the measured magnetization as a sum of different tissue magnetizations weighted by their fractional volume in a pixel. With the knowledge of the partial volume obtained from a high-resolution MRI image, it is possible to separate the individual tissue contributions by linear regression on a neighborhood of each pixel. We propose an extension of this algorithm by minimizing the total-variation of the tissue specific magnetization. This makes the algorithm more flexible to local changes in perfusion. We show that this method can be used to denoise ASL images without mixing the WM and GM signal.

Medical images often consist of low-contrast objects corrupted by random noise arising in the image acquisition process. Thus, image denoising is one of the fundamental tasks required by medical imaging analysis. In this work, we investigate an adaptive denoising scheme based on the nonlocal (NL)-means algorithm for medical imaging applications. In contrast with the traditional NL-means algorithm, the proposed adaptive NL-means (ANL-means) denoising scheme has three unique features. First, it employs the singular value decomposition (SVD) method and the K-means clustering (K-means) technique for robust classification of blocks in noisy images. Second, the local window is adaptively adjusted to match the local property of a block. Finally, a rotated block matching algorithm is adopted for better similarity matching. Experimental results from both additive white Gaussian noise (AWGN) and Rician noise are given to demonstrate the superior performance of the proposed ANL denoising technique over various image denoising benchmarks in term of both PSNR and perceptual quality comparison.

Patients suffering from cerebral ischemia or subarachnoid hemorrhage, undergo a 4D (3D+time) CT Perfusion (CTP) scan to assess the cerebral perfusion and a CT Angiography (CTA) scan to assess the vasculature. The aim of our research is to extract the vascular information from the CTP scan. This requires thin-slice CTP scans that suffer from a substantial amount of noise. Therefore noise reduction is an important prerequisite for further analysis. So far, the few noise filtering methods for 4D datasets proposed in literature deal with the temporal dimension as a 4th dimension similar to the 3 spatial dimensions, mixing temporal and spatial intensity information. We propose a bilateral noise reduction method based on time-intensity profile similarity (TIPS), which reduces noise while preserving temporal intensity information. TIPS was compared to 4D bilateral filtering on 10 patient CTP scans and, even though TIPS bilateral filtering is much faster, it results in better vessel visibility and higher image quality ranking (observer study) than 4D bilateral filtering.

We present an affinity-based optimization method for nearly-automatic vessels segmentation in CTA scans. The desired segmentation is modeled as a function that minimizes a quadratic affinity-based functional. The functional incorporates intensity and geometrical vessel shape information and a smoothing constraint. Given a few user-defined seeds, the minimum of the functional is obtained by solving a single set of linear equations. The binary segmentation is then obtained by applying a user-selected threshold. The advantages of our method are that it requires fewer initialization seeds, is robust, and yields better results than existing graph-based interactive segmentation methods. Experimental results on 20 vessel segments including the carotid arteries bifurcation and noisy parts of the carotid yield a mean symmetric surface error of 0.54mm (std=0.28).

We introduce a new 3D curved tubular intensity model in conjunction with a model fitting scheme for accurate segmentation and quantification of thin vessels in 3D tomographic images. The curved tubular model is formulated based on principles of the image formation process, and we have derived an analytic solution for the model function. In contrast to previous straight models, the new model allows to accurately represent curved tubular structures, to directly estimate the local curvature by model fitting, as well as to more accurately estimate the shape and other parameters of tubular structures. We have successfully applied our approach to 3D synthetic images as well as 3D MRA and 3D CTA vascular images of the human. It turned out that we achieved more accurate segmentation results in comparison to using a straight model.

The analysis of echocardiograms, whether visual or automated, is often hampered by ultrasound artifacts which obscure the moving myocardial wall. In this study, a probabilistic framework for tracking the endocardial surface in 3D ultrasound images is proposed, which distinguishes between visible and artifact-obscured myocardium. Motion estimation of visible myocardium relies more using a local, data-driven tracker, whereas tracking of obscured myocardium is assisted by a global, statistical model of cardiac motion. To make this distinction, the expectation-maximization algorithm is applied in a stationary and dynamic frame-of-reference. Evaluation on 35 three-dimensional echocardiographic sequences shows that this artifact-aware tracker gives better results than when no distinction is made. In conclusion, the proposed tracker is able to reduce the influence of artifacts, potentially improving quantitative analysis of clinical quality echocardiograms.

The performance of the k-nearest neighborhoods (k-NN) classifier is highly dependent on the distance metric used to identify the k nearest neighbors of the query points. The standard Euclidean distance is commonly used in practice. This paper investigates the performance of k-NN classifier with respect to different adaptive metrics in the context of medical imaging. We propose using adaptive metrics such that the structure of the data is better described, introducing some unsupervised learning knowledge in k-NN. We investigated four different metrics are estimated: a theoretical metric based on the assumption that images are drawn from Brownian Image Model (BIM), the normalized metric based on variance of the data, the empirical metric is based on the empirical covariance matrix of the unlabeled data, and an optimized metric obtained by minimizing the classification error. The spectral structure of the empirical covariance also leads to Principal Component Analysis (PCA) performed on it which results the subspace metrics. The metrics are evaluated on two data sets: lateral X-rays of the lumbar aortic/spine region, where we use k-NN for performing abdominal aorta calcification detection; and mammograms, where we use k-NN for breast cancer risk assessment. The results show that appropriate choice of metric can improve classification.

In oncological therapy monitoring, the estimation of tumor growth from consecutive CT scans is an important aspect in deciding whether the given treatment is adequate for the patient. This can be done by measuring and comparing the volume of a lesion in the scans based on a segmentation. However, simply counting the voxels within the segmentation mask can lead to significant differences in the volume, if the lesion has been segmented slightly differently by various readers or in different scans, due to the limited spatial resolution of CT and due to partial volume effects. We present a novel algorithm for measuring the volume of liver metastases and lymph nodes which considers partial volume effects at the surface of a lesion. Our algorithm is based on a spatial subdivision of the segmentation. We have evaluated the algorithm on a phantom and a multi-reader study. Our evaluations have shown that our algorithm allows determining the volume more accurately even for larger slice thicknesses. Moreover, it reduces inter-observer variability of volume measurements significantly. The calculation of the volume takes 2 seconds for 50³ voxels on a single 2.66GHz Intel Core2 CPU.

Diabetic Retinopathy (DR) is one of the leading causes of blindness in the industrialized world. Early detection is the key in providing effective treatment. However, the current number of trained eye care specialists is inadequate to screen the increasing number of diabetic patients. In recent years, automated and semi-automated systems to detect DR with color fundus images have been developed with encouraging, but not fully satisfactory results. In this study we present the initial results of a new technique for the detection and localization of microaneurysms, an early sign of DR. The algorithm is based on three steps: candidates selection, the actual microaneurysms detection and a final probability evaluation. We introduce the new Radon Cliff operator which is our main contribution to the field. Making use of the Radon transform, the operator is able to detect single noisy Gaussian-like circular structures regardless of their size or strength. The advantages over existing microaneurysms detectors are manifold: the size of the lesions can be unknown, it automatically distinguishes lesions from the vasculature and it provides a fair approach to microaneurysm localization even without post-processing the candidates with machine learning techniques, facilitating the training phase. The algorithm is evaluated on a publicly available dataset from the Retinopathy Online Challenge.

In this study, clinically produced multiphase CT volumetric data sets (pre-contrast, arterial and venous enhanced phase) are drawn upon to transcend the intrinsic limitations of single phase data sets for the robust and accurate segmentation of the liver in typically challenging cases. As an initial step, all other phase volumes are registered to either the arterial or venous phase volume by a symmetric nonlinear registration method using mutual information as similarity metric. Once registered, the multiphase CT volumes are pre-filtered to prepare for subsequent steps. Under the assumption that the intensity vectors of different organs follow the Gaussian Mixture model (GMM), expectation maximization (EM) is then used to classify the multiphase voxels into different clusters. The clusters for liver parenchyma, vessels and tumors are combined together and provide the initial liver mask that is used to generate initial zeros level set. Conversely, the voxels classified as non-liver will guide the speed image of the level sets in order to reduce leakage. Geodesic active contour level set using the gradient vector flow (GVF) derived from one of the enhanced phase volumes is then performed to further evolve the liver segmentation mask. Using EM clusters as the reference, the resulting liver mask is finally morphologically post-processed to add missing clusters and reduce leakage. The proposed method has been tested on the clinical data sets of ten patients with relatively complex and/or extensive liver cancer or metastases. A 95.8% dice similarity index when compared to expert manual segmentation demonstrates the high performance and the robustness of our proposed method - even for challenging cancer data sets - and confirms the potential of a more thorough computational exploitation of currently available clinical data sets.

A novel method is presented for definition of search lines in a variety of surface segmentation approaches. The method is inspired by properties of electric field direction lines and is applicable to general-purpose n-D shapebased image segmentation tasks. Its utility is demonstrated in graph construction and optimal segmentation of multiple mutually interacting objects. The properties of the electric field-based graph construction guarantee that inter-object graph connecting lines are non-intersecting and inherently covering the entire object-interaction space. When applied to inter-object cross-surface mapping, our approach generates one-to-one and all-to-all vertex correspondent pairs between the regions of mutual interaction. We demonstrate the benefits of the electric field approach in several examples ranging from relatively simple single-surface segmentation to complex multiobject multi-surface segmentation of femur-tibia cartilage. The performance of our approach is demonstrated in 60 MR images from the Osteoarthritis Initiative (OAI), in which our approach achieved a very good performance as judged by surface positioning errors (average of 0.29 and 0.59 mm for signed and unsigned cartilage positioning errors, respectively).

Wireless Capsule Endoscopy (WCE) integrates wireless transmission with image and video technology. It has been used to examine the small intestine non invasively. Medical specialists look for signicative events in the WCE video by direct visual inspection manually labelling, in tiring and up to one hour long sessions, clinical relevant frames. This limits the WCE usage. To automatically discriminate digestive organs such as esophagus, stomach, small intestine and colon is of great advantage. In this paper we propose to use textons for the automatic discrimination of abrupt changes within a video. In particular, we consider, as features, for each frame hue, saturation, value, high-frequency energy content and the responses to a bank of Gabor filters. The experiments have been conducted on ten video segments extracted from WCE videos, in which the signicative events have been previously labelled by experts. Results have shown that the proposed method may eliminate up to 70% of the frames from further investigations. The direct analysis of the doctors may hence be concentrated only on eventful frames. A graphical tool showing sudden changes in the textons frequencies for each frame is also proposed as a visual aid to find clinically relevant segments of the video.

While geodesic active contours (GAC) have become very popular tools for image segmentation, they are sensitive to model initialization. In order to get an accurate segmentation, the model typically needs to be initialized very close to the true object boundary. Apart from accuracy, automated initialization of the objects of interest is an important pre-requisite to being able to run the active contour model on very large images (such as those found in digitized histopathology). A second limitation of GAC model is that the edge detector function is based on gray scale gradients; color images typically being converted to gray scale prior to computing the gradient. For color images, however, the gray scale gradient results in broken edges and weak boundaries, since the other channels are not exploited for the gradient determination. In this paper we present a new geodesic active contour model that is driven by an accurate and rapid object initialization scheme-weighted mean shift normalized cuts (WNCut). WNCut draws its strength from the integration of two powerful segmentation strategies-mean shift clustering and normalized cuts. WNCut involves first defining a color swatch (typically a few pixels) from the object of interest. A multi-scale mean shift coupled normalized cuts algorithm then rapidly yields an initial accurate detection of all objects in the scene corresponding to the colors in the swatch. This detection result provides the initial boundary for GAC model. The edge-detector function of the GAC model employs a local structure tensor based color gradient, obtained by calculating the local min/max variations contributed from each color channel (e.g. R,G,B or H,S,V). Our color gradient based edge-detector function results in more prominent boundaries compared to classical gray scale gradient based function. We evaluate segmentation results of our new WNCut initialized color gradient based GAC (WNCut-CGAC) model against a popular region-based model (Chan & Vese) on a total of 60 digitized histopathology images. Across a total of 60 images, the WNCut-CGAC model yielded an average overlap, sensitivity, specificity, and positive predictive value of 73%, 83%, 97%, 84%, compared to the Chan & Vese model which had corresponding values of 64%, 75%, 95%, 72%. The rapid and accurate object initialization scheme (WNCut) and the color gradient make the WNCut-CGAC scheme, an ideal segmentation tool for very large, color imagery.

An atlas provides a reference anatomic structure and an associated coordinate system. An atlas may be used in a variety of applications, including segmentation and registration, and can be used to characterize anatomy across a population. We present a method for generating an atlas of the human retina from 500 color fundus image pairs. Using color fundus image pairs, we register image pairs to obtain a larger anatomic field of view. Key retinal anatomic features are selected for atlas landmarks: disk center, fovea, and main vessel arch. An atlas coordinate system is defined based on the statistics of the landmarks. Images from the population are warped into the atlas space to produce a statistical retinal atlas which can be used for automatic diagnosis, concise indexing, semantic blending, etc.

The topogram is a 2D projection image of human body formed using a Computed Tomography (CT) scanner. It could be used to delimitate the desired scan range for further precise 3D CT scan. In this paper, we present a robust and efficient system for automatically determining scan ranges and their associated anatomical landmarks for topogram images. The system could handle the cases when only about 50% of the desired regions are visible. The robustness of our system can be attributed to three key ingredients: 1. The detection is based on a hierarchical network; 2. Network optimization is based on sequentially optimizing a set of subnetworks; 3. The detection probability is further refined based on the detection context. Extensive experiments (including external testing) on over 1000 topogram images show that our approach works robustly and efficiently even on very challenging data.

Detecting pigmented network is a crucial step for melanoma diagnosis. In this paper, we present a novel graphbased pigment network detection method that can find and visualize round structures belonging to the pigment network. After finding sharp changes of the luminance image by an edge detection function, the resulting binary image is converted to a graph, and then all cyclic sub-graphs are detected. Theses cycles represent meshes that belong to the pigment network. Then, we create a new graph of the cyclic structures based on their distance. According to the density ratio of the new graph of the pigment network, the image is classified as "Absent" or "Present". Being Present means that a pigment network is detected in the skin lesion. Using this approach, we achieved an accuracy of 92.6% on five hundred unseen images.

We proposed an approach based on reconstructive sparse representations to segment tissues in optical images of the uterine cervix. Because of large variations in image appearance caused by the changing of the illumination and specular reflection, the color and texture features in optical images often overlap with each other and are not linearly separable. By leveraging sparse representations the data can be transformed to higher dimensions with sparse constraints and become more separated. K-SVD algorithm is employed to find sparse representations and corresponding dictionaries. The data can be reconstructed from its sparse representations and positive and/or negative dictionaries. Classification can be achieved based on comparing the reconstructive errors. In the experiments we applied our method to automatically segment the biomarker AcetoWhite (AW) regions in an archive of 60,000 images of the uterine cervix. Compared with other general methods, our approach showed lower space and time complexity and higher sensitivity.

Disorders of the heart valves constitute a considerable health problem and often require surgical intervention. Recently various approaches were published seeking to overcome the shortcomings of current clinical practice,that still relies on manually performed measurements for performance assessment. Clinical decisions are still based on generic information from clinical guidelines and publications and personal experience of clinicians. We present a framework for retrieval and decision support using learning based discriminative distance functions and visualization of patient similarity with relative neighborhood graphsbased on shape and derived features. We considered two learning based techniques, namely learning from equivalence constraints and the intrinsic Random Forest distance. The generic approach enables for learning arbitrary user-defined concepts of similarity depending on the application. This is demonstrated with the proposed applications, including automated diagnosis and interventional suitability classification, where classification rates of up to 88.9% and 85.9% could be observed on a set of valve models from 288 and 102 patients respectively.

In this paper, we present a shape based segmentation methodology for magnetic resonance prostate images. We first propose a new way to represent shapes via the hyperbolic tangent of the signed distance function. This effectively corrects the drawbacks of the signed distance function and yields very reasonable results for the shape registration and learning. Secondly, under a Bayesian statistical framework, instead of computing the posterior using a uniform prior, a directional distance map is introduced in order to incorporate a priori knowledge of image content as well as the estimated center of target object. Essentially, the image is modeled as a Finsler manifold and the metric is computed out of the directional derivative of the image. Then the directional distance map is computed to suppress the posterior remote from the object center. Thirdly, in the posterior image, a localized region based cost functional is designed to drive the shape based segmentation. Such cost functional utilizes the local regional information and is robust to both image noise and remote/irrelevant disturbances. With these three major components, the entire shape based segmentation procedure is provided as a complete open source pipeline and is applied to magnetic resonance image (MRI) prostate data.

We describe a simple and robust algorithm for estimating 3D shape given a number of silhouette points obtained from two or more viewpoints and a parametric model of the shape. Our algorithm minimizes (in the least squares sense) the distances from the lines obtained by unprojecting the silhouette points to 3D to their closest silhouette points on the 3D shape. The solution is found using an iterative approach. In each iteration, we locally approximate the least squares problem with a degree-4 polynomial function. The approximate problem is solved using a nonlinear conjugate gradient solver that takes advantage of its structure to perform exact and global line searches. We tested our algorithm by applying it to reconstruct patient-specific femur shapes from simulated biplanar X-ray images.

Shape description is an important step in image analysis. The skeleton is used as a simple, compact representation of a shape. A skeleton represents the line centered in the shape and must be homotopic and one point wide. Current skeletonization algorithms compute the skeleton over several image scans, using either thinning algorithms or distance transforms. The principle of thinning is to delete points as one goes along, preserving the topology of the shape. On the other hand, the maxima of the local distance transform identifies the skeleton and is an equivalent way to calculate the medial axis. However, with this method, the skeleton obtained is disconnected so it is required to connect all the points of the medial axis to produce the skeleton. In this study we introduce a translated distance transform and adapt an existing distance driven homotopic algorithm to perform skeletonization with a single scan and thus allow the processing of unbounded images. This method is applied, in our study, on micro scanner images of trabecular bones. We wish to characterize the bone micro architecture in order to quantify bone integrity.

In this article, we propose a unified statistical framework for image segmentation with shape prior information. The approach combines an explicitely parameterized point-based probabilistic statistical shape model (SSM) with a segmentation contour which is implicitly represented by the zero level set of a higher dimensional surface. These two aspects are unified in a Maximum a Posteriori (MAP) estimation where the level set is evolved to converge towards the boundary of the organ to be segmented based on the image information while taking into account the prior given by the SSM information. The optimization of the energy functional obtained by the MAP formulation leads to an alternate update of the level set and an update of the fitting of the SSM. We then adapt the probabilistic SSM for multi-shape modeling and extend the approach to multiple-structure segmentation by introducing a level set function for each structure. During segmentation, the evolution of the different level set functions is coupled by the multi-shape SSM. First experimental evaluations indicate that our method is well suited for the segmentation of topologically complex, non spheric and multiple-structure shapes. We demonstrate the effectiveness of the method by experiments on kidney segmentation as well as on hip joint segmentation in CT images.

We present a 3D extension and validation of an intra-operative registration framework that accommodates tissue resection. The framework is based on the bijective Demons method, but instead of regularizing with the traditional Gaussian smoother, we apply an anisotropic diffusion filter with the resection modeled as a diffusion sink. The diffusion sink prevents unwanted Demon forces that originates from the resected area from diffusing into the surrounding area. Another attractive property of the diffusion sink is the resulting continuous deformation field across the diffusion sink boundary, which allows us to move the boundary of the diffusion sink without changing values in the deformation field. The area of resection is estimated by a level-set method evolving in the space of image intensity disagreements in the intra-operative image domain. A product of using the bijective Demons method is that we can also provide an accurate estimate of the resected tissue in the preoperative image space. Validation of the proposed method was performed on a set of 25 synthetic images. Our experiments show a significant improvement in accommodating resection using the proposed method compared to two other Demons based methods.

The iterative closest point (ICP) algorithm is probably the most popular algorithm for fine registration of surfaces. Among its key properties are: a simple minimization scheme, proofs of convergence as well as the easiness to modify and improve it in many ways (e.g. use of fuzzy point correspondences, incorporation of a priori knowledge, extensions to non-linear deformations, speed-up strategies, etc.) while keeping the desirable properties of the original method. However, most ICP-like registration methods suffer from the fact that they only consider the distance between the surfaces to register in the criterion to minimize, and thus are highly dependent on how the surfaces are aligned in the first place. This explains why these methods are likely to be trapped in local minima and to lead to erroneous solutions. A solution to partly alleviate this problem would consist in adding higher-order information in the criterion to minimize (e.g. normals, curvatures, etc.), but previous works along these research tracks have led to computationally intractable minimization schemes. In this paper, we propose a new way to include the point unit normals in addition to the point coordinates to derive an ICP-like scheme for non-linear registration of surfaces, and we show how to keep the properties of the original ICP algorithm. Our algorithm rests on a simple formula showing how the unit normal changes when a surface undergoes a small deformation. The use of this formula in an ICP-like algorithm is made possible by adequate implementation choices, most notably the use of a local, differentiable, parametrization of the surfaces and a locally affine deformation model using this local parametrization. Then we experimentally show the strong added value of using the unit normals in a series of controlled experiments.

A major focus of computational anatomy is to extract the most relevant information to identify and characterize anatomical variability within a group of subjects as well as between different groups. The construction of atlases is central to this effort. An atlas is a deterministic or probabilistic model with intensity variance, structural, functional or biochemical information over a population. To date most algorithms to construct atlases have been based on a single subject assuming that the population is best described by a single atlas. However, we believe that in a population with a wide range of subjects multiple atlases may be more representative since they reveal the anatomical differences and similarities within the group. In this work, we propose to use the K-means clustering algorithm to partition a set of images into several subclasses, based on a joint distance which is composed of a distance quantifying the deformation between images and a dissimilarity measured from the registration residual. During clustering, the spatial transformations are averaged rather than images to form cluster centers, to ensure a crisp reference. At the end of this algorithm, the updated centers of the k clusters are our atlases. We demonstrate this algorithm on a subset of a public available database with whole brain volumes of subjects aged 18-96 years. The atlases constructed by this method capture the significant structural differences across the group.

Robust point matching (RPM) jointly estimates correspondences and non-rigid warps between unstructured point-clouds. RPM does not, however, utilize information of the topological structure or group memberships of the data it is matching. In numerous medical imaging applications, each extracted point can be assigned group membership attributes or labels based on segmentation, partitioning, or clustering operations. For example, points on the cortical surface of the brain can be grouped according to the four lobes. Estimated warps should enforce the topological structure of such point-sets, e.g. points belonging to the temporal lobe in the two point-sets should be mapped onto each other. We extend the RPM objective function to incorporate group membership labels by including a Label Entropy (LE) term. LE discourages mappings that transform points within a single group in one point-set onto points from multiple distinct groups in the other point-set. The resulting Labeled Point Matching (LPM) algorithm requires a very simple modification to the standard RPM update rules. We demonstrate the performance of LPM on coronary trees extracted from cardiac CT images. We partitioned the point sets into coronary sections without a priori anatomical context, yielding potentially disparate labelings (e.g. [1,2,3] → [a,b,c,d]). LPM simultaneously estimated label correspondences, point correspondences, and a non-linear warp. Non-matching branches were treated wholly through the standard RPM outlier process akin to non-matching points. Results show LPM produces warps that are more physically meaningful than RPM alone. In particular, LPM mitigates unrealistic branch crossings and results in more robust non-rigid warp estimates.

3D registration of brain MRI data is vital for many medical imaging applications. However, purely intensitybased approaches for inter-subject matching of brain structure are generally inaccurate in cortical regions, due to the highly complex network of sulci and gyri, which vary widely across subjects. Here we combine a surfacebased cortical registration with a 3D fluid one for the first time, enabling precise matching of cortical folds, but allowing large deformations in the enclosed brain volume, which guarantee diffeomorphisms. This greatly improves the matching of anatomy in cortical areas. The cortices are segmented and registered with the software Freesurfer. The deformation field is initially extended to the full 3D brain volume using a 3D harmonic mapping that preserves the matching between cortical surfaces. Finally, these deformation fields are used to initialize a 3D Riemannian fluid registration algorithm, that improves the alignment of subcortical brain regions. We validate this method on an MRI dataset from 92 healthy adult twins. Results are compared to those based on volumetric registration without surface constraints; the resulting mean templates resolve consistent anatomical features both subcortically and at the cortex, suggesting that the approach is well-suited for cross-subject integration of functional and anatomic data.

The Hamilton-Jacobi equation (HJE) appears widely in applied mathematics, physics, and optimal control theory. While its analytical solution is rarely available, the numerical solver is indispensable. In this work, firstly we propose a novel numerical method, based on the fast sweeping scheme, for the static HJE. Comparing with the original fast sweeping method, our algorithm speeds up the solution up to 8 times in 3D. The efficiency is due to incorporating the ideas of the fast marching into the fast sweeping. Essentially, the sweeping origin is selected so that the sweeping direction is more consistent with the information flow direction and the regions where the two directions are against are avoided. Moreover, the successive-overrelaxation nonlinear iterative method is used for faster convergence. Secondly, we provide a complete pipeline for brain tractography, in which the proposed solver is the key component for finding the optimal fiber tracts. Besides, the pipeline contains components from orientation distribution function estimation, multiple fiber extraction to the final fiber bundle volumetric segmentation, completing the process from DW-MRI image to segmented fiber bundles. The pipeline is integrated into the publicly available software 3D Slicer. The new solver has been tested and compared with the original scheme on various types of HJEs and the tractography pipeline was tested and performed consistently on all the 12 brain DW-MRI images.

Q-ball imaging (QBI) is a diffusion imaging technique that provides unique microstructural information and outperforms diffusion tensor imaging (DTI) in areas of fiber crossings. The generalized fractional anisotropy (GFA) is a widely accepted quantitative measure of fiber integrity for QBI. In this paper we demonstrate a major drawback of the GFA occuring in crossing fiber regions: the fiber integrity is heavily underestimated here. We present a new, directional anisotropy measure and compare it to the GFA using custom-designed fiber phantoms and Monte Carlo simulations. Furthermore, we present in-vivo data of the measure. With the expense of reduced CNR, the new measure allows a better quantification of fiber bundles in crossing regions and potentially yields additional information regarding non-dominant tracts within a voxel.

Diffusion tensor imaging (DTI) is widely used to characterize tissue micro-architecture and brain connectivity. Yet DTI suffers serious limitations in regions of crossing fibers because traditional tensor techniques cannot represent multiple, independent intra-voxel orientations. Compressed sensing has been proposed to resolve crossing fibers using a tensor mixture model (e.g., Crossing Fiber Angular Resolution of Intra-voxel structure, CFARI). Although similar in spirit to deconvolution approaches, CFARI uses sparsity to stabilize estimation with limited data rather than spatial consistency or limited model order. Here, we extend the CFARI approach to resolve crossing fibers through a strictly positive, parsimonious mixture model. Together with an optimized preconditioned conjugate gradient solver, estimation error and computational burden are greatly reduced over the initial presentation. Reliable estimates of intra-voxel orientations are demonstrated in simulation and in vivo using data representative of typical, low b-value (30 directions, 700 s/mm²) clinical DTI protocols. These sequences are achievable in 5 minutes at 3 T, and the whole brain CFARI analysis is tractable for routine analysis. With these improvements, CFARI provides a robust framework for identifying intra-voxel structure with conventional DTI protocols and shows great promise in helping to resolve the crossing fiber problem in current clinical imaging studies.

Recent studies reported the development of methods for rigid registration of 2D fetal brain imaging data to correct for unconstrained fetal and maternal motion, and allow the formation of a true 3D image of conventional fetal brain anatomy from conventional MRI. Diffusion tensor imaging provides additional valuable insight into the developing brain anatomy, however the correction of motion artifacts in clinical fetal diffusion imaging is still a challenging problem. This is due to the challenging problem of matching lower signal-to-noise ratio diffusion weighted EPI slice data to recover between-slice motion, compounded by the presence of possible geometric distortions in the EPI data. In addition, the problem of estimating a diffusion model (such as a tensor) on a regular grid that takes into account the inconsistent spatial and orientation sampling of the diffusion measurements needs to be solved in a robust way. Previous methods have used slice to volume registration within the diffusion dataset. In this work, we describe an alternative approach that makes use of an alignment of diffusion weighted EPI slices to a conventional structural MRI scan which provides a geometrically correct reference image. After spatial realignment of each diffusion slice, a tensor field representing the diffusion profile is estimated by weighted least squared fitting. By qualitative and quantitative evaluation of the results, we confirm the proposed algorithm successfully corrects the motion and reconstructs the diffusion tensor field.

The functional networks, extracted from fMRI images using independent component analysis, have been demonstrated informative for distinguishing brain states of cognitive functions and neurological diseases. In this paper, we propose a novel algorithm for discriminant analysis of functional networks encoded by spatial independent components. The functional networks of each individual are used as bases for a linear subspace, referred to as a functional connectivity pattern, which facilitates a comprehensive characterization of temporal signals of fMRI data. The functional connectivity patterns of different individuals are analyzed on the Grassmann manifold by adopting a principal angle based subspace distance. In conjunction with a support vector machine classifier, a forward component selection technique is proposed to select independent components for constructing the most discriminative functional connectivity pattern. The discriminant analysis method has been applied to an fMRI based schizophrenia study with 31 schizophrenia patients and 31 healthy individuals. The experimental results demonstrate that the proposed method not only achieves a promising classification performance for distinguishing schizophrenia patients from healthy controls, but also identifies discriminative functional networks that are informative for schizophrenia diagnosis.

Automatic image analysis of magnetic resonance (MR) images of the knee is simplified by bringing the knee into a reference position. While the knee is typically put into a reference position during image acquisition, this alignment will generally not be perfect. To correct for imperfections, we propose a two-step process of bone segmentation followed by elastic tissue deformation. The approach makes use of a fully-automatic segmentation of femur and tibia from T1 and T2* images. The segmentation algorithm is based on a continuous convex optimization problem, incorporating regional, and shape information. The regional terms are included from a probabilistic viewpoint, which readily allows the inclusion of shape information. Segmentation of the outer boundary of the cortical bone is encouraged by adding simple appearance-based information to the optimization problem. The resulting segmentation without the shape alignment step is globally optimal. Standard registration is problematic for knee alignment due to the distinct physical properties of the tissues constituting the knee (bone, muscle, etc.). We therefore develop an alternative alignment approach based on a simple elastic deformation model combined with strict enforcement of similarity transforms for femur and tibia based on the obtained segmentations.

Segmentation and representation of human brain cortex from Magnetic Resonance (MR) images is an important step for visualization and analysis in many neuro imaging applications. In this paper, we propose an automatic and fast algorithm to segment the brain cortex and to represent it as a geometric surface on which analysis can be carried out. The algorithm works on T1 weighted MR brain images with extracranial tissue removed. A fuzzy clustering algorithm with a parametric bias field model is applied to assign membership values of gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) to each voxel. The cortical boundaries, namely the WM-GM and GM-CSF boundary surfaces, are extracted as iso-surfaces of functions derived from these membership functions. The central surface (CS), which traces the peak values (or ridges) of the GM membership function, is then extracted using gradient vector diffusion. Our main contribution is to provide a generic, accurate, fast, yet fully-automatic approach to (i) produce a soft segmentation of the MR brain image with intensity field correction, (ii) extract both the boundary and the center of the cortex in a surface form, where the topology and geometry can be explicitly examined, and (iii) use the extracted surfaces to model the curvy, folding cortical volume, which allows an intuitive measurement of the thickness. As a demonstration, we compute cortical thickness from the surfaces and compare the results with what has been reported in the literature. The entire process from raw MR image to cortical surface reconstruction takes on average between five to ten minutes.