We address the problem of learning a binary classifier given partially labeled data where all labeled samples come from only one of the classes, commonly known as Positive Unlabeled (PU) learning. Classical methods such as clustering, out-of-distribution detection, and positive density estimation, while effective in low-dimensional scenarios, lose their efficacy as the dimensionality of data increases, because of the increasing complexity. This has led to the development of methods that address the problem in high-dimensional spaces; however, many of these methods are also impacted by the increased complexity inherent in high-dimensional data. The contribution of this paper is the learning of a neural network-based data representation by employing a loss function that enables the projection of unlabeled data into two distinct clusters - positive and negative - facilitating their identification through basic clustering techniques and mirroring the simplicity of the problem seen in low-dimensional settings. We further enhance this separation of unlabeled data clusters by implementing a vector quantization strategy. Our experimental results on benchmarking PU datasets validate the superiority of our method over existing state-of-the-art techniques. Additionally, we provide theoretical justification to support our cluster-based approach and algorithmic choices.

Segmentation is a critical step in analyzing the developing human fetal brain. There have been vast improvements in automatic segmentation methods in the past several years, and the Fetal Brain Tissue Annotation (FeTA) Challenge 2021 helped to establish an excellent standard of fetal brain segmentation. However, FeTA 2021 was a single center study, limiting real-world clinical applicability and acceptance. The multi-center FeTA Challenge 2022 focused on advancing the generalizability of fetal brain segmentation algorithms for magnetic resonance imaging (MRI). In FeTA 2022, the training dataset contained images and corresponding manually annotated multi-class labels from two imaging centers, and the testing data contained images from these two centers as well as two additional unseen centers. The multi-center data included different MR scanners, imaging parameters, and fetal brain super-resolution algorithms applied. 16 teams participated and 17 algorithms were evaluated. Here, the challenge results are presented, focusing on the generalizability of the submissions. Both in- and out-of-domain, the white matter and ventricles were segmented with the highest accuracy (Top Dice scores: 0.89, 0.87 respectively), while the most challenging structure remains the grey matter (Top Dice score: 0.75) due to anatomical complexity. The top 5 average Dices scores ranged from 0.81-0.82, the top 5 average percentile Hausdorff distance values ranged from 2.3-2.5mm, and the top 5 volumetric similarity scores ranged from 0.90-0.92. The FeTA Challenge 2022 was able to successfully evaluate and advance generalizability of multi-class fetal brain tissue segmentation algorithms for MRI and it continues to benchmark new algorithms.

Quantile Regression (QR) can be used to estimate aleatoric uncertainty in deep neural networks and can generate prediction intervals. Quantifying uncertainty is particularly important in critical applications such as clinical diagnosis, where a realistic assessment of uncertainty is essential in determining disease status and planning the appropriate treatment. The most common application of quantile regression models is in cases where the parametric likelihood cannot be specified. Although quantile regression is quite robust to outlier response observations, it can be sensitive to outlier covariate observations (features). Outlier features can compromise the performance of deep learning regression problems such as style translation, image reconstruction, and deep anomaly detection, potentially leading to misleading conclusions. To address this problem, we propose a robust solution for quantile regression that incorporates concepts from robust divergence. We compare the performance of our proposed method with (i) least trimmed quantile regression and (ii) robust regression based on the regularization of case-specific parameters in a simple real dataset in the presence of outlier. These methods have not been applied in a deep learning framework. We also demonstrate the applicability of the proposed method by applying it to a medical imaging translation task using diffusion models.

In the case of learning from positive and unlabeled data, the input data consist of (1) observations from the positive class and their corresponding labels and (2) unlabeled observations from both positive and negative classes. Generative Adversarial Networks (GANs) have been used to reduce the problem to the supervised setting with the advantage that supervised learning has state-of-the-art accuracy in classification tasks. In order to generate pseudo-negative observations, GANs are trained on positive and unlabeled observations with a modified loss. Using both positive and pseudo-negative observations leads to a supervised learning setting. The generation of pseudo-negative observations that are realistic enough to replace missing negative class samples is a bottleneck for current GAN-based algorithms. By including an additional classifier into the GAN architecture, we describe a novel GAN-based approach. The GAN discriminator instructs the generator to only produce samples that fall into the unlabeled data distribution, while a second classifier (observer) network monitors the GAN training to: (i) prevent the generated samples from falling into the positive distribution; and (ii) learn the features that are the key distinction between the positive and negative observations. Experiments on four image datasets demonstrate that our trained observer network performs better than existing techniques in discriminating between real unseen positive and negative samples.

Partial correlations (PCs) of functional magnetic resonance imaging (fMRI) time series play a principal role in revealing connectivity of brain networks. To explore nonlinear behavior of the blood-oxygen-level dependent signal, the present work postulates a kernel-based nonlinear connectivity model based on which it obtains topology revealing PCs. Instead of relying on a single predefined kernel, a data-driven approach is advocated to learn the combination of multiple kernel functions that optimizes the data fit. Synthetically generated data based on both a dynamic causal and a linear model are used to validate the proposed approach in resting-state fMRI scenarios, highlighting the gains in edge detection performance when compared with the popular linear PC method. Tests on real fMRI data demonstrate that connectivity patterns revealed by linear and nonlinear models are different.

We present a method that quantifies point-wise changes in surface morphology of the bones of the human wrist. The proposed method, referred to as Registration-based Bone Morphometry (RBM), consists of two steps: an atlas selection step and an atlas warping step. The atlas for individual wrist bones was selected based on the shortest ℓ2 distance to the ensemble of wrist bones from a database of a healthy population of subjects. The selected atlas was then warped to the corresponding bones of individuals in the population using a non-linear registration method based on regularized ℓ2 distance minimization. The displacement field thus calculated showed local differences in bone shape that then were used for the analysis of group differences. Our results indicate that RBM has potential to provide a standardized approach to shape analysis of bones of the human wrist. We demonstrate the performance of RBM for examining group differences in wrist bone shapes based on sex and between those of the right and left wrists in healthy individuals. We also present data to show the application of RBM for tracking bone erosion status in rheumatoid arthritis.

The estimation and analysis of kinetic parameters in dynamic positron emission tomography (PET) is frequently confounded by tissue heterogeneity and partial volume effects. We propose a new constrained model of dynamic PET to address these limitations. The proposed formulation incorporates an explicit mixture model in which each image voxel is represented as a mixture of different pure tissue types with distinct temporal dynamics. We use Cramér-Rao lower bounds to demonstrate that the use of prior information is important to stabilize parameter estimation with this model. As a result, we propose a constrained formulation of the estimation problem that we solve using a two-stage algorithm. In the first stage, a sparse signal processing method is applied to estimate the rate parameters for the different tissue compartments from the noisy PET time series. In the second stage, tissue fractions and the linear parameters of different time activity curves are estimated using a combination of spatial-regularity and fractional mixture constraints. A block coordinate descent algorithm is combined with a manifold search to robustly estimate these parameters. The method is evaluated with both simulated and experimental dynamic PET data.

We investigate using dual time-point PET data to perform Patlak modeling. This approach can be used for whole body dynamic PET studies in which we compute voxel-wise estimates of Patlak parameters using two frames of data for each bed position. Our approach directly uses list-mode arrival times for each event to estimate the Patlak parametric image. We use a penalized likelihood method in which the penalty function uses spatially variant weighting to ensure a count independent local impulse response. We evaluate performance of the method in comparison to fractional changes in SUV values (%DSUV) between the two frames using Cramer Rao analysis and Monte Carlo simulation. Receiver operating characteristic (ROC) curves are used to compare performance in differentiating tumors relative to background based on the dynamic data sets. Using area under the ROC curve as a performance metric, we show superior performance of Patlak relative to %DSUV over a range of dynamic data sets and parameters. These results suggest that Patlak analysis may be appropriate for analysis of dual time-point whole body PET data and could lead to superior detection of tumors relative to %DSUV metrics.

Time-of-flight (TOF) positron emission tomography (PET) scanners offer the potential for significantly improved signal-to-noise ratio (SNR) and lesion detectability in clinical PET. However, fully 3D TOF PET image reconstruction is a challenging task due to the huge data size. One solution to this problem is to rebin TOF data into a lower dimensional format. We have recently developed Fourier rebinning methods for mapping TOF data into non-TOF formats that retain substantial SNR advantages relative to sinograms acquired without TOF information. However, mappings for rebinning into non-TOF formats are not unique and optimization of rebinning methods has not been widely investigated. In this paper we address the question of optimal rebinning in order to make full use of TOF information. We focus on FORET-3D, which approximately rebins 3D TOF data into 3D non-TOF sinogram formats without requiring a Fourier transform in the axial direction. We optimize the weighting for FORET-3D to minimize the variance, resulting in H(2)-weighted FORET-3D, which turns out to be the best linear unbiased estimator (BLUE) under reasonable approximations and furthermore the uniformly minimum variance unbiased (UMVU) estimator under Gaussian noise assumptions. This implies that any information loss due to optimal rebinning is as a result only of the approximations used in deriving the rebinning equation and developing the optimal weighting. We demonstrate using simulated and real phantom TOF data that the optimal rebinning method achieves variance reduction and contrast recovery improvement compared to nonoptimized rebinning weightings. In our preliminary study using a simplified simulation setup, the performance of the optimal rebinning method was comparable to that of fully 3D TOF MAP.

We describe a nonparametric framework for incorporating information from co-registered anatomical images into positron emission tomographic (PET) image reconstruction through priors based on information theoretic similarity measures. We compare and evaluate the use of mutual information (MI) and joint entropy (JE) between feature vectors extracted from the anatomical and PET images as priors in PET reconstruction. Scale-space theory provides a framework for the analysis of images at different levels of detail, and we use this approach to define feature vectors that emphasize prominent boundaries in the anatomical and functional images, and attach less importance to detail and noise that is less likely to be correlated in the two images. Through simulations that model the best case scenario of perfect agreement between the anatomical and functional images, and a more realistic situation with a real magnetic resonance image and a PET phantom that has partial volumes and a smooth variation of intensities, we evaluate the performance of MI and JE based priors in comparison to a Gaussian quadratic prior, which does not use any anatomical information. We also apply this method to clinical brain scan data using F(18) Fallypride, a tracer that binds to dopamine receptors and therefore localizes mainly in the striatum. We present an efficient method of computing these priors and their derivatives based on fast Fourier transforms that reduce the complexity of their convolution-like expressions. Our results indicate that while sensitive to initialization and choice of hyperparameters, information theoretic priors can reconstruct images with higher contrast and superior quantitation than quadratic priors.

Manually labeled landmark sets are often required as inputs for landmark-based image registration. Identifying an optimal subset of landmarks from a training dataset may be useful in reducing the labor intensive task of manual labeling. In this paper, we present a new problem and a method to solve it: given a set of N landmarks, find the k(< N) best landmarks such that aligning these k landmarks that produce the best overall alignment of all N landmarks. The resulting procedure allows us to select a reduced number of landmarks to be labeled as a part of the registration procedure. We apply this methodology to the problem of registering cerebral cortical surfaces extracted from MRI data. We use manually traced sulcal curves as landmarks in performing inter-subject registration of these surfaces. To minimize the error metric, we analyze the correlation structure of the sulcal errors in the landmark points by modeling them as a multivariate Gaussian process. Selection of the optimal subset of sulcal curves is performed by computing the error variance for the subset of unconstrained landmarks conditioned on the constrained set. We show that the registration error predicted by our method closely matches the actual registration error. The method determines optimal curve subsets of any given size with minimal registration error.

Meiotic mapping of quantitative trait loci regulating expression (eQTLs) has allowed the construction of gene networks. However, the limited mapping resolution of these studies has meant that genotype data are largely ignored, leading to undirected networks that fail to capture regulatory hierarchies. Here we use high resolution mapping of copy number eQTLs (ceQTLs) in a mouse-hamster radiation hybrid (RH) panel to construct directed genetic networks in the mammalian cell. The RH network covering 20,145 mouse genes had significant overlap with, and similar topological structures to, existing biological networks. Upregulated edges in the RH network had significantly more overlap than downregulated. This suggests repressive relationships between genes are missed by existing approaches, perhaps because the corresponding proteins are not present in the cell at the same time and therefore unlikely to interact. Gene essentiality was positively correlated with connectivity and betweenness centrality in the RH network, strengthening the centrality-lethality principle in mammals. Consistent with their regulatory role, transcription factors had significantly more outgoing edges (regulating) than incoming (regulated) in the RH network, a feature hidden by conventional undirected networks. Directed RH genetic networks thus showed concordance with pre-existing networks while also yielding information inaccessible to current undirected approaches.

Neuroimaging data, such as 3-D maps of cortical thickness or neural activation, can often be analyzed more informatively with respect to the cortical surface rather than the entire volume of the brain. Any cortical surface-based analysis should be carried out using computations in the intrinsic geometry of the surface rather than using the metric of the ambient 3-D space. We present parameterization-based numerical methods for performing isotropic and anisotropic filtering on triangulated surface geometries. In contrast to existing FEM-based methods for triangulated geometries, our approach accounts for the metric of the surface. In order to discretize and numerically compute the isotropic and anisotropic geometric operators, we first parameterize the surface using a p-harmonic mapping. We then use this parameterization as our computational domain and account for the surface metric while carrying out isotropic and anisotropic filtering. To validate our method, we compare our numerical results to the analytical expression for isotropic diffusion on a spherical surface. We apply these methods to smoothing of mean curvature maps on the cortical surface, a step commonly required for analysis of gyrification or for registering surface-based maps across subjects.

Detection of small lesions in fluorodeoxyglucose (FDG) positron emission tomography (PET) is limited by image resolution and low signal to noise ratio. We have previously described a matched subspace detection method that uses the time activity curve to distinguish tumors from background in dynamic FDG PET. Applying this algorithm on a voxel by voxel basis throughout the dynamic image produces a test statistic image or "map" which on thresholding indicates the potential locations of secondary or metastatic tumors. In this paper, we describe a thresholding method that controls familywise error rate (FWER) for the matched subspace detection statistical map. The method involves three steps. First, the PET image is segmented into several homogeneous regions. Then, the statistical map is normalized to a zero mean unit variance Gaussian random field. Finally, the images are thresholded at a fixed FWER by estimating their spatial smoothness and applying a random field theory maximum statistic approach. We evaluate this thresholding method using digital phantoms generated from clinical dynamic images. We also present an application of the proposed approach to clinical PET data from a breast cancer patient with metastatic disease.

We describe a matched subspace detection algorithm to assist in the detection of small tumors in dynamic positron emission tomography (PET) images. The algorithm is designed to differentiate tumors from background using the time activity curves (TACs) that characterize the uptake of PET tracers. TACs are modeled using linear subspaces with additive Gaussian noise. Using TACs from a primary tumor region of interest (ROI) and one or more background ROIs, each identified by a human observer, two linear subspaces are identified. Applying a matched subspace detector to these identified subspaces on a voxel-by-voxel basis throughout the dynamic image produces a test statistic at each voxel which on thresholding indicates potential locations of secondary or metastatic tumors. The detector is derived for three cases: using a single TAC with white noise of unknown variance, using a single TAC with known noise covariance, and detection using multiple TACs within a small ROI with known noise covariance. The noise covariance is estimated for the reconstructed image from the observed sinogram data. To evaluate the proposed method, a simulation-based receiver operating characteristic (ROC) study for dynamic PET tumor detection is designed. The detector uses a dynamic sequence of frame-by-frame 2-D reconstructions as input. We compare the performance of the subspace detectors with that of a Hotelling observer applied to a single frame image and of the Patlak method applied to the dynamic data. We also show examples of the application of each detection approach to clinical PET data from a breast cancer patient with metastatic disease.

We present accurate and efficient methods for estimating the spatial resolution and noise properties of nonquadratically regularized image reconstruction for positron emission tomography (PET). It is well known that quadratic regularization tends to over-smooth sharp edges. Many types of edge-preserving nonquadratic penalties have been proposed to overcome this problem. However, there has been little research on the quantitative analysis of nonquadratic regularization due to its nonlinearity. In contrast, quadratically regularized estimators are approximately linear and are well understood in terms of resolution and variance properties. We derive new approximate expressions for the linearized local perturbation response (LLPR) and variance using the Taylor expansion with the remainder term. Although the expressions are implicit, we can use them to accurately predict resolution and variance for nonquadratic regularization where the conventional expressions based on the first-order Taylor truncation fail. They also motivate us to extend the use of a certainty-based modified penalty to nonquadratic regularization cases in order to achieve spatially uniform perturbation responses, analogous to uniform spatial resolution in quadratic regularization. Finally, we develop computationally efficient methods for predicting resolution and variance of nonquadratically regularized reconstruction and present simulations that illustrate the validity of these methods.

We describe a fast forward and back projector pair based on inverse Fourier rebinning for use in iterative image reconstruction for fully 3-D positron emission tomography (PET). The projector pair is used as part of a factored system matrix that takes into account detector-pair response by using shift-variant sinogram blur kernels, thereby combining the computational advantages of Fourier rebinning with iterative reconstruction using accurate system models. The forward projector consists of a 2-D projector, which maps 3-D images into 2-D direct sinograms, followed by exact inverse rebinning which maps the 2-D into fully 3-D sinograms. The back projector is implemented as the transpose of the forward projector and differs from the true exact rebinning operator in the sense that it does not require reprojection to compute missing lines of response (LORs). We compensate for two types of inaccuracies that arise in a cylindrical PET scanner when using inverse Fourier rebinning: 1) nonuniform radial sampling and 2) nonconstant oblique angles in the radial direction in a single oblique sinogram. We examine the effects of these corrections on sinogram accuracy and reconstructed image quality. We evaluate performance of the new projector pair for maximum a posteriori (MAP) reconstruction of simulated and in vivo data. The new projector results in only a small loss in resolution towards the edge of the field-of-view when compared to the fully 3-D geometric projector and requires an order of magnitude less computation.

In order to compare anatomical and functional brain imaging data across subjects, the images must first be registered to a common coordinate system in which anatomical features are aligned. Intensity-based volume registration methods can align subcortical structures well, but the variability in sulcal folding patterns typically results in misalignment of the cortical surface. Conversely, surface-based registration using sulcal features can produce excellent cortical alignment but the mapping between brains is restricted to the cortical surface. Here we describe a method for volumetric registration that also produces an accurate one-to-one point correspondence between cortical surfaces. This is achieved by first parameterizing and aligning the cortical surfaces using sulcal landmarks. We then use a constrained harmonic mapping to extend this surface correspondence to the entire cortical volume. Finally, this mapping is refined using an intensity-based warp. We demonstrate the utility of the method by applying it to T1-weighted magnetic resonance images (MRIs). We evaluate the performance of our proposed method relative to existing methods that use only intensity information; for this comparison we compute the intersubject alignment of expert-labeled subcortical structures after registration.

We describe a fast and globally convergent fully four-dimensional incremental gradient (4DIG) algorithm to estimate the continuous-time tracer density from list mode positron emission tomography (PET) data. Detection of 511-keV photon pairs produced by positron-electron annihilation is modeled as an inhomogeneous Poisson process whose rate function is parameterized using cubic B-splines. The rate functions are estimated by minimizing the cost function formed by the sum of the negative log-likelihood of arrival times, spatial and temporal roughness penalties, and a negativity penalty. We first derive a computable bound for the norm of the optimal temporal basis function coefficients. Based on this bound we then construct and prove convergence of an incremental gradient algorithm. Fully 4-D simulations demonstrate the substantially faster convergence behavior of the 4DIG algorithm relative to preconditioned conjugate gradient. Four-dimensional reconstructions of real data are also included to illustrate the performance of this method.

We derive computationally efficient methods for the estimation of the mean and variance properties of penalized likelihood dynamic positron emission tomography (PET) images. This allows us to predict the accuracy of reconstructed activity estimates and to compare reconstruction algorithms theoretically. We combine a bin-mode approach in which data is modeled as a collection of independent Poisson random variables at each spatiotemporal bin with the space-time separabilities in the imaging equation and penalties to derive rapidly computable analytic mean and variance approximations. We use these approximations to compare bias/variance properties of our dynamic PET image reconstruction algorithm with those of multiframe static PET reconstructions.

Randoms precorrected positron emission tomography (PET) data is formed as the difference of two Poisson random variables. Its exact probability mass function (PMF) is inconvenient for use in likelihood-based iterative image reconstruction as it contains an infinite summation. The shifted Poisson model is a tractable approximation to this PMF but requires that negative values are truncated, resulting in positively biased reconstructions in low count studies. Here we analyze the properties of the exact PMF and propose a simple but accurate approximation that allows negative valued data. We investigate the properties of this approximation and demonstrate its application to penalized maximum likelihood image reconstruction.

The Fisher information matrix (FIM) plays a key role in the analysis and applications of statistical image reconstruction methods based on Poisson data models. The elements of the FIM are a function of the reciprocal of the mean values of sinogram elements. Conventional plug-in FIM estimation methods do not work well at low counts, where the FIM estimate is highly sensitive to the reciprocal mean estimates at individual detector pairs. A generalized error look-up table (GELT) method is developed to estimate the reciprocal of the mean of the sinogram data. This approach is also extended to randoms precorrected data. Based on these techniques, an accurate FIM estimate is obtained for both Poisson and randoms precorrected data. As an application, the new GELT method is used to improve resolution uniformity and achieve near-uniform image resolution in low count situations.

We describe a method for computing a continuous time estimate of tracer density using list-mode positron emission tomography data. The rate function in each voxel is modeled as an inhomogeneous Poisson process whose rate function can be represented using a cubic B-spline basis. The rate functions are estimated by maximizing the likelihood of the arrival times of detected photon pairs over the control vertices of the spline, modified by quadratic spatial and temporal smoothness penalties and a penalty term to enforce nonnegativity. Randoms rate functions are estimated by assuming independence between the spatial and temporal randoms distributions. Similarly, scatter rate functions are estimated by assuming spatiotemporal independence and that the temporal distribution of the scatter is proportional to the temporal distribution of the trues. A quantitative evaluation was performed using simulated data and the method is also demonstrated in a human study using 11C-raclopride.

The human cerebral cortex is topologically equivalent to a sheet and can be considered topologically spherical if it is closed at the brain stem. Low-level segmentation of magnetic resonance (MR) imagery typically produces cerebral volumes whose tessellations are not topologically spherical. We present a novel algorithm that analyzes and constrains the topology of a volumetric object. Graphs are formed that represent the connectivity of voxel segments in the foreground and background of the image. These graphs are analyzed and minimal corrections to the volume are made prior to tessellation. We apply the algorithm to a simple test object and to cerebral white matter masks generated by a low-level tissue identification sequence. We tessellate the resulting objects using the marching cubes algorithm and verify their topology by computing their Euler characteristics. A key benefit of the algorithm is that it localizes the change to a volume to the specific areas of its topological defects.

We have evaluated the performance of two three-dimensional (3-D) reconstruction algorithms with data acquired from microPET, a high resolution tomograph dedicated to small animal imaging. The first was a linear filtered-backprojection algorithm (FBP) with reprojection of the missing data, and the second was a statistical maximum a posteriori probability algorithm (MAP). The two algorithms were evaluated in terms of their resolution performance, both in phantoms and in vivo. Sixty independent realizations of a phantom simulating the brain of a baby monkey were acquired, each containing three million counts. Each of these realizations was reconstructed independently with both algorithms. The ensemble of the 60 reconstructed realizations was used to estimate the standard deviation as a measure of the noise for each reconstruction algorithm. More detail was recovered in the MAP reconstruction without an increase in noise relative to FBP. Studies in a simple cylindrical compartment phantom demonstrated improved recovery of known activity ratios with MAP. Finally, in vivo studies also demonstrated a clear improvement in spatial resolution using the MAP algorithm. The quantitative accuracy of the MAP reconstruction was also evaluated by comparison with autoradiography and direct well counting of tissue samples and was shown to be superior.

We derive approximate analytical expressions for the local impulse response and covariance of images reconstructed from fully three-dimensional (3-D) positron emission tomography (PET) data using maximum a posteriori (MAP) estimation. These expressions explicitly account for the spatially variant detector response and sensitivity of a 3-D tomograph. The resulting spatially variant impulse response and covariance are computed using 3-D Fourier transforms. A truncated Gaussian distribution is used to account for the effect on the variance of the nonnegativity constraint used in MAP reconstruction. Using Monte Carlo simulations and phantom data from the microPET small animal scanner, we show that the approximations provide reasonably accurate estimates of contrast recovery and covariance of MAP reconstruction for priors with quadratic energy functions. We also describe how these analytical results can be used to achieve near-uniform contrast recovery throughout the reconstructed volume.

We examine the spatial resolution and variance properties of PET images reconstructed using maximum a posteriori (MAP) or penalized-likelihood methods. Resolution is characterized by the contrast recovery coefficient (CRC) of the local impulse response. Simplified approximate expressions are derived for the local impulse response CRC's and variances for each voxel. Using these results we propose a practical scheme for selecting spatially variant smoothing parameters to optimize lesion detectability through maximization of the local CRC-to-noise ratio in the reconstructed image.

The parameters of the prior, the hyperparameters, play an important role in Bayesian image estimation. Of particular importance for the case of Gibbs priors is the global hyperparameter, beta, which multiplies the Hamiltonian. Here we consider maximum likelihood (ML) estimation of beta from incomplete data, i.e., problems in which the image, which is drawn from a Gibbs prior, is observed indirectly through some degradation or blurring process. Important applications include image restoration and image reconstruction from projections. Exact ML estimation of beta from incomplete data is intractable for most image processing. Here we present an approximate ML estimator that is computed simultaneously with a maximum a posteriori (MAP) image estimate. The algorithm is based on a mean field approximation technique through which multidimensional Gibbs distributions are approximated by a separable function equal to a product of one-dimensional (1-D) densities. We show how this approach can be used to simplify the ML estimation problem. We also show how the Gibbs-Bogoliubov-Feynman (GBF) bound can be used to optimize the approximation for a restricted class of problems. We present the results of a Monte Carlo study that examines the bias and variance of this estimator when applied to image restoration.

We describe a new approach to imaging neural current sources from measurements of the magnetoencephalogram (MEG) associated with sensory, motor, or cognitive brain activation. Many previous approaches to this problem have concentrated on the use of weighted minimum norm (WMN) inverse methods. While these methods ensure a unique solution, they do not introduce information specific to the MEG inverse problem, often producing overly smoothed solutions and exhibiting severe sensitivity to noise. We describe a Bayesian formulation of the inverse problem in which a Gibbs prior is constructed to reflect the sparse focal nature of neural current sources associated with evoked response data. We demonstrate the method with simulated and experimental phantom data, comparing its performance with several WMN methods.

We describe a computerized method to automatically find and label the cortical surface in three-dimensional (3-D) magnetic resonance (MR) brain images. The approach we take is to model a prelabeled brain atlas as a physical object and give it elastic properties, allowing it to warp itself onto regions in a preprocessed image. Preprocessing consists of boundary-finding and a morphological procedure which automatically extracts the brain and sulci from an MR image and provides a smoothed representation of the brain surface to which the deformable model can rapidly converge. Our deformable models are energy-minimizing elastic surfaces that can accurately locate image features. The models are parameterized with 3-D bicubic B-spline surfaces. We design the energy function such that cortical fissure (sulci) points on the model are attracted to fissure points on the image and the remaining model points are attracted to the brain surface. A conjugate gradient method minimizes the energy function, allowing the model to automatically converge to the smoothed brain surface. Finally, labels are propagated from the deformed atlas onto the high-resolution brain surface.

Describes a new class of algorithms for active noise control (ANC) for use in environments in which impulsive noise is present. The well known filtered-X and filtered-U ANC algorithms are designed to minimize the variance of a measured error signal. For impulsive noise, which can be modeled using non-Gaussian stable processes, these standard approaches are not appropriate since the second order moments do not exist. The authors propose a new class of adaptive algorithms for ANC that are based on the minimization of a fractional lower order moment, p<2. By studying the effect of p on the convergence behavior of adaptive algorithms, they observe that superior performance is obtained by choosing p/spl ap//spl alpha/ where /spl alpha/<2 is a parameter reflecting the degree of impulsiveness of the noise. Applications of this approach to noise cancellation in a duct are presented.

In neuromagnetic source reconstruction, a functional map of neural activity is constructed from noninvasive magnetoencephalographic (MEG) measurements. The overall reconstruction problem is under-determined, so some form of source modeling must be applied. The authors review the two main classes of reconstruction techniques-parametric current dipole models and nonparametric distributed source reconstructions. Current dipole reconstructions use a physically plausible source model, but are limited to cases in which the neural currents are expected to be highly sparse and localized. Distributed source reconstructions can be applied to a wider variety of cases, but must incorporate an implicit source model in order to arrive at a single reconstruction. The authors examine distributed source reconstruction in a Bayesian framework to highlight the implicit nonphysical Gaussian assumptions of minimum norm based reconstruction algorithms. They conclude with a brief discussion of alternative non-Gaussian approaches.

The most widely used model for electroencephalography (EEG) and magnetoencephalography (MEG) assumes a quasi-static approximation of Maxwell's equations and a piecewise homogeneous conductor model. Both models contain an incremental field element that linearly relates an incremental source element (current dipole) to the field or voltage at a distant point. The explicit form of the field element is dependent on the head modeling assumptions and sensor configuration. Proper characterization of this incremental element is crucial to the inverse problem. The field element can be partitioned into the product of a vector dependent on sensor characteristics and a matrix kernel dependent only on head modeling assumptions. The authors present the matrix kernels for the general boundary element model (BEM) and for MEG spherical models. They show how these kernels are easily interchanged in a linear algebraic framework that includes sensor specifics such as orientation and gradiometer configuration. They then describe how this kernel is easily applied to "gain" or "transfer" matrices used in multiple dipole and source imaging models.

We describe an approximate ML estimator for the hyperparameters of a Gibbs prior which can be computed simultaneously with a maximum a posteriori (MAP) image estimate. The algorithm is based on a mean field approximation technique through which multidimensional Gibbs distributions are approximated by a separable function equal to a product of one dimensional densities. We show how this approach can be used to simplify the ML estimation problem. We also show how the Gibbs-Bogoliubov-Feynman bound can be used to optimize the approximation for a restricted class of problems.

We describe a method for automatically labelling regions of three-dimensional (3-D) Magnetic Resonance (MR) scans of human brains. Labelling consists of attaching anatomic names to particular regions of the cortical surface that appear in these images. The approach we take is to automatically match a deformable anatomical atlas model to preprocessed brain images, where preprocessing consists of 3-D Marr-Hildreth edge detection and morphological operations. These filtering operations automatically extract the brain and sulci from an MR image and provide a smoothed representation of the brain surface to which the deformable model can rapidly converge. The model itself is a 3-D B-spline surface whose control vertices are chosen to minimize a cost function that reflects the distance of the model from boundary-like features in the image. Minimization takes place using a conjugate gradient technique.>

The authors describe conjugate gradient algorithms for reconstruction of transmission and emission PET images. The reconstructions are based on a Bayesian formulation, where the data are modeled as a collection of independent Poisson random variables and the image is modeled using a Markov random field. A conjugate gradient algorithm is used to compute a maximum a posteriori (MAP) estimate of the image by maximizing over the posterior density. To ensure nonnegativity of the solution, a penalty function is used to convert the problem to one of unconstrained optimization. Preconditioners are used to enhance convergence rates. These methods generally achieve effective convergence in 15-25 iterations. Reconstructions are presented of an (18)FDG whole body scan from data collected using a Siemens/CTI ECAT931 whole body system. These results indicate significant improvements in emission image quality using the Bayesian approach, in comparison to filtered backprojection, particularly when reprojections of the MAP transmission image are used in place of the standard attenuation correction factors.

A method of computing the velocity field and pressure distribution from a sequence of ultrafast CT (UFCT) cardiac images is demonstrated. UFCT multi-slice cine imaging gives a series of tomographic slices covering the volume of the heart at a rate of 17 frames per second. The complete volume data set can be modeled using equations of continuum theory and through regularization, velocity vectors of both blood and tissue can be determined at each voxel in the volume. The authors present a technique to determine the pressure distribution throughout the volume of the left ventricle using the computed velocity field. A numerical algorithm is developed by discretizing the pressure Poisson equation (PPE), which Is based on the Navier-Stokes equation. The algorithm is evaluated using a mathematical phantom of known velocity and pressure-Couette flow. It is shown that the algorithm based on the PPE can reconstruct the pressure distribution using only the velocity data. Furthermore, the PPE is shown to be robust in the presence of noise. The velocity field and pressure distribution derived from a UFCT study of a patient are also presented.

A novel graph theoretic approach for data clustering is presented and its application to the image segmentation problem is demonstrated. The data to be clustered are represented by an undirected adjacency graph G with arc capacities assigned to reflect the similarity between the linked vertices. Clustering is achieved by removing arcs of G to form mutually exclusive subgraphs such that the largest inter-subgraph maximum flow is minimized. For graphs of moderate size ( approximately 2000 vertices), the optimal solution is obtained through partitioning a flow and cut equivalent tree of G, which can be efficiently constructed using the Gomory-Hu algorithm (1961). However for larger graphs this approach is impractical. New theorems for subgraph condensation are derived and are then used to develop a fast algorithm which hierarchically constructs and partitions a partially equivalent tree of much reduced size. This algorithm results in an optimal solution equivalent to that obtained by partitioning the complete equivalent tree and is able to handle very large graphs with several hundred thousand vertices. The new clustering algorithm is applied to the image segmentation problem. The segmentation is achieved by effectively searching for closed contours of edge elements (equivalent to minimum cuts in G), which consist mostly of strong edges, while rejecting contours containing isolated strong edges. This method is able to accurately locate region boundaries and at the same time guarantees the formation of closed edge contours.>

A highly accurate approximation is proposed for computing the joint likelihood for first-order Gauss-Markov random fields (GMRFs) defined on irregularly shaped lattices. The problem in computing the likelihood lies in evaluating the determinant of a very large matrix B. Its exact evaluation is limited to either very small irregular regions or a few regularly shaped regions. The approximation proposed here is based on an eigenanalysis of B.

A graph-theoretic approach for image segmentation is presented. The pixels of the image are represented by the vertices of an undirected adjacency graph G. All neighboring pairs of pixels are linked by arcs with capacities assigned to reflect the strength of an edge element between the linked vertices. Segmentation is achieved by removing arcs corresponding to selected minimum cuts of G to form mutually exclusive subgraphs such that the largest intersubgraph maximum flow is minimized. This is equivalent to partitioning the image using closed contours of edge elements, which consist mostly of strong edges. The method accurately locates region boundaries and at the same time rejects contours containing isolated strong edges. The minimum cuts in G can be computed from a partially cut-equivalent tree of G. A fast algorithm for constructing partially equivalent trees that can handle graphs with several hundred thousand vertices is developed.>

A widely used approach to restoration of images degraded by Poisson noise has been to represent the degraded image as a linear transformation of the intact image plus zero-mean, uncorrelated additive noise. The explicit Poisson model is retained and the authors formulate maximum a posteriori (MAP) image restoration using homogeneous, isotropic Markov random field priors. For this purpose, a deterministic algorithm based upon the expectation maximization (EM) method of optimization is derived. The performance of candidate priors is compared.>

An unsupervised hierarchical segmentation method is described, and its application to tissue classification in magnetic resonance (MR) images of the human brain is demonstrated. The images are modeled as a mosaic of homogeneous subimages where each subimage is modeled as a first-order Gauss-Markov random field (GMRF) with unknown parameters. The segmentation goal is to group the pixels into regions which, under a suitable hypothesis, are homogeneous GMRFs. The image is represented by a quadtree, and its segmentation is achieved by splitting and merging the image, followed by a step-wise maximum likelihood agglomerative clustering procedure. The difficulty of evaluating the likelihood for irregularly shaped regions is overcome using a highly accurate approximation for the determinant of the covariance matrix based on eigenanalysis.>

A nonlinear simplex algorithm is developed for solving a special class of linearly constrained optimization problems. The cost for this class of problem is a sequence of functionals indexed by the positive integers P such that for P=1 the problem reduces to a standard linear programming problem and for P to infinity the cost function selects the solution with the minimum number of nonzero elements in the solution vector. Two theorems, analogous to the fundamental theorem of linear programming, are given as the basis for the nonlinear simplex algorithm. Applications of this technique are discussed for system design and signal processing problems in which an optimally sparse solution vector is desired. Examples of such problems include FIR filter design, placement of beamformer arrays and seismic deconvolution.>

A general method for the optimization of convex or concave cost functions over the intersection of convex constraint sets is described for applications in signal reconstruction and restorarion. A unique continuous function is obtained, under certain conditions, by employing Fenchel's duality theorem to give a finite dimensional dual problem. This approach allows the solution of very complex constrained problems by separating the constraints into those related to the solution and those related to additive noise statistics. The method is applied to computed tomography with noisy data of which the noise covariance and bounds on the solution are known approximately. A fast implementation of the required optimization procedure is given and the resulting solution is shown to be significantly better than a suboptimal feasible solution.

The autoregressive and Prony methods for spectral estimation do not make full use of the statistics of the additive noise. These statistics are included in the following as constraints on the solution. The finite number of data samples imposes a limit on the resolution thus making it possible to approximate random signals by a set of deterministic signals which can be modelled in terms of a finite set of time-invariant parameters, θ. The non-Toeplitz covariance calculated from the data is modelled here using the θ parameters. Statistical constraints and covariance modelling are combined to produce non-linear methods for spectral estimation. The resulting higher resolution requires high computational complexity; this can often be substantially reduced by using knowledge-based techniques.

A unified theoretical treatment of constrained optimisation methods for tomographic and spectral estimation from discrete data is given. The solution is shown to be equivalent to the unconstrained optimisation of a dual functional in which the image or spectrum is modelled in terms of a Lagrange multiplier vector and the kernel of the constraint integrals. In order to obtain the best possible solution it is important to consider the effects of noise in the constraints. The problem is reformulated using the above models and the exact data is replaced with the noise statistics as constraints; this is solved using a penalty method. A very fast direct algorithm is also introduced which matches the noise variance provided the signal to noise ratio is approximately known.

Polar periodicities associated with the Radon transform and Lagrange multipliers are studies to obtain the sample/ recovery relationships. By considering a window function, which is zero outside the object boundaries and inside its transparent areas, a Fourier domain algorithm is introduced in which the Lagrange multipliers are determined using the window as a constraint. By expressing the backprojection image using the basis function of the projection transform a reconstruction model for fan beam projections is given. A very fast, iterative, space domain algorithm for divergent ray projections is introduced which gives the minimum energy reconstruction by exploiting a special matrix structure; this method matches the noise statistics and can utilise 'a priori' knowledge.