Image reconstruction from nonuniformly sampled frequency domain data is an important problem that arises in computed imaging. The current reconstruction techniques suffer from fundamental limitations in their model and implementation that result in blurred reconstruction and/or artifacts. Here, we present a new approach for solving this problem that relies on a more realistic model and involves an explicit measure for the reconstruction accuracy that is optimized iteratively. The image is assumed piecewise constant to impose practical display constraints using pixels. We express the mapping of these unknown pixel values to the available frequency domain values as a linear system. Even though the system matrix is shown to be dense and too large to solve for practical purposes, we observe that applying a simple orthogonal transformation to the rows of this matrix converts the matrix into a sparse format. The transformed system is subsequently solved using the conjugate gradient method. The proposed method is applied to reconstruct images of a numerical phantom as well as actual magnetic resonance images using spiral sampling. The results support the theory and show that the computational load of this method is similar to that of other techniques. This suggests its potential for practical use.

We investigate a new, provably convergent OSEM-like (ordered-subsets expectation-maximization) reconstruction algorithm for emission tomography. The new algorithm, which we term C-OSEM (complete-data OSEM), can be shown to monotonically increase the log-likelihood at each iteration. The familiar ML-EM reconstruction algorithm for emission tomography can be derived in a novel way. One may write a single objective function with complete, incomplete data and the reconstruction variables as in the EM approach. But in the objective function approach, there is no E-step. Instead, a suitable alternating descent on the complete data and then the reconstruction variables results in two update equations that can be shown to be equivalent to the familiar EM algorithm. Hence, minimizing this objective becomes equivalent to maximizing the likelihood. We derive our C-OSEM algorithm by modifying the above approach to update the complete data only along ordered subsets. The resulting update equation is quite different from OSEM, but still retains the speed-enhancing feature of the updates due to the limited backprojection facilitated by the ordered subsets. Despite this modification, we are able to show that the objective function decreases at each iteration, and (given a few more mild assumptions regarding the number of fixed points) conclude that the C-OSEM algorithm provides a monotonic convergence toward the maximum likelihood solution. We simulated noisy and noiseless emission projection data, and reconstructed them using the ML-EM, and the proposed C-OSEM with 4 subsets. We also reconstruct the data using the OSEM method. Anecdotal results show that the C-OSEM algorithm is much faster than ML-EM though slower than OSEM.

Based on Kunyansky's and our previous work, an efficient, analytical solution to the reconstruction problem of myocardial perfusion SPECT has been developed that allows simultaneous compensation for non-uniform attenuation, scatter, and system-dependent resolution variation, as well as suppression of signal-dependent Poisson noise. To avoid reconstructed images being corrupted by the presence of Poisson noise, a Karhunen-Loeve (K-L) domain adaptive Wiener filter is applied first to suppress the noise in the primary- and scatter-window measurements. The scatter contribution to the primary-energy-window measurements is then removed by our scatter estimation method, which is energy spectrum based, modified from the triple-energy-window acquisition protocol. The resolution variation is corrected by the depth-dependent deconvolution, which, being based on the central-ray approximation and the distance-frequency relation, deconvolves the scatter-free data with the accurate detector-response kernel in frequency domain. Finally, the deblurred projection data are analytically reconstructed with non-uniform attenuation by an algorithm based on Novikov's explicit inversion formula. The preliminary Monte Carlo simulation results using a realistic human thoracic phantom demonstrate that, for parallel-beam geometry, the proposed analytical reconstruction scheme is computationally comparable to filtered backprojection and quantitatively equivalent to iterative maximum a posteriori expectation-maximization reconstruction. Extension to other geometries is under progress.

We propose an alternating minimization (AM) image estimation algorithm for iteratively reconstructing transmission tomography images. The algorithm is based on a model that accounts for much of the underlying physics, including Poisson noise in the measured data, beam hardening of polyenergetic radiation, energy dependence of the attenuation coefficients and scatter. It is well-known that these nonlinear phenomena can cause severe artifacts throughout the image when high-density objects are present in soft tissue, especially when using the conventional technique of filtered back projection (FBP). If we assume no prior knowledge of the high-density object(s), our proposed algorithm yields much improved images in comparison to FBP, but retains significant streaking between the high-density regions. When we incorporate the knowledge of the attenuation and pose parameters of the high-density objects into the algorithm, our simulations yield images with greatly reduced artifacts. To accomplish this, we adapted the algorithm to perform a search at each iteration (or after every n iterations) to find the optimal pose of the object before updating the image. The final iteration returns pose values within 0.1 millimeters and 0.01 degrees of the actual location of the high-density structures.

Dual-energy (DE) X-ray computed tomography (CT) has shown promise for material characterization and for providing quantitatively accurate CT values in a variety of applications. However, DE-CT has not been used routinely in medicine to date, primarily due to dose considerations. Most methods for DE-CT have used the filtered backprojection method for image reconstruction, leading to suboptimal noise/dose properties. This paper describes a statistical (maximum-likelihood) method for dual-energy X-ray CT that accommodates a wide variety of potential system configurations and measurement noise models. Regularized methods (such as penalized-likelihood or Bayesian estimation) are straightforward extensions. One version of the algorithm monotonically decreases the negative log-likelihood cost function each iteration. An ordered-subsets variation of the algorithm provides a fast and practical version.

3D image reconstruction from cone-beam projections on a circular short-scan is analyzed using the tools of CB reconstruction theory. By definition, CB Data on a circular short-scan do not provide enough information for exact reconstruction outside the plane of the source path. However, circular short-scans are attractive in many applications due to practical constraints on the data acquisition geometry. This work attemps at determining the optimal tomographic capabilities of CB imaging with circular short-scan. Two new algorithms are suggested for reconstruction and tested on the FORBILD head phantom with a large cone angle. The results are very encouraging.

Since the recent introduction of multi-slice helical computed tomography (MHCT), new clinical applications have experienced tremendous growth in recent years. MHCT offers improved volume coverage, faster scan speed, more isotropic spatial resolution, and reduced x-ray tube loading. Similar to the single slice helical CT, the projection data collected in MHCT is inherently inconsistent due to the constant table motion. In addition, cone beam effects in MHCT produce additional complexity and image artifacts. Although the cone angle is quite smaller even for the 16-slice configuration, the impact on image artifacts cannot be ignored. Many reconstruction algorithms have been proposed and investigated recently to combat image artifacts associated with the MHCT data acquisition. In this paper, we propose a cone-angle dependent generalized weighting scheme for 16-slice helical CT that allows the production of MHCT images with only 2D backprojection. The cone-angle dependency of the algorithm suppresses image artifacts due to the cone beam effect and the generalized weighting portion enables interpolation be performed with conjugate samples for the 16-slice helical dataset. With the proposed algorithm, image artifacts are significantly reduced.

This paper describes an automated 3-D surface recovery algorithm for consistent and quantitative evaluation of the deformation in the ONH (optic nerve head). Additional measures, such as the changes in the volume of the cup and the disc as an improvement to the traditional cup to disc ratios, can thus be developed for longitudinal follow-up study of a patient. We propose an automated computerized technique for stereo pair registration and surface visualization of the ONH. Power cepstrum and zero mean cross correlation are embedded in the registration and a 3-D surface recovery technique is proposed. Preprocessing, as well as an overall registration, is performed upon stereo pairs. Then a coarse to fine feature matching strategy is used to reduce the ambiguity in finding the conjugate pair of the same point within the constraints of the epipolar plane. A cubic B-spline interpolation smooths the representation of the ONH obtained, while superimposition of features such as blood vessels is added. Studies show high correlation between traditional cup/disc measures derived from manual segmentation by ophthalmologists and computer generated cup/disc volume ratio. Such longitudinal studies over a large population of glaucoma patients are currently in progress for validation of the surface recovery algorithm.

The task of image segmentation implies estimation of the number and associated parameters of the classes within an image, and the class label for each image voxel. In this work, an over-segmentation of the data is first obtained using a Bayesian restoration algorithm. The method incorporates a novel spatial interaction prior, in which neighboring voxels can be classified differently so long as the distance between the centroids of their intensity distributions are within a certain extent. The corresponding functional is iteratively minimized using a series of local optimizations for the label field and a half-quadratic algorithm for the restoration. Redundant classes are then grouped in a second step by making use of information obtained in the initial restoration about the degree of affinity or interaction between the classes. The method is demonstrated on MRI images of the head.

During the last couple of years virtual endoscopic systems (VES) have emerged as standard tools that are nowadays close to be utilized in daily clinical practice. Such tools render hollow human structures, allowing a clinician to visualize their inside in an endoscopic-like paradigm. It is common practice that the camera of a virtual endoscope is attached to the centerline of the structure of interest, to facilitate navigation. This centerline has to be determined manually or automatically, prior to an investigation. While there exist techniques that can straightforwardly handle simple tube-like structures (e.g. colon, aorta), structures like the tracheobronchial tree still represent a challenge due to their complex branching. In these cases it is necessary to determine all branching points within the tree which is - because of the complexity - impractical to be accomplished in a manual manner. This paper presents a simultaneous segmentation/skeletonization algorithm that extracts all major airway branches and large parts of the minor distal branches (up to 7th order) using a front propagation approach. During the segmentation the algorithm keeps track of the centerline of the segmented structure and detects all branching points. This in turn allows the full reconstruction of the tracheobronchial tree.

The abdominal aorta is the most common site for an aneurysm, which may lead to hemorrhage and death, to develop. The aim of this study was to develop a semi-automated method to de-lineate the vessels and detect the center-line of these vessels to make measurements necessary for stent design from multi-detector computed tomograms. We developed a robust method of tracking the aortic vessel tree with branches from a user selected seed point along the vessel path using scale space approaches, central transformation measures, vessel direction findings, iterative corrections and a priori information in determining the vessel branches. Fifteen patients were scanned with contrast on Mx8000 CT scanner (Philips Medical Systems), with a 3.2 mm thickness, 1.5 mm slice spacing, and a stack of 512x512x320 volume data sets were reconstructed. The algorithm required an initial user input to locate the vessel seen in axial CT slice. Next, the automated image processing took approximately two minutes to compute the centerline and borders of the aortic vessel tree. The results between the manually and automatically generated vessel diameters were compared and statistics were computed. We observed our algorithm was consistent (less than 0.01 S.D) and similar (less than 0.1 S.D) to manual results.

This paper presents an extension of the theory and algorithms for fuzzy connectedness. In this framework, a strength of connectedness is assigned to every pair of image elements by finding the strongest connecting path between them. The strength of a path is the weakest affinity between successive pairs of elements along the path. Affinity specifies the degree to which elements hang together locally in the image. A fuzzy connected object containing a particular seed element is computed via dynamic programming. In all reported works so far, the minimum of affinities has been considered for path strength and the maximum of path strengths for fuzzy connectedness. The question thus remained all along as to whether there are other valid formulations for fuzzy connectedness. One of the main contributions of this paper is a theoretical investigation under reasonable axioms to establish that maximum of path strengths of minimum of affinities along each path is indeed the one and only valid choice. The second contribution here is to generalize the theory and algorithms of fuzzy connectedness to the multi-seeded case. The importance of multi-seeded fuzzy connectedness is illustrated with examples taken from several real medical imaging applications.

This paper describes the theory and algorithms of fuzzy distance transform (FDT). Fuzzy distance is defined as the length of the shortest path between two points. The length of a path in a fuzzy subset is defined as the integration of fuzzy membership values of the points along the path. The shortest path between two points is the one with the minimum length among all (infinitely many) paths between the two points. It is demonstrated that, unlike in the binary case, the shortest path in a fuzzy subset is not necessarily a straight-line segment. The support of a fuzzy subset is the set of points with nonzero membership values. It is shown that, for any fuzzy subset, fuzzy distance is a metric for the interior of its support. FDT is defined as the process on a fuzzy subset that assigns at each point the smallest fuzzy distance from the boundary of the support. The theoretical framework of FDT in continuous space is extended to digital spaces and a dynamic programming-based algorithm is presented for its computation. Several potential medical imaging applications are presented including the quantification of blood vessels and trabecular bone thickness in the regime of limited special resolution.

Segmentation of fluoroscopy images is useful for fluoroscopy-to-CT image registration. However, it is impossible to assign a unique tissue type to each pixel. Rather each pixel corresponds to an entire path of tissue types encountered along a ray from the X-ray source to the detector plate. Furthermore, there is an inherent many-to-one mapping between paths and pixel values. We address these issues by assigning to each pixel not a scalar value but a fuzzy vector of tissue probabilities. We perform this segmentation in a probabilistic way by first learning typical distributions of bone, air, and soft tissue that correspond to certain fluoroscopy image values and then assigning each value to a probability distribution over its most likely generating paths. We then evaluate this segmentation on ground truth patient data.

Accurate computation of the thickness of articular cartilage in 3D is crucial in diagnosis of joint diseases. The purpose of this research project is to develop an unsupervised method to produce three-dimensional (3D) thickness map of articular cartilage with magnetic resonance imaging (MRI). The method consists of two main parts, cartilage extraction and thickness map computation. The initial segmentation for cartilage extraction is achieved using a recently proposed algorithm which depends on region-growing. The regions produced during this process are labeled as cartilage or non-cartilage using a voting procedure which essentially depends on local 2-class clustering and makes use of prior knowledge about cartilage regions. Following cartilage extraction, femoral and tibial cartilages are separated by detecting the interface between them using a deformable model. After the separation, the cartilage surfaces are reconstructed as a triangular mesh and divided into two plates according to the relation between surface normal at each vertex and principal axes of the structure. For surface reconstruction, we propose an algorithm which incorporates a simple MR imaging model which allows surface representations with sub-voxel accuracy. Our thickness computation algorithm treats each plate separately as a deformable model while considering the other plate as the target surface towards which it is deformed. At the end of deformation, the thickness values at each vertex is defined as the distance between the locations at pre and post-deformation instances. The performance of the cartilage segmentation is compared to manual tracing. Also, the performance evaluation of the thickness computation algorithm on phantoms resulted in RMS errors on the order of 1%.

Optical coherence tomography (OCT) provides a non-contact and non-invasive means to visualize the corneal anatomy at micron scale resolution. We obtained corneal images from an arc-scanning (converging) OCT system operating at a wavelength of 830nm and a fan-shaped-scanning high-speed OCT system with an operating wavelength of 1310nm. Different scan protocols (arc/fan) and data acquisition rates, as well as wavelength dependent bio-tissue backscatter contrast and optical absorption, make the images acquired using the two systems different. We developed image-processing algorithms to automatically detect the air-tear interface, epithelium-Bowman's layer interface, laser in-situ keratomileusis (LASIK) flap interface, and the cornea-aqueous interface in both kinds of images. The overall segmentation scheme for 830nm and 1310nm OCT images was similar, although different strategies were adopted for specific processing approaches. Ultrasound pachymetry measurements of the corneal thickness and Placido-ring based corneal topography measurements of the corneal curvature were made on the same day as the OCT examination. Anterior/posterior corneal surface curvature measurement with OCT was also investigated. Results showed that automated segmentation of OCT images could evaluate anatomic outcome of LASIK surgery.

With the improvements in techniques for generating surface models from magnetic resonance (MR) images, it has recently become feasible to study the morphological characteristics of the human brain cortex in vivo. Studies of the entire surface are important for measuring global features, but analysis of specific cortical regions of interest provides a more detailed understanding of structure. We have previously developed a method for automatically segmenting regions of interest from the cortical surface using a watershed transform. Each segmented region corresponds to a cortical sulcus and is thus termed a sulcal region. In this work, we describe three important augmentations of this methodology. First, we describe a user interface that allows for the efficient labeling of the segmented sulcal regions called the Interactive Program for Sulcal Labeling (IPSL). Two additional augmentations of the methodology allow for even finer division of regions on the cortex. Both employ the fast marching technique to track curves of interest on the cortical surface. These curves are then used to separate segmented regions. Validation experiments indicate that the proposed methodology gives highly repeatable results.

In previous work, the authors presented a multi-stage procedure for the semi-automatic reconstruction of the cerebral cortex from magnetic resonance images. This method suffered from several disadvantages. First, the tissue classification algorithm used can be sensitive to noise within the image. Second, manual interaction was required for masking out undesired regions of the brain image, such as the ventricles and putamen. Third, iterated median filters were used to perform a topology correction on the initial cortical surface, resulting in an overly smoothed initial surface. Finally, the deformable surface used to converge to the cortex had difficulty capturing narrow gyri. In this work, all four disadvantages of the procedure have been addressed. A more robust tissue classification algorithm is employed and the manual masking step is replaced by an automatic method involving level set deformable models. Instead of iterated median filters, an algorithm developed specifically for topology correction is used. The last disadvantage is addressed using an algorithm that artificially separates adjacent sulcal banks. The new procedure is more automated but also more accurate than the previous one. Its utility is demonstrated by performing a preliminary study on data from the Baltimore Longitudinal Study of Aging.

Other researchers have proposed that the brain parenchymal fraction (or brain atrophy) may be a good surrogate measure for disease progression in patients with Multiple Sclerosis. This paper considers various factors influencing the measure of the brain parenchymal fraction obtained from dual spin-echo PD and T2-weighted head MRI scans. We investigate the robustness of the brain parenchymal fraction with respect to two factors: brain-mask border placement which determines the brain intra-dural volume, and brain scan incompleteness. We show that an automatic method for brain segmentation produces an atrophy measure which is fairly sensitive to the brain-mask placement. We also show that a robust, reproducible brain atrophy measure can be obtained from incomplete brain scans, using data in a centrally placed subvolume of the brain.

This paper proposes an image analysis system to segment multiple sclerosis lesions of magnetic resonance (MR) brain volumes consisting of 3 mm thick slices using three channels (images showing T1-, T2- and PD -weighted contrast). The method uses the statistical model of Markov Random Fields (MRF) both at low and high levels. The neighborhood system used in this MRF is defined in three types: (1) Voxel to voxel: a low-level heterogeneous neighborhood system is used to restore noisy images. (2) Voxel to segment: a fuzzy atlas, which indicates the probability distribution of each tissue type in the brain, is registered elastically with the MRF. It is used by the MRF as a-priori knowledge to correct miss-classified voxels. (3) Segment to segment: Remaining lesion candidates are processed by a feature based classifier that looks at unary and neighborhood information to eliminate more false positives. An expert's manual segmentation was compared with the algorithm.

In this study, subsets of MR slices were examined to assess their ability to optimally predict the total cerebral volume of gray matter, white matter and CSF. Patients underwent a clinical imaging protocol consisting of T1-, T2-, PD-, and FLAIR-weighted images after obtaining informed consent. MR imaging sets were registered, RF-corrected, and then analyzed with a hybrid neural network segmentation and classification algorithm to identify normal brain parenchyma. After processing the data, the correlation between the image subsets and the total cerebral volumes of gray matter, white matter and CSF were examined. The 29 subjects (18F, 11M) assessed in this study were 1.7 ? 18.7 (median = 5.2) years of age. The five subsets accounted for 5%, 15%, 24%, 56%, and 79% of the total cerebral volume. The predictive correlation for gray matter, white matter, and CSF in each of these subsets were: 5% (R= 0.94, 0.92, 0.91), 15% (R= 0.93, 0.95, 0.94), 24% (R= 0.92, 0.95, 0.94), 56% (R= 0.75, 0.95, 0.89), and 79% (R= 0.89, 0.98, 0.99) respectively. All subsets of slices examined were significantly correlated (p<0.001) with the total cerebral volume of gray matter, white matter, and CSF.

Splines, which were invented by Schoenberg more than fifty years ago, constitute an elegant framework for dealing with interpolation and discretization problems. They are widely used in computer-aided design and computer graphics, but have been neglected in medical imaging applications, mostly as a consequence of what one may call the bad press phenomenon. Thanks to some recent research efforts in signal processing and wavelet-related techniques, the virtues of splines have been revived in our community. There is now compelling evidence (several independent studies) that splines offer the best cost-performance tradeoff among available interpolation methods. In this presentation, we will argue that the spline representation is ideally suited for all processing tasks that require a continuous model of signals or images. We will show that most forms of spline fitting (interpolation, least squares approximation, smoothing splines) can be performed most efficiently using recursive digital filters. We will also have a look at their multiresolution properties which make them prime candidates for constructing wavelet bases and computing image pyramids. Typical application areas where these techniques can be useful are: image reconstruction from projection data, sampling grid conversion, geometric correction, visualization, rigid or elastic image registration, and feature extraction including edge detection and active contour models.

Clinical procedures that rely on biplane x-ray images for three-dimensional (3-D) information may be enhanced by three-dimensional reconstructions. However, the accuracy of reconstructed images is dependent on the uncertainty associated with the parameters that define the geometry of the camera system. In this paper, we use a numerical simulation to examine the effect of these uncertainties and to determine the limits required for adequate three-dimensional reconstruction. We then test our conclusions with images of a calibration phantom recorded using a clinical system. A set of reconstruction routines, developed for a cardiac mapping system, were used in this evaluation. The routines include procedures for correcting image distortion and for automatically locating catheter electrodes. Test images were created using a numerical simulation of a biplane x-ray projection system. The reconstruction routines were then applied using accurate and perturbed camera geometries and error maps were produced. Our results indicate that useful catheter reconstructions are possible with reasonable bounds on the uncertainty of camera geometry provided the locations of the camera isocenters are accurate. The results of this study provide a guide for the specification of camera geometry display systems and for researchers evaluating possible methodologies for determining camera geometry.

X-ray rotational angiography has recently gained increasing interest for computer-assisted quantitative analysis. It provides more accurate assessment of vascular diseases and precise inspection of complex structure of the arterial network via three-dimensional (3D) vascular reconstruction. The 3D spatial information can be obtained via a stereoscopic analysis of the two-dimensional (2D) projections of the opacified blood vessels. In this work, we focus on the problem of automatic 3D reconstruction of blood vessel networks for telediagnostic applications and therefore from uncalibrated X-ray rotational angiography image sequence. Three main issues are addressed: 1) automatic accurate subpixel vascular median axis network detection from non-subtracted 2D angiography images, 2) robust matching of the extracted features by using an original method based on statistical tests, and 3) three-dimensional reconstruction through epipolar geometry determination from uncalibrated 2D images. Our reconstruction method has the advantage to be independent of the angiography acquisition system. It is therefore interesting for telemedicine and specially for telediagnostic systems.

We have developed a technique based on epipolar geometry which allows reconstruction of vessel lumina containing asymmetries for arbitrary lumen orientation. Vessel centerlines and the imaging geometry are determined using previously described methods. Epipolar planes are established for pairs of centerline points and intensity profiles are extracted along epipolar lines. Extracted profiles are fitted to projection curves corresponding to elliptical cross-sections. The original vessel intensity profile is subtracted from the best-fit profile yielding the plaque profile. Corresponding best-fit profiles are used to reconstruct elliptical cross sections by backprojecting and identifying ray intersection points which lie within the model cross section. Assuming that the deformation of an otherwise elliptical lumen occurs as a result of plaque growing inward from the periphery, asymmetries are generated by deforming the cross sections to an extent consistent with the plaque profile in both views. The three-dimensional lumen may be obtained by combining individual cross sections in consecutive epipolar planes. The technique was evaluated using noiseless simulated angiograms. Reconstructed asymmetric lumen cross sections were found to be accurate to 95%, where accuracy was determined using area and distance criteria.

The purpose of this paper is to describe a framework for evaluating image segmentation algorithms. Image segmentation consists of object recognition and delineation. For evaluating segmentation methods, three factors - precision (reproducibility), accuracy (agreement with truth, validity), and efficiency (time taken) - need to be considered for both recognition and delineation. To assess precision, we need to choose a figure of merit, repeat segmentation considering all sources of variation, and determine variations in figure of merit via statistical analysis. It is impossible usually to establish true segmentation. Hence, to assess accuracy, we need to choose a surrogate of true segmentation and proceed as for precision. In determining accuracy, it may be important to consider different landmark areas of the structure to be segmented depending on the application. To assess efficiency, both the computational and the user time required for algorithm and operator training and for algorithm execution should be measured and analyzed. Precision, accuracy, and efficiency are interdependent. It is difficult to improve one factor without affecting others. Segmentation methods must be compared based on all three factors. The weight given to each factor depends on application.

In this work, we present a unique set of 3D MRI brain data that is appropriate for testing the intra and inter-site variability of image analysis methods. A single subject was scanned two times within a 24 hour time window each at five different MR sites over a period of six weeks using GE and Phillips 1.5 T scanners. The imaging protocol included T1 weighted, Proton Density and T2 weighted images. We applied three quantitative image analysis methods and analyzed their results via the coefficients of variability (COV) and the intra correlation coefficient. The tested methods include two multi-channel tissue segmentation techniques based on an anatomically guided manual seeding and an atlas-based seeding. The third tested method was a single-channel semi-automatic segmentation of the hippocampus. The results show that the outcome of image analysis methods varies significantly for images from different sites and scanners. With the exception of total brain volume, which shows consistent low variability across all images, the COV's were clearly larger between sites than within sites. Also, the COV's between sites with different scanner types are slightly larger than between sites with the same scanner type. The presented existence of a significant inter-site variability requires adaptations in image methods to produce repeatable measurements. This is especially of importance in multi-site clinical research.

We propose an evaluation process for segmentation which is made up of three different levels. It enables us to carry out the time consuming steps only for those segmentation methods for which a successful segmentation is foreseeable. In the first level the developer of a segmentation method does a coarse analysis of the usefulness of the individual segmentation methods by means of visual assessment of the results for few image examples. Methods which have been judged useful at the first level are investigated in a second evaluation step as to the stability of the segmentation results in case of slight deviations in the images. For the reproduction of the image formation process a multitude of realizations of a given region of interest are produced by means of the bootstrap technique. At the third level of the evaluation process the segmentation methods are tested for segmentation errors. The segmentation methods are judged by means of empirical discrepancy values, and the effectiveness of a method chosen for the respective task is finally estimated.

We describe the evaluation of a non-rigid image registration method for multi-modal data. The evaluation is made difficult by the absence of gold standard test data, for which the true transformation from one image to another is known. Different approaches have been used to deal with this deficiency, e.g., by using synthetically warped data, by comparison of anatomic regions of interest identified either manually or automatically, and by direct comparison of the registered data. Each of these approaches are limited and in this paper, we illustrate some of the problems that arise based on their application to the evaluation of our multi-modal non-rigid registration method.

For implementation of computer-aided diagnostic systems for chest radiographs, it is important to correctly identify the view position, i.e., posteroanterior (PA) or lateral view. Our purpose is to develop an advanced computerized method by using a template matching technique for correctly identifying either PA or lateral view, and to apply this method to approximately 48,000 PA and 16,000 lateral chest radiographs. To evaluate the similarity with templates, correlation values of a chest image with various templates were obtained and compared for determining whether a chest image is either a PA or a lateral view. By considering the variation in terms of patient's sizes, lung opacities, and lung sizes, we produced 24 templates of PA and lateral views. In the first step, two largest correlation values of an unknown case with 3 PA and 2 lateral templates for medium-size patients were compared for determining the view position. In the second step, the unidentifiable cases in the first step were re-examined by comparison of the correlation values with 11 to 19 templates for small and large patients. With the computerized method based on a template matching technique, 99.99%(=63,788/63,791) of chest images in the large database were correctly identified in terms of PA or lateral views.

High reproducibility of volumetric measurements is an important prerequisite for follow-up of small lung nodules in order to differentiate malignant from benign lesions in a lung cancer screening setting. This study was aimed to evaluate the measurement reproducibility of a new software tool for pulmonary nodule volumetry. In an ongoing study, 147 pulmonary nodules (size 1.6-17.5 mm) were examined with low-dose multidetector CT (Siemens Somatom Volume Zoom, 120 kVp, 20 mAs, detector collimation 4x1 mm, normalized pitch 1.75, slice thickness 1.25 mm, reconstruction increment 0.8 mm). Two consecutive low-dose scans covering the whole lung volume were performed within a few minutes. Between both scans, patients were asked to leave the CT scanner, and the second scan was planned independently from the first one. For all visually detected pulmonary nodules with a diameter <20 mm nodule volume was determined on both scans using a software prototype containing segmentation and volumetry algorithms. Results from both scans were compared. Nodule volume differences were determined as difference between the first and second measurement and ranged from 169 to 87%. The performance of the diagnostic test was measured using ROC analysis. For the detection of a volume doubling the area under curve (A_z) was 0.98, for a growth of 50% the A_z was 0.89. Further refinement of the segmentation algorithm should lead to more consistent measurements in ill-defined nodules. In conclusion, volumetric measurement of pulmonary nodules in multislice CT data sets is a reliable tool for the detection of growth in small pulmonary nodules.

Ultrasonic strain images describe the stiffness of soft tissues. Strain estimates are obtained from the spatial derivative of local displacements between echo fields acquired prior to and after applying a small compressive force. Multi-scale displacement measurements using correlation estimates yield the highest quality strain images but are computationally extensive. An approach to strain reconstruction is proposed that constrains the correlation search based on physical priors in order to accelerate and optimize the estimation of local displacements. The data correlation kernel is chosen large in regions of constant displacement to minimize noise and small in regions of varying displacement to maximize spatial resolution for strain. Multiple echo fields were acquired using a Siemens Elegra system with a 7.5 MHz linear array while slowly compressing phantoms and tissues. Gel phantoms with simulated stiff lesions and flow channels, as well as ex vivo muscle and in vivo breast tissues were examined. The new algorithm provided images with equal or lower noise as compared to the traditional algorithm. Adaptively limiting the search in smoothly compressed regions reduced computational time by a factor of 1.5 to 8 depending on the applied compression and complexity of motion.

Correction of motion artifacts in MRI due to interview in-plane 2-D rigid body translations is possible using only the raw data of two standard 2DFT images acquired of the same object with phase encode direction swapped. Previous techniques simply use multiple or orthogonal images to reduce artifacts by ghost interference or geometric averaging. The orthogonal k-space phase difference (ORKPHAD) provides an overdetermined system of linear equations that can be solved directly to compensate for the phase errors caused by the translation. This technique was used to correct images of a volunteer's moving wrist. For all slices of the motion corrupted data volumes, artifacts were dramatically reduced. Though the algorithm only accounts for interview translational motion, it is robust enough to correct real in vivo images that may also be corrupted by small amounts of rotational or out-of-plane motion. Experiments underway show the algorithm can tolerate bulk rotation of several degrees between orthogonal image pairs and that corrections using fractional NEX scans are possible. The current results and such ongoing advancements should make this correction technique practical for certain clinical scenarios vulnerable to in-plane translation.

In many medical imaging applications, due to the limited field of view of imaging devices, often, acquired images include only a part of a structure. In such situations, it is impossible to guarantee that the images will contain exactly the same physical extent of the structure at different scans, which leads to difficulties in registration and in many other tasks, such as the analysis of the morphology, architecture, and kinematics. To facilitate such analysis, we developed a general method, referred to as isoshaping, that generates structures of the same shape from segmented images. The basis for this method is to automatically find a set of key points, called shape centers, in the segmented partial anatomic structure such that these points are present in all images and that they represent the same physical location in the object, and then trim the structure using these points as reference. The application area considered here is the analysis of the morphology, architecture, and kinematics of the foot joints from MR images acquired at different joint positions, load conditions, and longitudinal time instances. The accuracy of the method is studied and it is quantitatively demonstrated that isoshaping improves the results of registration.

In the last decade in the field of diagnostic imaging, the problem of variability of reader skill as well as patient case difficulty has given rise to a multivariate approach to receiver operating characteristic (ROC) analysis. The multivariate approach is the so-called multiple-reader, multiple-case (MRMC) ROC paradigm in which every reader reads every patient case and, where possible, in each of two modalities under comparison. The present paper demonstrates the isomorphism between patient cases, image readers, and imaging modalities in diagnostic imaging and, respectively, test sets, training sets, and competing discriminant algorithms in the field of statistical pattern recognition (SPR). Thus, the MRMC paradigm can be brought directly across from imaging to SPR. Recent MRMC ROC analytical methods are demonstrated in the context of SPR for the task of analyzing the natural components of variance in that problem involving test sets, training sets, and competing discriminants. Monte Carlo trials reported here indicate that the conventional wisdom that the variance of measures of classifier accuracy comes mainly from the finite test set is only true when assessing a single algorithm in a very limited context. In particular, it is generally not true when comparing competing discriminant algorithms; in that case the variance is dominated by the finite training set.

The purpose of this study was to identify and characterize clusters in a heterogeneous breast cancer computer-aided diagnosis database. Identification of subgroups within the database could help elucidate clinical trends and facilitate future model building. Agglomerative hierarchical clustering and k-means clustering were used to identify clusters in a large, heterogeneous computer-aided diagnosis database based on mammographic findings (BI-RADS) and patient age. The clusters were examined in terms of their feature distributions. The clusters showed logical separation of distinct clinical subtypes such as architectural distortions, masses, and calcifications. Moreover, the common subtypes of masses and calcifications were stratified into clusters based on age groupings. The percent of the cases that were malignant was notably different among the clusters. Cluster analysis can provide a powerful tool in discerning the subgroups present in a large, heterogeneous computer-aided diagnosis database.

Statistical as well as adaptive clustering approaches are being currently used for both segmentation and vector quantization of medical images. However, a comparative evaluation of both approaches has rarely been done to identify the efficacy of such approaches to specific applications, for example, image segmentation and vector quantization. The rate distortion functions of three clustering algorithms, namely, the statistical based deterministic annealing, the adaptive fuzzy leader clustering algorithm, and LBG, have been computed for vector quantization using multi-scale vectors in the wavelet domain. Such comparative evaluation serves as a guide for proper selection of clustering algorithms for global codebook generation in vector quantization and for image segmentation.

We tested on an edge map computed from a local homogeneity measurement, which is a potential replacement for the traditional gradient-based edge map in level-set segmentation. In existing level-set methods, the gradient information is used as a stopping criteria for curve evolution, and also provides the attracting force to the zero level-set from the target boundary. However, in a discrete implementation, the gradient-based term can never fully stop the level-set evolution even for ideal edges, leakage is often unavoidable. Also the effective distance of the attracting force and blurring of edges become a trade-off in choosing the shape and support of the smoothing filter. The proposed homogeneity measurement provides easier and more robust edge estimation, and the possibility of fully stopping the level-set evolution. The homogeneity term decreasing from a homogenous region to the boundary, which dramatically increases the effective distance of the attracting force and also provides additional measurement of the overall approximation to the target boundary. Therefore, it provides a reliable criteria of adaptively changing the advent speed. By using this term, the leakage problem was avoided effectively in most cases compared to traditional level-set methods. The computation of the homogeneity is fast and its extension to the 3D case is straightforward.

Image segmentation is concerned with partitioning an image into non-overlapping, constituent regions, which are homogeneous with respect to certain features. In magnetic resonance imaging (MRI), the most discriminative and commonly used features are the image intensities themselves. However, due to noise, partial volume effects, natural and spurious intensity variations, intensity distributions of distinct tissues generally overlap, which makes segmentation difficult and less precise. Using multi-spectral MR images and mapping intensities into a multidimensional feature space may help in segmentation. To further facilitate segmentation, we map the intensities and second derivatives of multi-spectral images into a common multidimensional feature space. Integration of intensity and spatial information may yield complex clusters that cannot be described by Gaussian mixture models or by hyper-spherical shapes. For this reason we devise a novel segmentation method based on non-parametric valley-seeking clustering. The valleys are found by estimating feature density gradients. The proposed segmentation method, with and without spatial information, is tested on simulated and real, single- and multi-spectral, MR brain images. The segmentation results are highly consistent with the gold standard, especially when combined with a procedure for intensity non-uniformity correction, presented in MI 4684-177.

High-resolution multichannel textures are difficult to characterize with simple statistics and the high level of detail makes the selection of a particular contour using classical gradient-based methods not effective. We have developed a hybrid method that combines fuzzy connectedness and Voronoi diagram classification for the segmentation of color and multichannel objects. The multi-step classification process relies on homogeneity measures derived from moment statistics and histogram information. These color features have been optimized to best combine individual channel information in the classification process. The segmentation initialization requires only a set of interior and exterior seed points, minimizing user intervention and the influence of the initialization on the overall quality of the results. The method was tested on volumes from the Visible Human and on brain multi-protocol MRI data sets. The hybrid segmentation produced robust, rapid and finely detailed contours with good visual accuracy. The addition of quantized statistics and color histogram distances as classification features improved the robustness of the method with regards to initialization when compared to our original implementation.

An automatic segmentation method for detecting the prostate boundary in transrectal ultrasound (TRUS) images was developed. The TRUS images were preprocessed by using an adaptive directional filtering and an automatic attenuation compensation for noise removal and contrast enhancement. A directional search strategy was used to locate key-points on the prostate boundary. The prostate contour was interpolated from the key-points under the supervision of a morphological prostate boundary model, which had been trained from prior manual segmentation of a large number of TRUS images. A new prostate center was calculated based on the intermediate segmentation result. The algorithm is reiterated until the prostate boundary and center reach a stable state. The overall performance of the method was compared to manual segmentation of an expert radiologist. About 78% out of 282 TRUS images (excluding base and apex slices) from three types of ultrasound machine (Acuson, Siemens, and B&K) were correctly delineated. The segmentation error was 0.9 mm averaged on 30 selected images, 10 for each type of machine. The computation time for a typical series of TRUS images is approximately 1 minute on a Pentium-II computer.

This paper presents a three-dimensional (3-D) tissue analysis method and its applications in partial volume correction and change analysis. The method uses a stochastic model-based approach and consists of two steps: (1) unsupervised tissue quantification and (2) 3-D segmentation. Firstly, the MR image volume is modeled by the standard finite normal mixture (SFNM) distribution. It has been shown that the SFNM converges to the true distribution when the pixel images are asymptotically independent. Secondly, the tissue quantification is achieved through (1) model selection by minimum description length (MDL) criterion; (2) parameter initialization by optimal histogram quantization and (3) parameter estimation by a fast EM algorithm using the global 3-D histogram rather than conventionally the raw data. Finally, we develop a 3-D segmentation method using the maximum likelihood (ML) classification and contextual Bayesian relaxation labeling (CBRL). The CBRL is developed to obtain a consistent labeling solution, based on localized SFNM formulation by using neighborhood contextual regularities. The method has been applied to partial volume correction for PET brain images and change analysis for MR breast images.

Active Appearance Models (AAMs) are useful for the segmentation of cardiac MR images since they exploit prior knowledge about the cardiac shape and image appearance. However, traditional AAMs only process 2D images, not taking into account the 3D data inherent to MR. This paper presents a novel, true 3D Active Appearance Model that models the intrinsic 3D shape and image appearance of the left ventricle in cardiac MR data. In 3D-AAM, shape and appearance of the Left Ventricle (LV) is modeled from a set of expert drawn contours. The contours are then resampled to a manually defined set of landmark points, and subsequently aligned. Appearance variations in both shape and texture are captured using Principal Component Analysis (PCA) on the training set. Segmentation is achieved by minimizing the model appearance-to-target differences by adjusting the model eigen-coefficients using a gradient descent approach. The clinical potential of the 3D-AAM is demonstrated in short-axis cardiac magnetic resonance (MR) images. The method's performance was assessed by comparison with manually-identified independent standards in 56 clinical MR sequences. The method showed good agreement with the independent standards using quantitative indices such as border positioning errors, endo- and epicardial volumes, and left ventricular mass. The 3D AAM method shows high promise for successful segmentation of three-dimensional images in MR.

In this paper we propose an appearance-based method to heart detection by principal component analysis (PCA). In contrast to conventional methods of PCA-based training, there is no brightness and contrast normalization since such normalization is usually based on maximum and minimum intensity values and is very sensitive to noises. We propose to integrate the normalization procedure into the detection phase. This is achieved by projecting the intensity-transformed image (with unknown scale and shift parameters) onto the eigen-images and minimizing the error of fit. This leads to a set of equations on both the intensity transformation parameters and the projection coefficients. By using the least-squares method, these equations can be easily solved for the scale and shift parameters. After an initial detection of heart positions is conducted, robust fitting of the heart trajectory is used to correct any detection errors. Besides, we also propose an eigen-image re-orthonormalization method for multiple resolution detection without extra training on multiple scales.

A novel 3-D Active Appearance Model (3-D AAM) is applied to fully automated endocardial contour detection in 2-D + time (2DT) 4-chamber ultrasound sequences, without knowledge of cardiac phase (ED/ES frames). 2DT appearance of the heart is modeled in 3-D by converting the stack of 2-D time slices into a 3-D voxel space. In a training set, an expert defines corresponding endocardial contour points for one complete cardiac cycle (ED to ED). 2DT shape is represented as a 3-D surface. Image appearance is modeled as a vector of voxel intensities in a volume-patch spanned by the 3-D surface. Principal Component Analysis extracts eigenvariations of 3-D shape and appearance, capturing typical cardiac motion patterns. 3-D AAM segments the image volume by minimizing 3-D model-to-target intensity differences, adjusting eigenvariation coefficients and 3-D pose using gradient descent minimization. This provides time-continuous border localization for one beat in both time and space. The method was used on 3-beat sequences from 129 patients split randomly into a training (65) and a test set (64). An independent expert manually drew all endocardial contours. 3-D AAM converged well in 89% of test cases. Average absolute temporal error was 37.0 msec, spatial error 3.35 mm, comparable to human inter-observer variabilities.

An automated method for the segmentation of thrombus in abdominal aortic aneurysms from CTA data is presented. The method is based on Active Shape Model (ASM) fitting in sequential slices, using the contour obtained in one slice as the initialisation in the adjacent slice. The optimal fit is defined by maximum correlation of grey value profiles around the contour in successive slices, in contrast to the original ASM scheme as proposed by Cootes and Taylor, where the correlation with profiles from training data is maximised. An extension to the proposed approach prevents the inclusion of low-intensity tissue and allows the model to refine to nearby edges. The applied shape models contain either one or two image slices, the latter explicitly restricting the shape change from slice to slice. To evaluate the proposed methods a leave-one-out experiment was performed, using six datasets containing 274 slices to segment. Both adapted ASM schemes yield significantly better results than the original scheme (p<0.0001). The extended slice correlation fit of a one-slice model showed best overall performance. Using one manually delineated image slice as a reference, on average a number of 29 slices could be automatically segmented with an accuracy within the bounds of manual inter-observer variability.

Knowledge about the location of the diaphragm dome surface, which separates the lungs and the heart from the abdominal cavity, is of vital importance for applications like automated segmentation of adjacent organs (e.g., liver) or functional analysis of the respiratory cycle. We present a new 3D Active Appearance Model (AAM) approach to segmentation of the top layer of the diaphragm dome. The 3D AAM consists of three parts: a 2D closed curve (reference curve), an elevation image and texture layers. The first two parts combined represent 3D shape information and the third part image intensity of the diaphragm dome and the surrounding layers. Differences in height between dome voxels and a reference plane are stored in the elevation image. The reference curve is generated by a parallel projection of the diaphragm dome outline in the axial direction. Landmark point placement is only done on the (2D) reference curve, which can be seen as the bounding curve of the elevation image. Matching is based on a gradient-descent optimization process and uses image intensity appearance around the actual dome shape. Results achieved in 60 computer generated phantom data sets show a high degree of accuracy (positioning error -0.07+/-1.29 mm). Validation using real CT data sets yielded a positioning error of -0.16+/-2.95 mm. Additional training and testing on in-vivo CT image data is ongoing.

A 3D surface model for the segmentation of organic structures in CT-and MR-datasets has been developed. The training dataset for the surface model is computed from semiautomatically generated voxelsets. Triangulated meshes of the voxelsets representing the objects' surfaces are generated. The surface model is able to learn the shape variations in the training dataset by a principal component analysis of the information provided by the points forming the triangulated surface-meshes. Furthermore, the image information at the mesh points is also added into gray value models describing the gray value distribution at this particular surface section. The optimization of the model is performed by iteratively moving the surface points of the model towards image structures fitting to the gray value models. During the optimization process the models shape information assures that the surface stays plausible. A 3D-model of the spleen consisting of 10 objects, and a kidney model generated from 7 left kidneys have been developed. The models have been tested on 3 unknown spleen- and 3 unknown kidney-datasets. The total cover between the model and the organs varied between 65% and 75%, which is a respectable result in the face of the small training datasets the models were generated from.

The aim of this work is to explore the performance of active appearance models (AAMs) in reconstruction and interpretation of bones in hand radiographs. AAM is a generative approach that unifies image segmentation and image understanding. Initial locations for the AAM search are generated by an exhaustive filtering method. A series of AAMs for smaller groups of bones are used. It is found that AAM successful reconstructs 99% of metacarpals, proximal and medial phalanges and the distal 3 cm of radius and ulna. The rms accuracy is better than 240 microns (point-to-curve). The generative property is used (1) to define a measure of fit that allows the models to self-evaluate and chose between the multiple found solutions, (2) to overcome obstacles in the image in the form of rings by predicting the missing part, and (3) to detect anomalies, e.g. rheumatoid arthritis. The shape scores are used as biometrics to check the identity of patients in a longitudinal study. The conclusion is that AAM provides a highly efficient and unified framework for various tasks in diagnosis and assessment of bone related disorders.

This paper describes a novel approach to register 3D computed tomography (CT) or magnetic resonance (MR) images to a set of 2D X-ray images. Such a registration may be a valuable tool for intraoperative determination of the precise position and orientation of some anatomy of interest, defined in preoperative images. The registration is based solely on the information present in 2D and 3D images. It does not require fiducial markers, X-ray image segmentation, or construction of digitally reconstructed radiographs. The originality of the approach is in using normals to bone surfaces, preoperatively defined in 3D MR or CT data, and gradients of intraoperative X-ray images, which are back-projected towards the X-ray source. The registration is then concerned with finding that rigid transformation of a CT or MR volume, which provides the best match between surface normals and back projected gradients, considering their amplitudes and orientations. The method is tested on a lumbar spine phantom. Gold standard registration is obtained by fidicual markers attached to the phantom. Volumes of interest, containing single vertebrae, are registered to different pairs of X-ray images from different starting positions, chosen randomly and uniformly around the gold standard position. Target registration errors and rotation errors are in order of 0.3 mm and 0.35 degrees for the CT to X-ray registration and 1.3 mm and 1.5 degrees for MR to X-ray registration. The registration is shown to be fast and accurate.

This paper focuses on the problem of ill-posedness of deformable point set registration and we propose a new approach to restrict the solution space using shape information. The basic elements of the investigated kind of registration algorithm are a cost functional, an optimization strategy and a motion model. The motion model determines the kind of motions and deformations that are allowed and how they are restricted. The motion model itself is mainly determined by the kind of parameterized transformation used to express the motion/deformation. Here, we observe that matching with more degrees of freedom (the parameters of the transformation) than necessary can introduce mismatches due to a higher sensitivity to noise or by destroying local shape information. In this paper we propose a cost functional which is robust to noise and we introduce a new method to specify a shape adapted deformation model based on thin-plate splines and initial control point placing using point clustering. We show that these initial positions have a strong impact on the match and we define them as cluster centers where we cluster on one of the point sets (weighting each point of this set with its distance to the other point set). Our experiments with known ground truth show that the shape adapted model recovers constantly very accurately corresponding points. In our evaluation with more than 1200 single experiments we showed that, compared to a conventional octree based scheme, we could save more than 60% of degrees of freedom while preserving matching quality.

We are investigating interventional MRI guided radio- frequency (RF) thermal ablation for the minimally invasive treatment of prostate cancer. Among many potential applications of registration, we wish to compare registered MR images acquired before and immediately after RF ablation in order to determine whether a tumor is adequately treated. Warping registration is desired to correct for potential deformations of the pelvic region and movement of the prostate. We created a two-step, three-dimensional (3D) registration algorithm using mutual information and thin plate spline (TPS) warping for MR images. First, automatic rigid body registration was used to capture the global transformation. Second, local warping registration was applied. Interactively placed control points were automatically optimized by maximizing the mutual information of corresponding voxels in small volumes of interest and by using a 3D TPS to express the deformation throughout the image volume. Images were acquired from healthy volunteers in different conditions simulating potential applications. A variety of evaluation methods showed that warping consistently improved registration for volume pairs whenever patient position or condition was purposely changed between acquisitions. A TPS transformation based on 180 control points generated excellent warping throughout the pelvis following rigid body registration. The prostate centroid displacement for a typical volume pair was reduced from 3.4 mm to 0.6 mm when warping was added.

The generation of patient specific meshes for Finite Element Methods (FEM) with application to brain deformation is a time consuming process, but is essential for the modeling of intra-operative deformation of the brain during neurosurgery using FEM techniques. We present an automatic method for the generation of FEM meshes fitting patient data. The method is based on non-rigid registration of patient MR images to an atlas brain image, followed by deformation of a high-quality mesh of this atlas brain. We demonstrate the technique on brain MRI images from 12 patients undergoing neurosurgery. We show that the FEM meshes generated by our technique are of good quality. We then demonstrate the utility of these FEM meshes by simulating simple neuro-surgical scenarios on example patients, and show that the deformations predicted by our brain model match the observed deformations. The meshes generated by our technique are of good quality, and are suitable for modeling the types of deformation observed during neurosurgery. The deformations predicted by a simple loading scenario match well with those observed following the actual surgery. This paper does not attempt an exhaustive study of brain deformation, but does provide an essential tool for such a study - a method of rapidly generating Finite Element Meshes fitting individual subject brains.