Abstract. Many inverse problems and signal processing problems involve low-rank regularizers based on the nuclear norm. Commonly, proximal gradient methods (PGMs) are adopted to solve this type of nonsmooth problems, as they can offer fast and guaranteed convergence. However, PGM methods cannot be simply applied in settings where low-rank models are imposed locally on overlapping patches; therefore, heuristic approaches have been proposed that lack convergence guarantees. In this work we propose to replace the nuclear norm with a smooth approximation in which a Huber-type function is applied to each singular value. By providing a theoretical framework based on singular value function theory, we show that important properties can be established for the proposed regularizer, such as convexity, differentiability, and Lipschitz continuity of the gradient. Moreover, we provide a closed-form expression for the regularizer gradient, enabling the use of standard iterative gradient-based optimization algorithms (e.g., the nonlinear conjugate gradient) that can easily address the case of overlapping patches and have well-known convergence guarantees. In addition, we provide a novel step-size selection strategy based on a quadratic majorizer of the line search function that leverages the Huber characteristics of the proposed regularizer. Finally, we assess the proposed optimization framework by providing empirical results in dynamic magnetic resonance imaging (MRI) reconstruction in the context of local low-rank models with overlapping patches.

Deep learning (DL) methods can reconstruct highly accelerated magnetic resonance imaging (MRI) scans, but they rely on application-specific large training datasets and often generalize poorly to out-of-distribution data. Self-supervised deep learning algorithms perform scan-specific reconstructions, but still require complicated hyperparameter tuning based on the acquisition and often offer limited acceleration. This work develops a bilevel-optimized implicit neural representation (INR) approach for scan-specific MRI reconstruction. The method explicitly formulates the undersampled MRI reconstruction problem as a bilevel optimization problem and automatically optimizes the multidimensional hyperparameters of the reconstruction method for a given acquisition protocol, enabling a tailored reconstruction without training data. The proposed algorithm uses Gaussian process regression to optimize INR hyperparameters, accommodating various acquisitions. The INR includes a trainable positional encoder for high-dimensional feature embedding and a small multilayer perceptron for decoding. The bilevel optimization is computationally efficient, requiring only a few minutes per typical 2D Cartesian scan. On the scanner hardware, the subsequent scan-specific reconstruction-using offline-optimized hyperparameters-is completed within seconds, while achieving comparable or improved image quality compared to previous model-based and self-supervised learning methods.

Data assimilation (DA) integrates observations with a dynamical model to estimate states of PDE-governed systems. Model-driven methods (e.g., Kalman Filter, Particle Filter) presuppose full knowledge of the true dynamics, which is not always satisfied in practice, while purely data-driven solvers learn a deterministic mapping between observations and states and therefore miss the intrinsic stochasticity of real processes. Recently, score-based diffusion models have shown promise for DA by learning a global diffusion prior to represent stochastic dynamics. However, their one-shot generation lacks stepwise physical consistency and struggles with complex stochastic processes. To address these issues, we propose FlowDAS, a generative DA framework that employs stochastic interpolants to learn state transition dynamics through step-by-step stochastic updates. By incorporating observations into each transition, FlowDAS can produce stable, measurement-consistent forecasts. Experiments on Lorenz-63, Navier–Stokes super-resolution/sparse-observation scenarios, and large-scale weather forecasting—where dynamics are partly or wholly unknown—show that FlowDAS surpasses model-driven methods, neural operators, and score-based baselines in accuracy and physical plausibility. Our implementation is available at https://github.com/umjiayx/FlowDAS.

Diffusion models can learn strong image priors from underlying data distribution and use them to solve inverse problems, but the training process is computationally expensive and requires lots of data. Such bottlenecks prevent most existing works from being feasible for high-dimensional and high-resolution data such as 3D images. This paper proposes a method to learn an efficient data prior for the entire image by training diffusion models only on patches of images. Specifically, we propose a patch-based position-aware diffusion inverse solver, called PaDIS, where we obtain the score function of the whole image through scores of patches and their positional encoding and utilize this as the prior for solving inverse problems. First of all, we show that this diffusion model achieves an improved memory efficiency and data efficiency while still maintaining the capability to generate entire images via positional encoding. Additionally, the proposed PaDIS model is highly flexible and can be plugged in with different diffusion inverse solvers (DIS). We demonstrate that the proposed PaDIS approach enables solving various inverse problems in both natural and medical image domains, including CT reconstruction, deblurring, and superresolution, given only patch-based priors. Notably, PaDIS outperforms previous DIS methods trained on entire image priors in the case of limited training data, demonstrating the data efficiency of our proposed approach by learning patch-based prior.

Diffusion models face significant challenges when employed for large-scale medical image reconstruction in real practice such as 3D Computed Tomography (CT). Due to the demanding memory, time, and data requirements, it is difficult to train a diffusion model directly on the entire volume of high-dimensional data to obtain an efficient 3D diffusion prior. Existing works utilizing diffusion priors on single 2D image slice with hand-crafted cross-slice regularization would sacrifice the z-axis consistency, which results in severe artifacts along the z-axis. In this work, we propose a novel framework that enables learning the 3D image prior through position-aware 3D-patch diffusion score blending for reconstructing large-scale 3D medical images. To the best of our knowledge, we are the first to utilize a 3D-patch diffusion prior for 3D medical image reconstruction. Extensive experiments on sparse view and limited angle CT reconstruction show that our DiffusionBlend method significantly outperforms previous methods and achieves state-of-the-art performance on real-world CT reconstruction problems with high-dimensional 3D image (i.e., $256 \times 256 \times 500$). Our algorithm also comes with better or comparable computational efficiency than previous state-of-the-art methods. Code is available at https://github.com/efzero/DiffusionBlend.

Oscillating Steady-State Imaging (OSSI) is a recently developed fMRI acquisition method that can provide 2 to 3 times higher SNR than standard fMRI approaches. However, because the OSSI signal exhibits a nonlinear oscillation pattern, one must acquire and combine nc(e.g., 10) OSSI images to get an image that is free of oscillation for fMRI, and fully sampled acquisitions would compromise temporal resolution. To improve temporal resolution and accurately model the nonlinearity of OSSI signals, instead of using subspace models that are not well suited for the data, we build the MR physics for OSSI signal generation as a regularizer for the undersampled reconstruction. Our proposed physics-based manifold model turns the disadvantages of OSSI acquisition into advantages and enables joint reconstruction and quantification. OSSI manifold model (OSSIMM) outperforms subspace models and reconstructs high-resolution fMRI images with a factor of 12 acceleration and without spatial or temporal smoothing. Furthermore, OSSIMM can dynamically quantify important physics parameters, including R*2maps, with a temporal resolution of 150 ms.

Mixtures of probabilistic principal component analysis (MPPCA) is a well-known mixture model extension of principal component analysis (PCA). Similar to PCA, MPPCA assumes the data samples in each mixture contain homoscedastic noise. However, datasets with heterogeneous noise across samples are becoming increasingly common, as larger datasets are generated by collecting samples from several sources with varying noise profiles. The performance of MPPCA is suboptimal for data with heteroscedastic noise across samples. This paper proposes a heteroscedastic mixtures of probabilistic PCA technique (HeMPPCAT) that uses a gen-eralized expectation-maximization (GEM) algorithm to jointly estimate the unknown underlying factors, means, and noise variances under a heteroscedastic noise setting. Simulation results illustrate the improved factor estimates and clustering accuracies of HeMPPCAT compared to MPPCA.

Iterative neural networks (INN) are rapidly gaining attention for solving inverse problems in imaging, image processing, and computer vision. INNs combine regression NNs and an iterative model-based image reconstruction (MBIR) algorithm, often leading to both good generalization capability and outperforming reconstruction quality over existing MBIR optimization models. This paper proposes the first fast and convergent INN architecture, Momentum-Net, by generalizing a block-wise MBIR algorithm that uses momentum and majorizers with regression NNs. For fast MBIR, Momentum-Net uses momentum terms in extrapolation modules, and noniterative MBIR modules at each iteration by using majorizers, where each iteration of Momentum-Net consists of three core modules: image refining, extrapolation, and MBIR. Momentum-Net guarantees convergence to a fixed-point for general differentiable (non)convex MBIR functions (or data-fit terms) and convex feasible sets, under two asymptomatic conditions. To consider data-fit variations across training and testing samples, we also propose a regularization parameter selection scheme based on the "spectral spread" of majorization matrices. Numerical experiments for light-field photography using a focal stack and sparse-view computational tomography demonstrate that, given identical regression NN architectures, Momentum-Net significantly improves MBIR speed and accuracy over several existing INNs; it significantly improves reconstruction quality compared to a state-of-the-art MBIR method in each application.

Accurate scatter estimation is important in quantitative SPECT for improving image contrast and accuracy. With a large number of photon histories, Monte-Carlo (MC) simulation can yield accurate scatter estimation, but is computationally expensive. Recent deep learning-based approaches can yield accurate scatter estimates quickly, yet full MC simulation is still required to generate scatter estimates as ground truth labels for all training data. Here we propose a physics-guided weakly supervised training framework for fast and accurate scatter estimation in quantitative SPECT by using a 100× shorter MC simulation as weak labels and enhancing them with deep neural networks. Our weakly supervised approach also allows quick fine-tuning of the trained network to any new test data for further improved performance with an additional short MC simulation (weak label) for patient-specific scatter modelling. Our method was trained with 18 XCAT phantoms with diverse anatomies / activities and then was evaluated on 6 XCAT phantoms, 4 realistic virtual patient phantoms, 1 torso phantom and 3 clinical scans from 2 patients for 177Lu SPECT with single / dual photopeaks (113, 208 keV). Our proposed weakly supervised method yielded comparable performance to the supervised counterpart in phantom experiments, but with significantly reduced computation in labeling. Our proposed method with patient-specific fine-tuning achieved more accurate scatter estimates than the supervised method in clinical scans. Our method with physics-guided weak supervision enables accurate deep scatter estimation in quantitative SPECT, while requiring much lower computation in labeling, enabling patient-specific fine-tuning capability in testing.

Optimizing k-space sampling trajectories is a promising yet challenging topic for fast magnetic resonance imaging (MRI). This work proposes to optimize a reconstruction method and sampling trajectories jointly concerning image reconstruction quality in a supervised learning manner. We parameterize trajectories with quadratic B-spline kernels to reduce the number of parameters and apply multi-scale optimization, which may help to avoid sub-optimal local minima. The algorithm includes an efficient non-Cartesian unrolled neural network-based reconstruction and an accurate approximation for backpropagation through the non-uniform fast Fourier transform (NUFFT) operator to accurately reconstruct and back-propagate multi-coil non-Cartesian data. Penalties on slew rate and gradient amplitude enforce hardware constraints. Sampling and reconstruction are trained jointly using large public datasets. To correct for possible eddy-current effects introduced by the curved trajectory, we use a pencil-beam trajectory mapping technique. In both simulations and in- vivo experiments, the learned trajectory demonstrates significantly improved image quality compared to previous model-based and learning-based trajectory optimization methods for 10× acceleration factors. Though trained with neural network-based reconstruction, the proposed trajectory also leads to improved image quality with compressed sensing-based reconstruction.

The recent trend in regularization methods for inverse problems is to replace handcrafted sparsifying operators with data-driven approaches. Although using such machine learning techniques often improves image reconstruction methods, the results can depend significantly on the learning methodology. This paper compares two supervised learning methods. First, the paper considers a transform learning approach and, to learn the transform, introduces a variant on the Procrustes method for wide matrices with orthogonal rows. Second, we consider a bilevel convolutional filter learning approach. Numerical experiments show the learned transform performs worse for denoising than both the handcrafted finite difference transform and the learned filters, which perform similarly. Our results motivate the use of bilevel learning.

This paper discusses algorithms for phase retrieval where the measurements follow independent Poisson distributions. We developed an optimization problem based on maximum likelihood estimation (MLE) for the Poisson model and applied Wirtinger flow algorithm to solve it. Simulation results with a random Gaussian sensing matrix and Poisson measurement noise demonstrated that the Wirtinger flow algorithm based on the Poisson model produced higher quality reconstructions than when algorithms derived from Gaussian noise models (Wirtinger flow, Gerchberg Saxton) are applied to such data, with significantly improved computational efficiency.

The Fast Fourier Transform (FFT) is a fundamental algorithm in signal processing; significant efforts have been made to improve its performance using software optimizations and specialized hardware accelerators. Computational imaging modalities, such as MRI, often rely on the Non-uniform Fast Fourier Transform (NuFFT), a variant of the FFT for processing data acquired from non-uniform sampling patterns. The most time-consuming step of the NuFFT algorithm is “gridding;” wherein non-uniform samples are interpolated to allow a uniform FFT to be computed over the data. Each non-uniform sample affects a window of non-contiguous memory locations, resulting in poor cache and memory bandwidth utilization. As a result, gridding can account for more than 99.6% of the NuFFT computation time, while the FFT requires less than 0.4%. We present Slice-and-Dice, a novel approach to the NuFFT's gridding step that eliminates the presorting operations required by prior methods and maps more efficiently to hardware. Our GPU implementation achieves gridding speedups of over 250× and 16× vs prior state-of-the-art CPU and GPU implementations, respectively. We achieve further speedup and energy efficiency gains by implementing Slice-and-Dice in hardware with JIGSAW, a streaming hardware accelerator for non-uniform data gridding. JIGSAW uses stall-free fixed-point pipelines to process M non-uniform samples in approximately M cycles, irrespective of sampling pattern-yielding speedups of over 1500× the CPU baseline and 36× the state-of-the-art GPU implementation, consuming $\sim 200\mathrm{m}\mathrm{W}$ power and $\sim 12\mathrm{m}\mathrm{m}^{2}$ area in 16 nm technology. Slice-and-Dice GPU and JIGSAW ASIC implementations achieve unprecedented end-to-end NuFFT speedups of 8× and 36× compared to the state-of-the-art GPU implementation, respectively.

Despite continuing research into inter-GPU communication mechanisms, extracting performance from multi-GPU systems remains a significant challenge. Inter-GPU communication via bulk DMA-based transfers exposes data transfer latency on the GPU’s critical execution path because these large transfers are logically interleaved between compute kernels. Conversely, fine-grained peer-to-peer memory accesses during kernel execution lead to memory stalls that can exceed the GPUs’ ability to cover these operations via multi-threading. Worse yet, these sub-cacheline transfers are highly inefficient on current inter-GPU interconnects. To remedy these issues, we propose PROACT, a system enabling remote memory transfers with the programmability and pipeline advantages of peer-to-peer stores, while achieving interconnect efficiency that rivals bulk DMA transfers. Combining compile-time instrumentation with fine-grain tracking of data block readiness within each GPU, PROACT enables interconnect-friendly data transfers while hiding the transfer latency via pipelining during kernel execution. This work describes both hardware and software implementations of PROACT and demonstrates the effectiveness of a PROACT software prototype on three generations of GPU hardware and interconnects. Achieving near-ideal interconnect efficiency, PROACT realizes a mean speedup of 3.0× over single-GPU performance for 4-GPU systems, capturing 83% of available performance opportunity. On a 16-GPU NVIDIA DGX-2 system, we demonstrate an 11.0× average strong-scaling speedup over single-GPU performance, 5.3× better than a bulk DMA-based approach.

Digital breast tomosynthesis (DBT) is a quasi-three-dimensional imaging modality that can reduce false negatives and false positives in mass lesion detection caused by overlapping breast tissue in conventional two-dimensional (2D) mammography. The patient dose of a DBT scan is similar to that of a single 2D mammogram, while acquisition of each projection view adds detector readout noise. The noise is propagated to the reconstructed DBT volume, possibly obscuring subtle signs of breast cancer such as microcalcifications (MCs). This study developed a deep convolutional neural network (DCNN) framework for denoising DBT images with a focus on improving the conspicuity of MCs as well as preserving the ill-defined margins of spiculated masses and normal tissue textures. We trained the DCNN using a weighted combination of mean squared error (MSE) loss and adversarial loss. We configured a dedicated x-ray imaging simulator in combination with digital breast phantoms to generate realistic in silico DBT data for training. We compared the DCNN training between using digital phantoms and using real physical phantoms. The proposed denoising method improved the contrast-to-noise ratio (CNR) and detectability index (d') of the simulated MCs in the validation phantom DBTs. These performance measures improved with increasing training target dose and training sample size. Promising denoising results were observed on the transferability of the digital-phantom-trained denoiser to DBT reconstructed with different techniques and on a small independent test set of human subject DBT images.

This paper examines a combined supervised-unsupervised framework involving dictionary-based blind learning and deep supervised learning for MR image reconstruction from under-sampled k-space data. A major focus of the work is to investigate the possible synergy of learned features in traditional shallow reconstruction using adaptive sparsity-based priors and deep prior-based reconstruction. Specifically, we propose a framework that uses an unrolled network to refine a blind dictionary learning-based reconstruction. We compare the proposed method with strictly supervised deep learning-based reconstruction approaches on several datasets of varying sizes and anatomies. We also compare the proposed method to alternative approaches for combining dictionary-based methods with supervised learning in MR image reconstruction. The improvements yielded by the proposed framework suggest that the blind dictionary-based approach preserves fine image details that the supervised approach can iteratively refine, suggesting that the features learned using the two methods are complementary.

This paper proposes a new method for joint design of radiofrequency (RF) and gradient waveforms in Magnetic Resonance Imaging (MRI), and applies it to the design of 3D spatially tailored saturation and inversion pulses. The joint design of both waveforms is characterized by the ODE Bloch equations, to which there is no known direct solution. Existing approaches therefore typically rely on simplified problem formulations based on, e.g., the small-tip approximation or constraining the gradient waveforms to particular shapes, and often apply only to specific objective functions for a narrow set of design goals (e.g., ignoring hardware constraints). This paper develops and exploits an auto-differentiable Bloch simulator to directly compute Jacobians of the (Bloch-simulated) excitation pattern with respect to RF and gradient waveforms. This approach is compatible with arbitrary sub-differentiable loss functions, and optimizes the RF and gradients directly without restricting the waveform shapes. For computational efficiency, we derive and implement explicit Bloch simulator Jacobians (approximately halving computation time and memory usage). To enforce hardware limits (peak RF, gradient, and slew rate), we use a change of variables that makes the 3D pulse design problem effectively unconstrained; we then optimize the resulting problem directly using the proposed auto-differentiation framework. We demonstrate our approach with two kinds of 3D excitation pulses that cannot be easily designed with conventional approaches: Outer-volume saturation (90° flip angle), and inner-volume inversion.

Light-field (LF) reconstruction from focal stack images has diverse applications including face recognition, autonomous driving, and 3D reconstruction in virtual reality. It is a large-scale ill-conditioned inverse problem and typically requires regularized iterative algorithms to solve, which can be slow. This paper proposes a non-iterative LF reconstruction and depth estimation method based on three sequential convolutional neural networks (CNNs). The first CNN estimates an all-in-focus image from focal stack images. The second CNN estimates 4D ray depth from the estimated all-in-focus image via the first CNN, and focal stack images. The third CNN refines a Lambertian LF that is rendered using the all-in-focus image and ray depth estimated by the first and second CNNs, respectively. Numerical experiments show that the proposed CNN-based method achieves significantly more accurate and/or faster LF reconstruction, compared to a state-of-the-art sequential CNN using a single image, conventional model-based image reconstruction from a focal stack, and direct regression CNN from a focal stack.

The field of medical image reconstruction has seen roughly four types of methods. The first type tended to be analytical methods, such as filtered backprojection (FBP) for X-ray computed tomography (CT) and the inverse Fourier transform for magnetic resonance imaging (MRI), based on simple mathematical models for the imaging systems. These methods are typically fast, but have suboptimal properties such as poor resolution-noise tradeoff for CT. A second type is iterative reconstruction methods based on more complete models for the imaging system physics and, where appropriate, models for the sensor statistics. These iterative methods improved image quality by reducing noise and artifacts. The U.S. Food and Drug Administration (FDA)-approved methods among these have been based on relatively simple regularization models. A third type of methods has been designed to accommodate modified data acquisition methods, such as reduced sampling in MRI and CT to reduce scan time or radiation dose. These methods typically involve mathematical image models involving assumptions such as sparsity or low rank. A fourth type of methods replaces mathematically designed models of signals and systems with data-driven or adaptive models inspired by the field of machine learning. This article focuses on the two most recent trends in medical image reconstruction: methods based on sparsity or low-rank models and data-driven methods based on machine learning techniques.

Convolutional operator learning is gaining attention in many signal processing and computer vision applications. Learning kernels has mostly relied on so-called patch-domain approaches that extract and store many overlapping patches across training signals. Due to memory demands, patch-domain methods have limitations when learning kernels from large datasets - particularly with multi-layered structures, e.g., convolutional neural networks - or when applying the learned kernels to high-dimensional signal recovery problems. The so-called convolution approach does not store many overlapping patches, and thus overcomes the memory problems particularly with careful algorithmic designs; it has been studied within the "synthesis" signal model, e.g., convolutional dictionary learning. This paper proposes a new convolutional analysis operator learning (CAOL) framework that learns an analysis sparsifying regularizer with the convolution perspective, and develops a new convergent Block Proximal Extrapolated Gradient method using a Majorizer (BPEG-M) to solve the corresponding block multi-nonconvex problems. To learn diverse filters within the CAOL framework, this paper introduces an orthogonality constraint that enforces a tight-frame filter condition, and a regularizer that promotes diversity between filters. Numerical experiments show that, with sharp majorizers, BPEG-M significantly accelerates the CAOL convergence rate compared to the state-of-the-art block proximal gradient (BPG) method. Numerical experiments for sparse-view computational tomography show that a convolutional sparsifying regularizer learned via CAOL significantly improves reconstruction quality compared to a conventional edge-preserving regularizer. Using more and wider kernels in a learned regularizer better preserves edges in reconstructed images.

CT images are often affected by beam-hardening artifacts due to the polychromatic nature of the X-ray spectra. These artifacts appear in the image as cupping in homogeneous areas and as dark bands between dense regions such as bones. This paper proposes a simplified statistical reconstruction method for X-ray CT based on Poisson statistics that accounts for the non-linearities caused by beam hardening. The main advantages of the proposed method over previous algorithms are that it avoids the preliminary segmentation step, which can be tricky, especially for low-dose scans, and it does not require knowledge of the whole source spectrum, which is often unknown. Each voxel attenuation is modeled as a mixture of bone and soft tissue by defining density-dependent tissue fractions and maintaining one unknown per voxel. We approximate the energy-dependent attenuation corresponding to different combinations of bone and soft tissues, the so-called beam-hardening function, with the 1D function corresponding to water plus two parameters that can be tuned empirically. Results on both simulated data with Poisson sinogram noise and two rodent studies acquired with the ARGUS/CT system showed a beam hardening reduction (both cupping and dark bands) similar to analytical reconstruction followed by post-processing techniques but with reduced noise and streaks in cases with a low number of projections, as expected for statistical image reconstruction.

Quantitative yttrium-90 (Y-90) SPECT imaging is challenging due to the nature of Y-90, an almost pure beta emitter that is associated with a continuous spectrum of bremsstrahlung photons that have a relatively low yield. This paper proposes joint spectral reconstruction (JSR), a novel bremsstrahlung SPECT reconstruction method that uses multiple narrow acquisition windows with accurate multi-band forward modeling to cover a wide range of the energy spectrum. Theoretical analyses using Fisher information and Monte-Carlo (MC) simulation with a digital phantom show that the proposed JSR model with multiple acquisition windows has better performance in terms of covariance (precision) than previous methods using multi-band forward modeling with a single acquisition window, or using a single-band forward modeling with a single acquisition window. We also propose an energy-window subset (ES) algorithm for JSR to achieve fast empirical convergence and maximum-likelihood based initialization for all reconstruction methods to improve quantification accuracy in early iterations. For both MC simulation with a digital phantom and experimental study with a physical multi-sphere phantom, our proposed JSR-ES, a fast algorithm for JSR with ES, yielded higher recovery coefficients (RCs) on hot spheres over all iterations and sphere sizes than all the other evaluated methods, due to fast empirical convergence. In experimental study, for the smallest hot sphere (diameter 1.6cm), at the 20th iteration the increase in RCs with JSR-ES was 66 and 31% compared with single wide and narrow band forward models, respectively. JSR-ES also yielded lower residual count error (RCE) on a cold sphere over all iterations than other methods for MC simulation with known scatter, but led to greater RCE compared with single narrow band forward model at higher iterations for experimental study when using estimated scatter.

The goals of fMRI acquisition include high spatial and temporal resolutions with a high signal to noise ratio (SNR). Oscillating Steady-State Imaging (OSSI) is a new fMRI acquisition method that provides large oscillating signals with the potential for high SNR, but does so at the expense of spatial and temporal resolutions. The unique oscillation pattern of OSSI images makes it well suited for high-dimensional modeling. We propose a patch-tensor low-rank model to exploit the local spatial-temporal low-rankness of OSSI images. We also develop a practical sparse sampling scheme with improved sampling incoherence for OSSI. With an alternating direction method of multipliers (ADMM) based algorithm, we improve OSSI spatial and temporal resolutions with a factor of 12 acquisition acceleration and 1.3 mm isotropic spatial resolution in prospectively undersampled experiments. The proposed model yields high temporal SNR with more activation than other low-rank methods. Compared to the standard grad- ient echo (GRE) imaging with the same spatial-temporal resolution, 3D OSSI tensor model reconstruction demonstrates 2 times higher temporal SNR with 2 times more functional activation.

Dual-energy computed tomography (DECT) imaging plays an important role in advanced imaging applications due to its material decomposition capability. Image-domain decomposition operates directly on CT images using linear matrix inversion, but the decomposed material images can be severely degraded by noise and artifacts. This paper proposes a new method dubbed DECT-MULTRA for image-domain DECT material decomposition that combines conventional penalized weighted-least squares (PWLS) estimation with regularization based on a mixed union of learned transforms (MULTRA) model. Our proposed approach pre-learns a union of common-material sparsifying transforms from patches extracted from all the basis materials, and a union of cross-material sparsifying transforms from multi-material patches. The common-material transforms capture the common properties among different material images, while the cross-material transforms capture the cross-dependencies. The proposed PWLS formulation is optimized efficiently by alternating between an image update step and a sparse coding and clustering step, with both of these steps having closed-form solutions. The effectiveness of our method is validated with both XCAT phantom and clinical head data. The results demonstrate that our proposed method provides superior material image quality and decomposition accuracy compared to other competing methods.

Image reconstruction in low-count PET is particularly challenging because gammas from natural radioactivity in Lu-based crystals cause high random fractions that lower the measurement signal-to-noise-ratio (SNR). In model-based image reconstruction (MBIR), using more iterations of an unregularized method may increase the noise, so incorporating regularization into the image reconstruction is desirable to control the noise. New regularization methods based on learned convolutional operators are emerging in MBIR. We modify the architecture of an iterative neural network, BCD-Net, for PET MBIR, and demonstrate the efficacy of the trained BCD-Net using XCAT phantom data that simulates the low true coincidence count-rates with high random fractions typical for Y-90 PET patient imaging after Y-90 microsphere radioembolization. Numerical results show that the proposed BCD-Net significantly improves CNR and RMSE of the reconstructed images compared to MBIR methods using non-trained regularizers, total variation (TV) and non-local means (NLM). Moreover, BCD-Net successfully generalizes to test data that differs from the training data. Improvements were also demonstrated for the clinically relevant phantom measurement data where we used training and testing datasets having very different activity distributions and count-levels.

Low-dose CT image reconstruction has been a popular research topic in recent years. A typical reconstruction method based on post-log measurements is called penalized weighted-least squares (PWLS). Due to the underlying limitations of the post-log statistical model, the PWLS reconstruction quality is often degraded in low-dose scans. This paper investigates a shifted-Poisson (SP) model based likelihood function that uses the pre-log raw measurements that better represents the measurement statistics, together with a data-driven regularizer exploiting a Union of Learned TRAnsforms (SPULTRA). Both the SP induced data-fidelity term and the regularizer in the proposed framework are nonconvex. The proposed SPULTRA algorithm uses quadratic surrogate functions for the SP induced data-fidelity term. Each iteration involves a quadratic subproblem for updating the image, and a sparse coding and clustering subproblem that has a closed-form solution. The SPULTRA algorithm has a similar computational cost per iteration as its recent counterpart PWLS-ULTRA that uses post-log measurements, and it provides better image reconstruction quality than PWLS-ULTRA, especially in low-dose scans.

Convolutional analysis operator learning (CAOL) enables the unsupervised training of (hierarchical) convolutional sparsifying operators or autoencoders from large datasets. One can use many training images for CAOL, but a precise understanding of the impact of doing so has remained an open question. This letter presents a series of results that lend insight into the impact of dataset size on the filter update in CAOL. The first result is a general deterministic bound on errors in the estimated filters, and is followed by a bound on the expected errors as the number of training samples increases. The second result provides a high probability analogue. The bounds depend on properties of the training data, and we investigate their empirical values with real data. Taken together, these results provide evidence for the potential benefit of using more training data in CAOL.

Convolutional dictionary learning (CDL or sparsifying CDL) has many applications in image processing and computer vision. There has been growing interest in developing efficient algorithms for CDL, mostly relying on the augmented Lagrangian (AL) method or the variant alternating direction method of multipliers (ADMM). When their parameters are properly tuned, AL methods have shown fast convergence in CDL. However, the parameter tuning process is not trivial due to its data dependence and, in practice, the convergence of AL methods depends on the AL parameters for nonconvex CDL problems. To moderate these problems, this paper proposes a new practically feasible and convergent Block Proximal Gradient method using a Majorizer (BPG-M) for CDL. The BPG-M-based CDL is investigated with different block updating schemes and majorization matrix designs, and further accelerated by incorporating some momentum coefficient formulas and restarting techniques. All of the methods investigated incorporate a boundary artifacts removal (or, more generally, sampling) operator in the learning model. Numerical experiments show that, without needing any parameter tuning process, the proposed BPG-M approach converges more stably to desirable solutions of lower objective values than the existing state-of-the-art ADMM algorithm and its memory-efficient variant do. Compared with the ADMM approaches, the BPG-M method using a multi-block updating scheme is particularly useful in single-threaded CDL algorithm handling large data sets, due to its lower memory requirement and no polynomial computational complexity. Image denoising experiments show that, for relatively strong additive white Gaussian noise, the filters learned by BPG-M-based CDL outperform those trained by the ADMM approach.

This paper introduces a fast, general method for dictionary-free parameter estimation in quantitative magnetic resonance imaging (QMRI) parameter estimation via regression with kernels (PERK). PERK first uses prior distributions and the nonlinear MR signal model to simulate many parameter-measurement pairs. Inspired by machine learning, PERK then takes these parameter-measurement pairs as labeled training points and learns from them a nonlinear regression function using kernel functions and convex optimization. PERK admits a simple implementation as per-voxel nonlinear lifting of MRI measurements followed by linear minimum mean-squared error regression. We demonstrate PERK for $ {\textit {T}_{1}}, {\textit {T}_{2}}$ estimation, a well-studied application where it is simple to compare PERK estimates against dictionary-based grid search estimates and iterative optimization estimates. Numerical simulations as well as single-slice phantom and in vivo experiments demonstrate that PERK and other tested methods produce comparable $ {\textit {T}_{1}}, {\textit {T}_{2}}$ estimates in white and gray matter, but PERK is consistently at least $140\times $ faster. This acceleration factor may increase by several orders of magnitude for full-volume QMRI estimation problems involving more latent parameters per voxel.

Penalized least-squares iterative image reconstruction algorithms used for spatial resolution-limited imaging, such as functional magnetic resonance imaging (fMRI), commonly use a quadratic roughness penalty to regularize the reconstructed images. When used for complex-valued images, the conventional roughness penalty regularizes the real and imaginary parts equally. However, these imaging methods sometimes benefit from separate penalties for each part. The spatial smoothness from the roughness penalty on the reconstructed image is dictated by the regularization parameter(s). One method to set the parameter to a desired smoothness level is to evaluate the full width at half maximum of the reconstruction method's local impulse response. Previous work has shown that when using the conventional quadratic roughness penalty, one can approximate the local impulse response using an FFT-based calculation. However, that acceleration method cannot be applied directly for separate real and imaginary regularization. This paper proposes a fast and stable calculation for this case that also uses FFT-based calculations to approximate the local impulse responses of the real and imaginary parts. This approach is demonstrated with a quadratic image reconstruction of fMRI data that uses separate roughness penalties for the real and imaginary parts.

Over past several years, machine learning, or more generally artificial intelligence, has generated overwhelming research interest and attracted unprecedented public attention. As tomographic imaging researchers, we share the excitement from our imaging perspective [item 1) in the Appendix], and organized this special issue dedicated to the theme of "Machine learning for image reconstruction." This special issue is a sister issue of the special issue published in May 2016 of this journal with the theme "Deep learning in medical imaging" [item 2) in the Appendix]. While the previous special issue targeted medical image processing/analysis, this special issue focuses on data-driven tomographic reconstruction. These two special issues are highly complementary, since image reconstruction and image analysis are two of the main pillars for medical imaging. Together we cover the whole workflow of medical imaging: from tomographic raw data/features to reconstructed images and then extracted diagnostic features/readings.

This paper describes a new image reconstruction method for digital breast tomosynthesis (DBT). The new method incorporates detector blur into the forward model. The detector blur in DBT causes correlation in the measurement noise. By making a few approximations that are reasonable for breast imaging, we formulated a regularized quadratic optimization problem with a data-fit term that incorporates models for detector blur and correlated noise (DBCN). We derived a computationally efficient separable quadratic surrogate (SQS) algorithm to solve the optimization problem that has a non-diagonal noise covariance matrix. We evaluated the SQS-DBCN method by reconstructing DBT scans of breast phantoms and human subjects. The contrast-to-noise ratio and sharpness of microcalcifications were analyzed and compared with those by the simultaneous algebraic reconstruction technique. The quality of soft tissue lesions and parenchymal patterns was examined. The results demonstrate the potential to improve the image quality of reconstructed DBT images by incorporating the system physics model. This paper is a first step toward model-based iterative reconstruction for DBT.

The development of computed tomography (CT) image reconstruction methods that significantly reduce patient radiation exposure, while maintaining high image quality is an important area of research in low-dose CT imaging. We propose a new penalized weighted least squares (PWLS) reconstruction method that exploits regularization based on an efficient Union of Learned TRAnsforms (PWLS-ULTRA). The union of square transforms is pre-learned from numerous image patches extracted from a dataset of CT images or volumes. The proposed PWLS-based cost function is optimized by alternating between a CT image reconstruction step, and a sparse coding and clustering step. The CT image reconstruction step is accelerated by a relaxed linearized augmented Lagrangian method with ordered-subsets that reduces the number of forward and back projections. Simulations with 2-D and 3-D axial CT scans of the extended cardiac-torso phantom and 3-D helical chest and abdomen scans show that for both normal-dose and low-dose levels, the proposed method significantly improves the quality of reconstructed images compared to PWLS reconstruction with a nonadaptive edge-preserving regularizer. PWLS with regularization based on a union of learned transforms leads to better image reconstructions than using a single learned square transform. We also incorporate patch-based weights in PWLS-ULTRA that enhance image quality and help improve image resolution uniformity. The proposed approach achieves comparable or better image quality compared to learned overcomplete synthesis dictionaries, but importantly, is much faster (computationally more efficient).

We study accelerated dual gradient-based methods for image denoising/deblurring problems based on the total variation (TV) model. For the TV-based denoising problem, combining the dual approach and Nesterov's fast gradient projection (FGP) method has been found effective. The corresponding denoising method minimizes the dual function with FGP's optimal rate O(1/k2) where k denotes the number of iterations, and guarantees a rate O(1/k) for the primal function decrease. Considering that the dual projected gradient decrease is closely related to the primal function decrease, this paper proposes new accelerated gradient projection methods that decrease the projected gradient norm with a fast rate O(1/k1.5) and that are as efficient as FGP. The proposed approach also decreases the primal function with a faster rate O(1/k1.5). We provide preliminary results on image denoising/deblurring problems with a TV regularizer, where the fast and efficient denoising solver is iteratively used for solving a deblurring problem as the inner proximal update of a fast iterative shrinkage/thresholding algorithm (FISTA).

Rapid, reliable quantification of MR relaxation parameters T1and T2is desirable for many clinical applications. Steady-state sequences such as Spoiled Gradient-Recalled Echo (SPGR) and Dual-Echo Steady-State (DESS) are fast and well-suited for relaxometry because the signals they produce are quite sensitive to T1and T2variation. However, T1, T2estimation with these sequences typically requires multiple scans with varied sets of acquisition parameters. This paper describes a systematic framework for selecting scan types (e.g., combinations of SPGR and DESS scans) and optimizing their respective parameters (e.g., flip angles and repetition times). The method is based on a Cramér-Rao Bound (CRB)-inspired min-max optimization that finds scan parameter combinations that robustly enable precise object parameter estimation. We apply this technique to optimize combinations of SPGR and DESS scans for T1, T2relaxometry in white matter (WM) and grey matter (GM) regions of the human brain at 3T field strength. Phantom accuracy experiments show that SPGR/DESS scan combinations are in excellent agreement with reference measurements. Phantom precision experiments show that trends in T1,T2pooled sample standard deviations reflect CRB-based predictions. In vivo experiments show that in WM and GM, T1and T2estimates from a pair of optimized DESS scans exhibit precision (but not necessarily accuracy) comparable to that of optimized combinations of SPGR and DESS scans. To our knowledge, T1maps from DESS acquisitions alone are new. This example application illustrates that scan optimization may help reveal new parameter mapping techniques from combinations of established pulse sequences.

Sparsity-based approaches have been popular in many applications in image processing and imaging. Compressed sensing exploits the sparsity of images in a transform domain or dictionary to improve image recovery fromundersampledmeasurements. In the context of inverse problems in dynamic imaging, recent research has demonstrated the promise of sparsity and low-rank techniques. For example, the patches of the underlying data are modeled as sparse in an adaptive dictionary domain, and the resulting image and dictionary estimation from undersampled measurements is called dictionary-blind compressed sensing, or the dynamic image sequence is modeled as a sum of low-rank and sparse (in some transform domain) components (L+S model) that are estimated from limited measurements. In this work, we investigate a data-adaptive extension of the L+S model, dubbed LASSI, where the temporal image sequence is decomposed into a low-rank component and a component whose spatiotemporal (3D) patches are sparse in some adaptive dictionary domain. We investigate various formulations and efficient methods for jointly estimating the underlying dynamic signal components and the spatiotemporal dictionary from limited measurements. We also obtain efficient sparsity penalized dictionary-blind compressed sensing methods as special cases of our LASSI approaches. Our numerical experiments demonstrate the promising performance of LASSI schemes for dynamicmagnetic resonance image reconstruction from limited k-t space data compared to recent methods such as k-t SLR and L+S, and compared to the proposed dictionary-blind compressed sensing method.

Predicting noise properties of iteratively reconstructed CT images is useful for analyzing reconstruction methods; for example, local noise power spectrum (NPS) predictions may be used to quantify the detectability of an image feature, to design regularization methods, or to determine dynamic tube current adjustment during a CT scan. This paper presents a method for fast prediction of reconstructed image variance and local NPS for statistical reconstruction methods using quadratic or locally quadratic regularization. Previous methods either require impractical computation times to generate an approximate map of the variance of each reconstructed voxel, or are restricted to specific CT geometries. Our method can produce a variance map of the entire image, for locally shift-invariant CT geometries with sufficiently fine angular sampling, using a computation time comparable to a single back-projection. The method requires only the projection data to be used in the reconstruction, not a reconstruction itself, and is reasonably accurate except near image edges where edge-preserving regularization behaves highly nonlinearly. We evaluate the accuracy of our method using reconstructions of both simulated CT data and real CT scans of a thorax phantom.

We propose a new method for the joint design of k-space trajectory and RF pulse in 3D small-tip tailored excitation. Designing time-varying RF and gradient waveforms for a desired 3D target excitation pattern in MRI poses a non-linear, non-convex, constrained optimization problem with relatively large problem size that is difficult to solve directly. Existing joint pulse design approaches are therefore typically restricted to predefined trajectory types such as EPI or stack-of-spirals that intrinsically satisfy the gradient maximum and slew rate constraints and reduce the problem size (dimensionality) dramatically, but lead to suboptimal excitation accuracy for a given pulse duration. Here we use a 2nd-order B-spline basis that can be fitted to an arbitrary k-space trajectory, and allows the gradient constraints to be implemented efficiently. We show that this allows the joint optimization problem to be solved with quite general k-space trajectories. Starting from an arbitrary initial trajectory, we first approximate the trajectory using B-spline basis, and then optimize the corresponding coefficients. We evaluate our method in simulation using four different k-space initializations: stack-of-spirals, SPINS, KT-points, and a new method based on KT-points. In all cases, our approach leads to substantial improvement in excitation accuracy for a given pulse duration. We also validated our method for inner-volume excitation using phantom experiments. The computation is fast enough for online applications.

Statistical image reconstruction (SIR) methods are studied extensively for X-ray computed tomography (CT) due to the potential of acquiring CT scans with reduced X-ray dose while maintaining image quality. However, the longer reconstruction time of SIR methods hinders their use in X-ray CT in practice. To accelerate statistical methods, many optimization techniques have been investigated. Over-relaxation is a common technique to speed up convergence of iterative algorithms. For instance, using a relaxation parameter that is close to two in alternating direction method of multipliers (ADMM) has been shown to speed up convergence significantly. This paper proposes a relaxed linearized augmented Lagrangian (AL) method that shows theoretical faster convergence rate with over-relaxation and applies the proposed relaxed linearized AL method to X-ray CT image reconstruction problems. Experimental results with both simulated and real CT scan data show that the proposed relaxed algorithm (with ordered-subsets [OS] acceleration) is about twice as fast as the existing unrelaxed fast algorithms, with negligible computation and memory overhead.

First-order methods are used widely for large scale optimization problems in signal/image processing and machine learning, because their computation depends mildly on the problem dimension. Nesterov's fast gradient method (FGM) has the optimal convergence rate among first-order methods for smooth convex minimization; its extension to non-smooth case, the fast iterative shrinkage-thresholding algorithm (FISTA), also satisfies the optimal rate; thus both algorithms have gained great interest. We recently introduced a new optimized gradient method (OGM) (for smooth convex functions) having a theoretical convergence speed that is 2× faster than Nesterov's FGM. This paper further discusses the convergence analysis of OGM and explores its fast convergence on an image restoration problem using a smoothed total variation (TV) regularizer. In addition, we empirically investigate the extension of OGM to nonsmooth convex minimization for image restoration with l1-sparsity regularization.

Resolving the 3-D structure of translucent objects is of fundamental interest in combustion research. Recent progress in light-field imaging systems provides a single-shot, single-access-point imaging tool for volumetric image reconstruction. However, when a translucent object being imaged has low contrast edges, the conventional block-matching depth reconstruction algorithm used in commercial plenoptic cameras performs poorly. To solve this problem, this paper investigates reconstructing translucent objects using a model-based image reconstruction (MBIR) method. Preliminary simulation results illustrate that useful 3D object information can be resolved using MBIR methods. However, severe z-direction dispersion is observed in the reconstructed images due to the limited angular variation of incident light rays.

Image denoising is a fundamental operation in image processing, and its applications range from the direct (photographic enhancement) to the technical (as a subproblem in image reconstruction algorithms). In many applications, the number of pixels has continued to grow, while the serial execution speed of computational hardware has begun to stall. New image processing algorithms must exploit the power offered by massively parallel architectures like graphics processing units (GPUs). This paper describes a family of image denoising algorithms well-suited to the GPU. The algorithms iteratively perform a set of independent, parallel 1D pixel-update subproblems. To match GPU memory limitations, they perform these pixel updates in-place and only store the noisy data, denoised image, and problem parameters. The algorithms can handle a wide range of edge-preserving roughness penalties, including differentiable convex penalties and anisotropic total variation. Both algorithms use the majorize-minimize framework to solve the 1D pixel update subproblem. Results from a large 2D image denoising problem and a 3D medical imaging denoising problem demonstrate that the proposed algorithms converge rapidly in terms of both iteration and run-time.

Statistical image reconstruction methods for X-ray computed tomography (CT) provide improved spatial resolution and noise properties over conventional filtered back-projection (FBP) reconstruction, along with other potential advantages such as reduced patient dose and artifacts. Conventional regularized image reconstruction leads to spatially variant spatial resolution and noise characteristics because of interactions between the system models and the regularization. Previous regularization design methods aiming to solve such issues mostly rely on circulant approximations of the Fisher information matrix that are very inaccurate for undersampled geometries like short-scan cone-beam CT. This paper extends the regularization method proposed in to 3-D cone-beam CT by introducing a hypothetical scanning geometry that helps address the sampling properties. The proposed regularization designs were compared with the original method in with both phantom simulation and clinical reconstruction in 3-D axial X-ray CT. The proposed regularization methods yield improved spatial resolution or noise uniformity in statistical image reconstruction for short-scan axial cone-beam CT.

Statistical X-ray computed tomography (CT) reconstruction can improve image quality from reduced dose scans, but requires very long computation time. Ordered subsets (OS) methods have been widely used for research in X-ray CT statistical image reconstruction (and are used in clinical PET and SPECT reconstruction). In particular, OS methods based on separable quadratic surrogates (OS-SQS) are massively parallelizable and are well suited to modern computing architectures, but the number of iterations required for convergence should be reduced for better practical use. This paper introduces OS-SQS-momentum algorithms that combine Nesterov's momentum techniques with OS-SQS methods, greatly improving convergence speed in early iterations. If the number of subsets is too large, the OS-SQS-momentum methods can be unstable, so we propose diminishing step sizes that stabilize the method while preserving the very fast convergence behavior. Experiments with simulated and real 3D CT scan data illustrate the performance of the proposed algorithms.

Sparsity-promoting regularization is useful for combining compressed sensing assumptions with parallel MRI for reducing scan time while preserving image quality. Variable splitting algorithms are the current state-of-the-art algorithms for SENSE-type MR image reconstruction with sparsity-promoting regularization. These methods are very general and have been observed to work with almost any regularizer; however, the tuning of associated convergence parameters is a commonly-cited hindrance in their adoption. Conversely, majorize-minimize algorithms based on a single Lipschitz constant have been observed to be slow in shift-variant applications such as SENSE-type MR image reconstruction since the associated Lipschitz constants are loose bounds for the shift-variant behavior. This paper bridges the gap between the Lipschitz constant and the shift-variant aspects of SENSE-type MR imaging by introducing majorizing matrices in the range of the regularizer matrix. The proposed majorize-minimize methods (called BARISTA) converge faster than state-of-the-art variable splitting algorithms when combined with momentum acceleration and adaptive momentum restarting. Furthermore, the tuning parameters associated with the proposed methods are unitless convergence tolerances that are easier to choose than the constraint penalty parameters required by variable splitting algorithms.

Augmented Lagrangian (AL) methods for solving convex optimization problems with linear constraints are attractive for imaging applications with composite cost functions due to the empirical fast convergence rate under weak conditions. However, for problems such as X-ray computed tomography (CT) image reconstruction, where the inner least-squares problem is challenging and requires iterations, AL methods can be slow. This paper focuses on solving regularized (weighted) least-squares problems using a linearized variant of AL methods that replaces the quadratic AL penalty term in the scaled augmented Lagrangian with its separable quadratic surrogate function, leading to a simpler ordered-subsets (OS) accelerable splitting-based algorithm, OS-LALM. To further accelerate the proposed algorithm, we use a second-order recursive system analysis to design a deterministic downward continuation approach that avoids tedious parameter tuning and provides fast convergence. Experimental results show that the proposed algorithm significantly accelerates the convergence of X-ray CT image reconstruction with negligible overhead and can reduce OS artifacts when using many subsets.

Algorithms with Lipschitz bounds such as ISTA and FISTA are useful for solving optimization problems with sparsity-promoting regularizers. However, they can be slow in applications that involve shift-variant system matrices. One example of such an application is MRI with multiple sensitivity coils. We propose a reconstruction algorithm for wavelet regularized SENSE MR image reconstruction that exploits the spatial localization of the wavelet basis and the shift-variant behavior of the MR system matrix to accelerate algorithm convergence. Our results indicate that the proposed method is faster than state-of-the-art variable splitting algorithms in terms of convergence speed for a SENSE-type reconstruction problem even when the variable splitting methods are tuned carefully. Unlike variable splitting methods, the proposed method requires no convergence parameter tuning.

Fast and accurate quantification of spin-spin relaxation parameter T2is of importance for clinical MRI applications. Classical spin echo (SE) sequences yield straightforward T2estimates, but require undesirably long scans. By contrast, steady-state sequences such as the Dual-Echo Steady-State (DESS) sequence are considerably faster, but produce signals that depend on more complex functions of both desired and nuisance parameters. Conventional method-of-moments estimators exhibit systematic error because of the approximations used to bypass nuisance parameter estimation. To improve T2mapping accuracy, we propose a novel, model-based approach to this nonlinear estimation problem. We use a fast scan to estimate nuisance parameters M0*and T1, and then use the exact DESS signal model for regularized T2estimation from DESS data, with minimal approximations. MR brain simulation results show that the proposed approach substantially improves T2estimation accuracy and precision, compared to conventional method-of-moments estimators.

The Bloch-Siegert (BS) method is a fast method for finding RF transmit fields (|B1+| maps) in MRI based on the phase shift induced by an off-resonance pulse. The phase calculation can be very noisy in regions with low image magnitude, leading to inaccurate estimations using the conventional method. This paper proposes a penalized likelihood (PL) estimator of the |B1+| map in the image domain, and a fast optimization transfer algorithm for solving the problem. We validated our proposed method with both simulation and experimental data sets.

We propose a reconstruction method for the phase retrieval problem prevalent in optics, crystallography, and other imaging applications. Our approach uses signal sparsity to provide robust reconstruction, even in the presence of outliers. Our method is multi-layered, involving multiple random initial conditions, convex majorization, variable splitting, and alternating directions method of multipliers (ADMM)-based implementation. Monte Carlo simulations demonstrate that our algorithm can correctly and robustly detect sparse signals from full and undersampled sets of squared-magnitude-only measurements, corrupted by additive noise or outliers.

The ordered subset expectation maximization (OSEM) algorithm approximates the gradient of a likelihood function using a subset of projections instead of using all projections so that fast image reconstruction is possible for emission and transmission tomography such as SPECT, PET, and CT. However, OSEM does not significantly accelerate reconstruction with computationally expensive regularizers such as patch-based nonlocal (NL) regularizers, because the regularizer gradient is evaluated for every subset. We propose to use variable splitting to separate the likelihood term and the regularizer term for penalized emission tomographic image reconstruction problem and to optimize it using the alternating direction method of multiplier (ADMM). We also propose a fast algorithm to optimize the ADMM parameter based on convergence rate analysis. This new scheme enables more sub-iterations related to the likelihood term. We evaluated our ADMM for 3-D SPECT image reconstruction with a patch-based NL regularizer that uses the Fair potential function. Our proposed ADMM improved the speed of convergence substantially compared to other existing methods such as gradient descent, EM, and OSEM using De Pierro's approach, and the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm.

Spectral computed tomography (CT) provides information on material characterization and quantification because of its ability to separate different basis materials. Dual-energy (DE) CT provides two sets of measurements at two different source energies. In principle, two materials can be accurately decomposed from DECT measurements. However, many clinical and industrial applications require three or more material images. For triple-material decomposition, a third constraint, such as volume conservation, mass conservation or both, is required to solve three sets of unknowns from two sets of measurements. The recently proposed flexible image-domain (ID) multi-material decomposition) method assumes each pixel contains at most three materials out of several possible materials and decomposes a mixture pixel by pixel. We propose a penalized-likelihood (PL) method with edge-preserving regularizers for each material to reconstruct multi-material images using a similar constraint from sinogram data. We develop an optimization transfer method with a series of pixel-wise separable quadratic surrogate (PWSQS) functions to monotonically decrease the complicated PL cost function. The PWSQS algorithm separates pixels to allow simultaneous update of all pixels, but keeps the basis materials coupled to allow faster convergence rate than our previous proposed material- and pixel-wise SQS algorithms. Comparing with the ID method using 2-D fan-beam simulations, the PL method greatly reduced noise, streak and cross-talk artifacts in the reconstructed basis component images, and achieved much smaller root mean square errors.

SPIRiT (iterative self-consistent parallel imaging reconstruction), and its sparsity-regularized variant L1-SPIRiT, are compatible with both Cartesian and non-Cartesian magnetic resonance imaging sampling trajectories. However, the non-Cartesian framework is more expensive computationally, involving a nonuniform Fourier transform with a nontrivial Gram matrix. We propose a novel implementation of the regularized reconstruction problem using variable splitting, alternating minimization of the augmented Lagrangian, and careful preconditioning. Our new method based on the alternating direction method of multipliers converges much faster than existing methods because of the preconditioners' heightened effectiveness. We demonstrate such rapid convergence substantially improves image quality for a fixed computation time. Our framework is a step forward towards rapid non-Cartesian L1-SPIRiT reconstructions.

Parallel excitation requires fast and accurate B1 map estimation. Bloch-Siegert (BS) B1 mapping is very fast and accurate over a large dynamic range. When applied to multi-coil systems, however, this phase-based method may produce low signal-to-noise ratio estimates in low magnitude regions due to localized excitation patterns of parallel excitation systems. Also, the imaging time increases with the number of coils. In this work, we first propose to modify the standard BS B1 mapping sequence so that it avoids the scans required by previous B1 phase estimation methods. A regularized method is then proposed to jointly estimate the magnitude and phase of multi-coil B1 maps from BS B1 mapping data, improving estimation quality by using the prior knowledge of the smoothness of B1 magnitude and phase. Lastly, we use Cramer-Rao lower bound analysis to optimize the coil combinations, to improve the quality of the raw data for B1 estimation. The proposed methods are demonstrated by simulations and phantom experiments.

Statistical image reconstruction methods provide improved image quality in low-dose X-ray CT. However, the long computation time of iterative algorithms limits their clinical use. Ordered subsets algorithms based on separable quadratic surrogates (OS-SQS) are attractive as they are simple and amenable for massive parallelization in modern computing architecture, but require many iterations to converge. Here, we further accelerate OS algorithms by using momentum techniques. We use real patient CT scan to illustrate that the proposed algorithms converge rapidly compared to previous OS algorithms.

To reduce blur in noisy images, regularized image restoration methods have been proposed that use nonquadratic regularizers (like l1 regularization or total-variation) that suppress noise while preserving edges in the image. Most of these methods assume a circulant blur (periodic convolution with a blurring kernel) that can lead to wraparound artifacts along the boundaries of the image due to the implied periodicity of the circulant model. Using a noncirculant model could prevent these artifacts at the cost of increased computational complexity. In this paper, we propose to use a circulant blur model combined with a masking operator that prevents wraparound artifacts. The resulting model is noncirculant, so we propose an efficient algorithm using variable splitting and augmented Lagrangian (AL) strategies. Our variable splitting scheme, when combined with the AL framework and alternating minimization, leads to simple linear systems that can be solved noniteratively using fast Fourier transforms (FFTs), eliminating the need for more expensive conjugate gradient-type solvers. The proposed method can also efficiently tackle a variety of convex regularizers, including edge-preserving (e.g., total-variation) and sparsity promoting (e.g., l1-norm) regularizers. Simulation results show fast convergence of the proposed method, along with improved image quality at the boundaries where the circulant model is inaccurate.

Several magnetic resonance parallel imaging techniques require explicit estimates of the receive coil sensitivity profiles. These estimates must be accurate over both the object and its surrounding regions to avoid generating artifacts in the reconstructed images. Regularized estimation methods that involve minimizing a cost function containing both a data-fit term and a regularization term provide robust sensitivity estimates. However, these methods can be computationally expensive when dealing with large problems. In this paper, we propose an iterative algorithm based on variable splitting and the augmented Lagrangian method that estimates the coil sensitivity profile by minimizing a quadratic cost function. Our method, ADMM-Circ, reformulates the finite differencing matrix in the regularization term to enable exact alternating minimization steps. We also present a faster variant of this algorithm using intermediate updating of the associated Lagrange multipliers. Numerical experiments with simulated and real data sets indicate that our proposed method converges approximately twice as fast as the preconditioned conjugate gradient method over the entire field-of-view. These concepts may accelerate other quadratic optimization problems.

Motion-compensated image reconstruction (MCIR) methods incorporate motion models to improve image quality in the presence of motion. MCIR methods differ in terms of how they use motion information and they have been well studied separately. However, there have been less theoretical comparisions of different MCIR methods. This paper compares the theoretical noise properties of three popular MCIR methods assuming known nonrigid motion. We show the relationship among three MCIR methods-motion-compensated temporal regularization (MTR), the parametric motion model (PMM), and post-reconstruction motion correction (PMC)-for penalized weighted least square cases. These analyses show that PMM and MTR are matrix-weighted sums of all registered image frames, while PMC is a scalar-weighted sum. We further investigate the noise properties of MCIR methods with Poisson models and quadratic regularizers by deriving accurate and fast variance prediction formulas using an "analytical approach." These theoretical noise analyses show that the variances of PMM and MTR are lower than or comparable to the variance of PMC due to the statistical weighting. These analyses also facilitate comparisons of the noise properties of different MCIR methods, including the effects of different quadratic regularizers, the influence of the motion through its Jacobian determinant, and the effect of assuming that total activity is preserved. Two-dimensional positron emission tomography simulations demonstrate the theoretical results.

Compensating for the collimator-detector response (CDR) in SPECT is important for accurate quantification. The CDR consists of both a geometric response and a septal penetration and collimator scatter response. The geometric response can be modeled analytically and is often used for modeling the whole CDR if the geometric response dominates. However, for radionuclides that emit medium or high-energy photons such as I-131, the septal penetration and collimator scatter response is significant and its modeling in the CDR correction is important for accurate quantification. There are two main methods for modeling the depth-dependent CDR so as to include both the geometric response and the septal penetration and collimator scatter response. One is to fit a Gaussian plus exponential function that is rotationally invariant to the measured point source response at several source-detector distances. However, a rotationally-invariant exponential function cannot represent the star-shaped septal penetration tails in detail. Another is to perform Monte-Carlo (MC) simulations to generate the depth-dependent point spread functions (PSFs) for all necessary distances. However, MC simulations, which require careful modeling of the SPECT detector components, can be challenging and accurate results may not be available for all of the different SPECT scanners in clinics. In this paper, we propose an alternative approach to CDR modeling. We use a Gaussian function plus a 2-D B-spline PSF template and fit the model to measurements of an I-131 point source at several distances. The proposed PSF-template-based approach is nearly non-parametric, captures the characteristics of the septal penetration tails, and minimizes the difference between the fitted and measured CDR at the distances of interest. The new model is applied to I-131 SPECT reconstructions of experimental phantom measurements, a patient study, and a MC patient simulation study employing the XCAT phantom. The proposed model yields up to a 16.5 and 10.8% higher recovery coefficient compared to the results with the conventional Gaussian model and the Gaussian plus exponential model, respectively.

Statistical image reconstruction algorithms in X-ray computed tomography (CT) provide improved image quality for reduced dose levels but require substantial computation time. Iterative algorithms that converge in few iterations and that are amenable to massive parallelization are favorable in multiprocessor implementations. The separable quadratic surrogate (SQS) algorithm is desirable as it is simple and updates all voxels simultaneously. However, the standard SQS algorithm requires many iterations to converge. This paper proposes an extension of the SQS algorithm that leads to spatially nonuniform updates. The nonuniform (NU) SQS encourages larger step sizes for the voxels that are expected to change more between the current and the final image, accelerating convergence, while the derivation of NU-SQS guarantees monotonic descent. Ordered subsets (OS) algorithms can also accelerate SQS, provided suitable "subset balance" conditions hold. These conditions can fail in 3-D helical cone-beam CT due to incomplete sampling outside the axial region-of-interest (ROI). This paper proposes a modified OS algorithm that is more stable outside the ROI in helical CT. We use CT scans to demonstrate that the proposed NU-OS-SQS algorithm handles the helical geometry better than the conventional OS methods and "converges" in less than half the time of ordinary OS-SQS.

Magnetic resonance image (MRI) reconstruction from undersampled k-space data requires regularization to reduce noise and aliasing artifacts. Proper application of regularization however requires appropriate selection of associated regularization parameters. In this work, we develop a data-driven regularization parameter adjustment scheme that minimizes an estimate [based on the principle of Stein's unbiased risk estimate (SURE)] of a suitable weighted squared-error measure in k-space. To compute this SURE-type estimate, we propose a Monte-Carlo scheme that extends our previous approach to inverse problems (e.g., MRI reconstruction) involving complex-valued images. Our approach depends only on the output of a given reconstruction algorithm and does not require knowledge of its internal workings, so it is capable of tackling a wide variety of reconstruction algorithms and nonquadratic regularizers including total variation and those based on the l1-norm. Experiments with simulated and real MR data indicate that the proposed approach is capable of providing near mean squared-error optimal regularization parameters for single-coil undersampled non-Cartesian MRI reconstruction.

Image restoration is a well studied problem and there are several proposed methods for deblurring and denoising. Recently, there is increasing interest in iterative schemes that employ non-quadratic regularizers, especially edge-preserving like Total Variation (TV) and sparsity promoting like l1regularization. Most methods make simplifying assumptions concerning the system model and the most common one is the use of a circulant blurring model because it facilitates using the FFT. In this work we focus on a more realistic non-circulant blurring model and apply existing algorithms for image restoration with non-quadratic regularization, tailored to work with our non-circulant model.

Regularization is an effective means of reducing noise and artifacts in MR image reconstruction from undersampled k-space data. Proper application of regularization demands appropriate selection of the associated regularization parameter. Generalized cross-validation (GCV) is a popular parameter tuning technique especially for linear reconstruction methods, but its application to nonlinear iterative MRI reconstruction is more involved as it demands the evaluation of the Jacobian matrix of the reconstruction algorithm with respect to complex-valued data. We derive analytical expressions for recursively updating this Jacobian matrix for an iterative reweighted least-squares reconstruction algorithm. Our method can also be used to calculate a predicted risk estimate (PSURE) for MRI based on Stein's principle. We demonstrate with simulations and experiments with real data that regularization parameter selection based on GCV and PSURE provides near-MSE-optimal results for nonlinear MRI reconstruction from undersampled k-space data using ℓ1-regularization.

Regularized iterative reconstruction algorithms for imaging inverse problems require selection of appropriate regularization parameter values. We focus on the challenging problem of tuning regularization parameters for nonlinear algorithms for the case of additive (possibly complex) Gaussian noise. Generalized cross-validation (GCV) and (weighted) mean-squared error (MSE) approaches (based on Steinfs Unbiased Risk Estimate. SURE) need the Jacobian matrix of the nonlinear reconstruction operator (representative of the iterative algorithm) with respect to the data. We derive the desired Jacobian matrix for two types of nonlinear iterative algorithms: a fast variant of the standard iterative reweighted least-squares method and the contemporary split-Bregman algorithm, both of which can accommodate a wide variety of analysis- and synthesis-type regularizers. The proposed approach iteratively computes two weighted SURE-type measures: Predicted-SURE and Projected-SURE (that require knowledge of noise variance Ð2), and GCV (that does not need Ð2) for these algorithms. We apply the methods to image restoration and to magnetic resonance image (MRI) reconstruction using total variation (TV) and an analysis-type .1-regularization. We demonstrate through simulations and experiments with real data that minimizing Predicted-SURE and Projected-SURE consistently lead to near-MSE-optimal reconstructions. We also observed that minimizing GCV yields reconstruction results that are near-MSE-optimal for image restoration and slightly suboptimal for MRI. Theoretical derivations in this work related to Jacobian matrix evaluations can be extended, in principle, to other types of regularizers and reconstruction algorithms.

Many motion-compensated image reconstruction (MCIR) methods have been proposed to correct for subject motion in medical imaging. MCIR methods incorporate motion models to improve image quality by reducing motion artifacts and noise. This paper analyzes the spatial resolution properties of MCIR methods and shows that nonrigid local motion can lead to nonuniform and anisotropic spatial resolution for conventional quadratic regularizers. This undesirable property is akin to the known effects of interactions between heteroscedastic log-likelihoods (e.g., Poisson likelihood) and quadratic regularizers. This effect may lead to quantification errors in small or narrow structures (such as small lesions or rings) of reconstructed images. This paper proposes novel spatial regularization design methods for three different MCIR methods that account for known nonrigid motion. We develop MCIR regularization designs that provide approximately uniform and isotropic spatial resolution and that match a user-specified target spatial resolution. Two-dimensional PET simulations demonstrate the performance and benefits of the proposed spatial regularization design methods.

Statistical image reconstruction using penalized weighted least-squares (PWLS) criteria can improve image-quality in X-ray computed tomography (CT). However, the huge dynamic range of the statistical weights leads to a highly shift-variant inverse problem making it difficult to precondition and accelerate existing iterative algorithms that attack the statistical model directly. We propose to alleviate the problem by using a variable-splitting scheme that separates the shift-variant and ("nearly") invariant components of the statistical data model and also decouples the regularization term. This leads to an equivalent constrained problem that we tackle using the classical method-of-multipliers framework with alternating minimization. The specific form of our splitting yields an alternating direction method of multipliers (ADMM) algorithm with an inner-step involving a "nearly" shift-invariant linear system that is suitable for FFT-based preconditioning using cone-type filters. The proposed method can efficiently handle a variety of convex regularization criteria including smooth edge-preserving regularizers and nonsmooth sparsity-promoting ones based on the l(1)-norm and total variation. Numerical experiments with synthetic and real in vivo human data illustrate that cone-filter preconditioners accelerate the proposed ADMM resulting in fast convergence of ADMM compared to conventional (nonlinear conjugate gradient, ordered subsets) and state-of-the-art (MFISTA, split-Bregman) algorithms that are applicable for CT.

Compressed sensing (CS) has been used for accelerating magnetic resonance imaging acquisitions, but its use in applications with rapid spatial phase variations is challenging, e.g., proton resonance frequency shift (PRF-shift) thermometry and velocity mapping. Previously, an iterative MRI reconstruction with separate magnitude and phase regularization was proposed for applications where magnitude and phase maps are both of interest, but it requires fully sampled data and unwrapped phase maps. In this paper, CS is combined into this framework to reconstruct magnitude and phase images accurately from undersampled data. Moreover, new phase regularization terms are proposed to accommodate phase wrapping and to reconstruct images with encoded phase variations, e.g., PRF-shift thermometry and velocity mapping. The proposed method is demonstrated with simulated thermometry data and in vivo velocity mapping data and compared to conventional phase corrected CS.

X-ray computed tomography (CT) images could be improved using iterative image reconstruction if the 3D cone-beam forward- and back-projection computations can be accelerated significantly. We investigated the feasibility of a field-programmable gate array (FPGA) implementation of the separable footprint (SF) forward projector. A 16-bit fixed-point quantization introduces negligible numerical errors without affecting the perceptual image quality. The SF-based 3D cone-beam projector can be efficiently parallelized and its memory bandwidth reduced by exploiting projection geometry and data locality. We demonstrate a fully pipelined, 75-way parallel hardware architecture of the SF forward projector on a Xilinx Virtex-5 FPGA that can complete one forward projection of a 320×320×61 object over 3,625 views in 6.3 seconds.

In simultaneous dual-isotope myocardial perfusion SPECT (MPS) imaging, data are simultaneously acquired to determine the distributions of two radioactive isotopes. The goal of this work was to develop penalized maximum likelihood (PML) algorithms for a novel cross-tracer prior that exploits the fact that the two images reconstructed from simultaneous dual-isotope MPS projection data are perfectly registered in space. We first formulated the simultaneous dual-isotope MPS reconstruction problem as a joint estimation problem. A cross-tracer prior that couples voxel values on both images was then proposed. We developed an iterative algorithm to reconstruct the MPS images that converges to the maximum a posteriori solution for this prior based on separable surrogate functions. To accelerate the convergence, we developed a fast algorithm for the cross-tracer prior based on the complete data OS-EM (COSEM) framework. The proposed algorithm was compared qualitatively and quantitatively to a single-tracer version of the prior that did not include the cross-tracer term. Quantitative evaluations included comparisons of mean and standard deviation images as well as assessment of image fidelity using the mean square error. We also evaluated the cross tracer prior using a three-class observer study with respect to the three-class MPS diagnostic task, i.e., classifying patients as having either no defect, reversible defect, or fixed defects. For this study, a comparison with conventional ordered subsets-expectation maximization (OS-EM) reconstruction with postfiltering was performed. The comparisons to the single-tracer prior demonstrated similar resolution for areas of the image with large intensity changes and reduced noise in uniform regions. The cross-tracer prior was also superior to the single-tracer version in terms of restoring image fidelity. Results of the three-class observer study showed that the proposed cross-tracer prior and the convergent algorithms improved the image quality of dual-isotope MPS images compared to OS-EM.

Magnetic resonance image (MRI) reconstruction using SENSitivity Encoding (SENSE) requires regularization to suppress noise and aliasing effects. Edge-preserving and sparsity-based regularization criteria can improve image quality, but they demand computation-intensive nonlinear optimization. In this paper, we present novel methods for regularized MRI reconstruction from undersampled sensitivity encoded data--SENSE-reconstruction--using the augmented Lagrangian (AL) framework for solving large-scale constrained optimization problems. We first formulate regularized SENSE-reconstruction as an unconstrained optimization task and then convert it to a set of (equivalent) constrained problems using variable splitting. We then attack these constrained versions in an AL framework using an alternating minimization method, leading to algorithms that can be implemented easily. The proposed methods are applicable to a general class of regularizers that includes popular edge-preserving (e.g., total-variation) and sparsity-promoting (e.g., l(1)-norm of wavelet coefficients) criteria and combinations thereof. Numerical experiments with synthetic and in vivo human data illustrate that the proposed AL algorithms converge faster than both general-purpose optimization algorithms such as nonlinear conjugate gradient (NCG) and state-of-the-art MFISTA.

Iterative methods for 3D image reconstruction have the potential to improve image quality over conventional filtered back projection (FBP) in X-ray computed tomography (CT). However, the computation burden of 3D cone-beam forward and back-projectors is one of the greatest challenges facing practical adoption of iterative methods for X-ray CT. Moreover, projector accuracy is also important for iterative methods. This paper describes two new separable footprint (SF) projector methods that approximate the voxel footprint functions as 2D separable functions. Because of the separability of these footprint functions, calculating their integrals over a detector cell is greatly simplified and can be implemented efficiently. The SF-TR projector uses trapezoid functions in the transaxial direction and rectangular functions in the axial direction, whereas the SF-TT projector uses trapezoid functions in both directions. Simulations and experiments showed that both SF projector methods are more accurate than the distance-driven (DD) projector, which is a current state-of-the-art method in the field. The SF-TT projector is more accurate than the SF-TR projector for rays associated with large cone angles. The SF-TR projector has similar computation speed with the DD projector and the SF-TT projector is about two times slower.

Nonrigid image registration methods using intensity-based similarity metrics are becoming increasingly common tools to estimate many types of deformations. Nonrigid warps can be very flexible with a large number of parameters and gradient optimization schemes are widely used to estimate them. However, for large datasets, the computation of the gradient of the similarity metric with respect to these many parameters becomes very time consuming. Using a small random subset of image voxels to approximate the gradient can reduce computation time. This work focuses on the use of importance sampling to reduce the variance of this gradient approximation. The proposed importance sampling framework is based on an edge-dependent adaptive sampling distribution designed for use with intensity-based registration algorithms. We compare the performance of registration based on stochastic approximations with and without importance sampling to that using deterministic gradient descent. Empirical results, on simulated magnetic resonance brain data and real computed tomography inhale-exhale lung data from eight subjects, show that a combination of stochastic approximation methods and importance sampling accelerates the registration process while preserving accuracy.

Large-tip-angle multidimensional radio-frequency (RF) pulse design is a difficult problem, due to the nonlinear response of magnetization to applied RF at large tip-angles. In parallel excitation, multidimensional RF pulse design is further complicated by the possibility for transmit field patterns to change between subjects, requiring pulses to be designed rapidly while a subject lies in the scanner. To accelerate pulse design, we introduce a fast version of the optimal control method for large-tip-angle parallel excitation. The new method is based on a novel approach to analytically linearizing the Bloch equation about a large-tip-angle RF pulse, which results in an approximate linear model for the perturbations created by adding a small-tip-angle pulse to a large-tip-angle pulse. The linear model can be evaluated rapidly using nonuniform fast Fourier transforms, and we apply it iteratively to produce a sequence of pulse updates that improve excitation accuracy. We achieve drastic reductions in design time and memory requirements compared to conventional optimal control, while producing pulses of similar accuracy. The new method can also compensate for nonidealities such as main field inhomogeneties.

Dual-energy (DE) X-ray computed tomography (CT) has been found useful in various applications. In medical imaging, one promising application is using low-dose DECT for attenuation correction in positron emission tomography (PET). Existing approaches to sinogram material decomposition ignore noise characteristics and are based on logarithmic transforms, producing noisy component sinogram estimates for low-dose DECT. In this paper, we propose two novel sinogram restoration methods based on statistical models: penalized weighted least square (PWLS) and penalized likelihood (PL), yielding less noisy component sinogram estimates for low-dose DECT than classical methods. The proposed methods consequently provide more precise attenuation correction of the PET emission images than do previous methods for sinogram material decomposition with DECT. We report simulations that compare the proposed techniques and existing approaches.

Statistical methods for tomographic image reconstruction have improved noise and spatial resolution properties that may improve image quality in X-ray computed tomography (CT). Penalized weighted least squares (PWLS) methods using conventional quadratic regularization lead to nonuniform and anisotropic spatial resolution due to interactions between the weighting, which is necessary for good noise properties, and the regularizer. Previously, we addressed this problem for parallel-beam emission tomography using matrix algebra methods to design data-dependent, shift-variant regularizers that improve resolution uniformity. This paper develops a fast angular integral mostly analytical (AIMA) regularization design method for 2-D fan-beam X-ray CT imaging, for which parallel-beam tomography is a special case. Simulation results demonstrate that the new method for regularization design requires very modest computation and leads to nearly uniform and isotropic spatial resolution in transmission tomography when using quadratic regularization.

In fast magnetic resonance (MR) imaging with long readout times, such as echo-planar imaging (EPI) and spiral scans, it is important to correct for the effects of field inhomogeneity to reduce image distortion and blurring. Such corrections require an accurate field map, a map of the off-resonance frequency at each voxel. Standard field map estimation methods yield noisy field maps, particularly in image regions with low spin density. This paper describes regularized methods for field map estimation from two or more MR scans having different echo times. These methods exploit the fact that field maps are generally smooth functions. The methods use algorithms that decrease monotonically a regularized least-squares cost function, even though the problem is highly nonlinear. Results show that the proposed regularized methods significantly improve the quality of field map estimates relative to conventional unregularized methods.

Blood oxygen level dependent (BOLD) functional magnetic resonance imaging (fMRI) is conventionally done by reconstructing T(2)(*)-weighted images. However, since the images are unitless they are nonquantifiable in terms of important physiological parameters. An alternative approach is to reconstruct R(2)(*) maps which are quantifiable and have comparable BOLD contrast as T(2)(*)-weighted images. However, conventional R(2)(*) mapping involves long readouts and ignores relaxation during readout. Another problem with fMRI imaging is temporal drift/fluctuations in off-resonance. Conventionally, a field map is collected at the start of the fMRI study to correct for off-resonance, ignoring any temporal changes. Here, we propose a new fast regularized iterative algorithm that jointly reconstructs R(2)(*) and field maps for all time frames in fMRI data. To accelerate the algorithm we linearize the MR signal model, enabling the use of fast regularized iterative reconstruction methods. The regularizer was designed to account for the different resolution properties of both R(2)(*) and field maps and provide uniform spatial resolution. For fMRI data with the same temporal frame rate as data collected for T(2)(*)-weighted imaging the resulting R(2)(*) maps performed comparably to T(2)(*)-weighted images in activation detection while also correcting for spatially global and local temporal changes in off-resonance.

The majorize-minimize (MM) optimization technique has received considerable attention in signal and image processing applications, as well as in statistics literature. At each iteration of an MM algorithm, one constructs a tangent majorant function that majorizes the given cost function and is equal to it at the current iterate. The next iterate is obtained by minimizing this tangent majorant function, resulting in a sequence of iterates that reduces the cost function monotonically. A well-known special case of MM methods are expectation-maximization algorithms. In this paper, we expand on previous analyses of MM, due to Fessler and Hero, that allowed the tangent majorants to be constructed in iteration-dependent ways. Also, this paper overcomes an error in one of those earlier analyses. There are three main aspects in which our analysis builds upon previous work. First, our treatment relaxes many assumptions related to the structure of the cost function, feasible set, and tangent majorants. For example, the cost function can be nonconvex and the feasible set for the problem can be any convex set. Second, we propose convergence conditions, based on upper curvature bounds, that can be easier to verify than more standard continuity conditions. Furthermore, these conditions allow for considerable design freedom in the iteration-dependent behavior of the algorithm. Finally, we give an original characterization of the local region of convergence of MM algorithms based on connected (e.g., convex) tangent majorants. For such algorithms, cost function minimizers will locally attract the iterates over larger neighborhoods than typically is guaranteed with other methods. This expanded treatment widens the scope of the MM algorithm designs that can be considered for signal and image processing applications, allows us to verify the convergent behavior of previously published algorithms, and gives a fuller understanding overall of how these algorithms behave.

Respiratory motion remains a significant source of errors in treatment planning for the thorax and upper abdomen. Recently, we proposed a method to estimate two-dimensional (2-D) object motion from a sequence of slowly rotating X-ray projection views, which we called deformation from orbiting views (DOVs). In this method, we model the motion as a time varying deformation of a static prior of the anatomy. We then optimize the parameters of the motion model by maximizing the similarity between the modeled and actual projection views. This paper extends the method to full three-dimensional (3-D) motion and cone-beam projection views. We address several practical issues for using a cone-beam computed tomography (CBCT) scanner that is integrated in a radiotherapy system, such as the effects of Compton scatter and the limited gantry rotation for one breathing cycle. We also present simulation and phantom results to illustrate the performance of this method.

Accurate predictions of image variances can be useful for reconstruction algorithm analysis and for the design of regularization methods. Computing the predicted variance at every pixel using matrix-based approximations [1] is impractical. Even most recently adopted methods that are based on local discrete Fourier approximations are impractical since they would require a forward and backprojection and two fast Fourier transform (FFT) calculations for every pixel, particularly for shift-variant systems like fan-beam tomography. This paper describes new "analytical" approaches to predicting the approximate variance maps of 2-D images that are reconstructed by penalized-likelihood estimation with quadratic regularization in fan-beam geometries. The simplest of the proposed analytical approaches requires computation equivalent to one backprojection and some summations, so it is computationally practical even for the data sizes in X-ray computed tomography (CT). Simulation results show that it gives accurate predictions of the variance maps. The parallel-beam geometry is a simple special case of the fan-beam analysis. The analysis is also applicable to 2-D positron emission tomography (PET).

No convergent ordered subsets (OS) type image reconstruction algorithms for transmission tomography have been proposed to date. In contrast, in emission tomography, there are two known families of convergent OS algorithms: methods that use relaxation parameters, and methods based on the incremental expectation-maximization (EM) approach. This paper generalizes the incremental EM approach by introducing a general framework, "incremental optimization transfer." The proposed algorithms accelerate convergence speeds and ensure global convergence without requiring relaxation parameters. The general optimization transfer framework allows the use of a very broad family of surrogate functions, enabling the development of new algorithms. This paper provides the first convergent OS-type algorithm for (nonconcave) penalized-likelihood (PL) transmission image reconstruction by using separable paraboloidal surrogates (SPS) which yield closed-form maximization steps. We found it is very effective to achieve fast convergence rates by starting with an OS algorithm with a large number of subsets and switching to the new "transmission incremental optimization transfer (TRIOT)" algorithm. Results show that TRIOT is faster in increasing the PL objective than nonincremental ordinary SPS and even OS-SPS yet is convergent.

In order to reconstruct attenuation maps with improved spatial resolution and quantitative accuracy, we developed an approximate method of incorporating system resolution compensation (RC) in the ordered-subset transmission (OSTR) algorithm for transmission reconstruction. Our method approximately models the blur caused by the finite intrinsic detector resolution, the nonideal source collimation and detector collimation. We derived the formulation using the optimization transfer principle as in the derivation of the OSTR algorithm. The formulation includes one forward-blur step and one back-blur step, which do not severely slow down reconstruction. The formulation could be applicable to various transmission geometries, such as point-source, line-source, and sheet-source systems. Through computer simulations of the MCAT phantom and transmission measurements of the air-filled Data Spectrum Deluxe single photo emission computed tomography (SPECT) Phantom on a system which employed a cone-beam geometry and a system which employed a scanning-line-source geometry, we showed that incorporation of RC increased spatial resolution and improved the quantitative accuracy of reconstruction. In simulation studies, attenuation maps reconstructed with RC correction improved the quantitative accuracy of emission reconstruction.

Fourier-based forward and back-projection methods can reduce computation in iterative tomographic image reconstruction. Recently, an optimized nonuniform fast Fourier transform (NUFFT) approach was shown to yield accurate parallel-beam projections. In this paper, we extend the NUFFT approach to describe an O (N2 log N) projector/backprojector pair for fan-beam transmission tomography. Simulations and experiments with real CT data show that fan-beam Fourier-based forward and back-projection methods can reduce computation for iterative reconstruction while still providing accuracy comparable to their O (N3) space-based counterparts.

We consider the task of detecting a statistically varying signal of known location on a statistically varying background in a reconstructed tomographic image. We analyze the performance of linear observer models in this task. We show that, if one chooses a suitable reconstruction method, a broad family of linear observers can exactly achieve the optimal detection performance attainable with any combination of a linear observer and linear reconstructor. This conclusion encompasses several well-known observer models from the literature, including models with a frequency-selective channel mechanism and certain types of internal noise. Interestingly, the "optimal" reconstruction methods are unregularized and in some cases quite unconventional. These results suggest that, for the purposes of designing regularized reconstruction methods that optimize lesion detectability, known-location tasks are of limited use.

A new image reconstruction method to correct for the effects of magnetic field inhomogeneity in non-Cartesian sampled magnetic resonance imaging (MRI) is proposed. The conjugate phase reconstruction method, which corrects for phase accumulation due to applied gradients and magnetic field inhomogeneity, has been commonly used for this case. This can lead to incomplete correction, in part, due to the presence of gradients in the field inhomogeneity function. Based on local distortions to the k-space trajectory from these gradients, a spatially variant sample density compensation function is introduced as part of the conjugate phase reconstruction. This method was applied to both simulated and experimental spiral imaging data and shown to produce more accurate image reconstructions. Two approaches for fast implementation that allow the use of fast Fourier transforms are also described. The proposed method is shown to produce fast and accurate image reconstructions for spiral sampled MRI.

For reconstructing a complex object wavefront from digital holograms, we propose a new penalized-likelihood approach based on the measurement statistics and edge-preserving regularization. The log-likelihood is complicated since the measurements are related to the magnitude of the complex beam. We use optimization transfer to derive a new simplified iterative algorithm that monotonically decreases the cost function. Unlike the conventional FFT-based holographic reconstruction method, the new approach uses all of the measured data and can be applied to holograms with any (known) reference beam pattern. Simulation results demonstrate the potential for improved image quality.

Most positron emission tomography (PET) emission scans are corrected for accidental coincidence (AC) events by real-time subtraction of delayed-window coincidences, leaving only the randoms-precorrected data available for image reconstruction. The real-time randoms precorrection compensates in mean for AC events but destroys the Poisson statistics. The exact log-likelihood for randoms-precorrected data is inconvenient, so practical approximations are needed for maximum likelihood or penalized-likelihood image reconstruction. Conventional approximations involve setting negative sinogram values to zero, which can induce positive systematic biases, particularly for scans with low counts per ray. We propose new likelihood approximations that allow negative sinogram values without requiring zero-thresholding. With negative sinogram values, the log-likelihood functions can be nonconcave, complicating maximization; nevertheless, we develop monotonic algorithms for the new models by modifying the separable paraboloidal surrogates and the maximum-likelihood expectation-maximization (ML-EM) methods. These algorithms ascend to local maximizers of the objective function. Analysis and simulation results show that the new shifted Poisson (SP) model is nearly free of systematic bias yet keeps low variance. Despite its simpler implementation, the new SP performs comparably to the saddle-point model which has shown the best performance (as to systematic bias and variance) in randoms-precorrected PET emission reconstruction.

The ordinary sample correlation coefficient is a popular similarity measure for aligning images from the same or similar modalities. However, this measure can be sensitive to the presence of "outlier" objects that appear in one image but not the other, such as surgical instruments, the patient table, etc., which can lead to biased registrations. This paper describes an intensity-based image registration technique that uses a robust correlation coefficient as a similarity measure. Relative to the ordinary sample correlation coefficient, the proposed similarity measure reduces the influence of outliers. We also compared the performance of the proposed method with the mutual information-based method. The robust correlation-based method should be useful for image registration in radiotherapy (KeV to MeV X-ray images) and image-guided surgery applications. We have investigated the properties of the proposed method by theoretical analysis, computer simulations, a phantom experiment, and with functional magnetic resonance imaging data.

Iterative image reconstruction algorithms play an increasingly important role in modern tomographic systems, especially in emission tomography. With the fast increase of the sizes of the tomographic data, reduction of the computation demands of the reconstruction algorithms is of great importance. Fourier-based forward and back-projection methods have the potential to considerably reduce the computation time in iterative reconstruction. Additional substantial speed-up of those approaches can be obtained utilizing powerful and cheap off-the-shelf fast Fourier transform (FFT) processing hardware. The Fourier reconstruction approaches are based on the relationship between the Fourier transform of the image and Fourier transformation of the parallel-ray projections. The critical two steps are the estimations of the samples of the projection transform, on the central section through the origin of Fourier space, from the samples of the transform of the image, and vice versa for back-projection. Interpolation errors are a limitation of Fourier-based reconstruction methods. We have applied min-max optimized Kaiser-Bessel interpolation within the nonuniform FFT (NUFFT) framework and devised ways of incorporation of resolution models into the Fourier-based iterative approaches. Numerical and computer simulation results show that the min-max NUFFT approach provides substantially lower approximation errors in tomographic forward and back-projection than conventional interpolation methods. Our studies have further confirmed that Fourier-based projectors using the NUFFT approach provide accurate approximations to their space-based counterparts but with about ten times faster computation, and that they are viable candidates for fast iterative image reconstruction.

Imaging systems that form estimates using a statistical approach generally yield images with nonuniform resolution properties. That is, the reconstructed images possess resolution properties marked by space-variant and/or anisotropic responses. We have previously developed a space-variant penalty for penalized-likelihood (PL) reconstruction that yields nearly uniform resolution properties. We demonstrated how to calculate this penalty efficiently and apply it to an idealized positron emission tomography (PET) system whose geometric response is space-invariant. In this paper, we demonstrate the efficient calculation and application of this penalty to space-variant systems. (The method is most appropriate when the system matrix has been precalculated.) We apply the penalty to a large field of view PET system where crystal penetration effects make the geometric response space-variant, and to a two-dimensional single photon emission computed tomography system whose detector responses are modeled by a depth-dependent Gaussian with linearly varying full-width at half-maximum. We perform a simulation study comparing reconstructions using our proposed PL approach with other reconstruction methods and demonstrate the relative resolution uniformity, and discuss tradeoffs among estimators that yield nearly uniform resolution. We observe similar noise performance for the PL and post-smoothed maximum-likelihood (ML) approaches with carefully matched resolution, so choosing one estimator over another should be made on other factors like computational complexity and convergence rates of the iterative reconstruction. Additionally, because the postsmoothed ML and the proposed PL approach can outperform one another in terms of resolution uniformity depending on the desired reconstruction resolution, we present and discuss a hybrid approach adopting both a penalty and post-smoothing.

Resolution and covariance predictors have been derived previously for penalized-likelihood estimators. These predictors can provide accurate approximations to the local resolution properties and covariance functions for tomographic systems given a good estimate of the mean measurements. Although these predictors may be evaluated iteratively, circulant approximations are often made for practical computation times. However, when numerous evaluations are made repeatedly (as in penalty design or calculation of variance images), these predictors still require large amounts of computing time. In Stayman and Fessler (2000), we discussed methods for precomputing a large portion of the predictor for shift-invariant system geometries. In this paper, we generalize the efficient procedure discussed in Stayman and Fessler (2000) to shift-variant single photon emission computed tomography (SPECT) systems. This generalization relies on a new attenuation approximation and several observations on the symmetries in SPECT systems. These new general procedures apply to both two-dimensional and fully three-dimensional (3-D) SPECT models, that may be either precomputed and stored, or written in procedural form. We demonstrate the high accuracy of the predictions based on these methods using a simulated anthropomorphic phantom and fully 3-D SPECT system. The evaluation of these predictors requires significantly less computation time than traditional prediction techniques, once the system geometry specific precomputations have been made.

We present two types of globally convergent relaxed ordered subsets (OS) algorithms for penalized-likelihood image reconstruction in emission tomography: modified block sequential regularized expectation-maximization (BSREM) and relaxed OS separable paraboloidal surrogates (OS-SPS). The global convergence proof of the existing BSREM (De Pierro and Yamagishi, 2001) required a few a posteriori assumptions. By modifying the scaling functions of BSREM, we are able to prove the convergence of the modified BSREM under realistic assumptions. Our modification also makes stepsize selection more convenient. In addition, we introduce relaxation into the OS-SPS algorithm (Erdoğan and Fessler, 1999) that otherwise would converge to a limit cycle. We prove the global convergence of diagonally scaled incremental gradient methods of which the relaxed OS-SPS is a special case; main results of the proofs are from (Nedić and Bertsekas, 2001) and (Correa and Lemaréchal, 1993). Simulation results showed that both new algorithms achieve global convergence yet retain the fast initial convergence speed of conventional unrelaxed ordered subsets algorithms.

Regularization is desirable for image reconstruction in emission tomography. A powerful regularization method is the penalized-likelihood (PL) reconstruction algorithm (or equivalently, maximum a posteriori reconstruction), where the sum of the likelihood and a noise suppressing penalty term (or Bayesian prior) is optimized. Usually, this approach yields position-dependent resolution and bias. However, for some applications in emission tomography, a shift-invariant point spread function would be advantageous. Recently, a new method has been proposed, in which the penalty term is tuned in every pixel to impose a uniform local impulse response. In this paper, an alternative way to tune the penalty term is presented. We performed positron emission tomography and single photon emission computed tomography simulations to compare the performance of the new method to that of the postsmoothed maximum-likelihood (ML) approach, using the impulse response of the former method as the postsmoothing filter for the latter. For this experiment, the noise properties of the PL algorithm were not superior to those of postsmoothed ML reconstruction.

In magnetic resonance imaging, magnetic field inhomogeneities cause distortions in images that are reconstructed by conventional fast Fourier trasform (FFT) methods. Several noniterative image reconstruction methods are used currently to compensate for field inhomogeneities, but these methods assume that the field map that characterizes the off-resonance frequencies is spatially smooth. Recently, iterative methods have been proposed that can circumvent this assumption and provide improved compensation for off-resonance effects. However, straightforward implementations of such iterative methods suffer from inconveniently long computation times. This paper describes a tool for accelerating iterative reconstruction of field-corrected MR images: a novel time-segmented approximation to the MR signal equation. We use a min-max formulation to derive the temporal interpolator. Speedups of around 60 were achieved by combining this temporal interpolator with a nonuniform fast Fourier transform with normalized root mean squared approximation errors of 0.07%. The proposed method provides fast, accurate, field-corrected image reconstruction even when the field map is not smooth.

Iterative coordinate ascent algorithms have been shown to be useful for image recovery, but are poorly suited to parallel computing due to their sequential nature. This paper presents a new fast converging parallelizable algorithm for image recovery that can be applied to a very broad class of objective functions. This method is based on paraboloidal surrogate functions and a concavity technique. The paraboloidal surrogates simplify the optimization problem. The idea of the concavity technique is to partition pixels into subsets that can be updated in parallel to reduce the computation time. For fast convergence, pixels within each subset are updated sequentially using a coordinate ascent algorithm. The proposed algorithm is guaranteed to monotonically increase the objective function and intrinsically accommodates nonnegativity constraints. A global convergence proof is summarized. Simulation results show that the proposed algorithm requires less elapsed time for convergence than iterative coordinate ascent algorithms. With four parallel processors, the proposed algorithm yields a speedup factor of 3.77 relative to single processor coordinate ascent algorithms for a three-dimensional (3-D) confocal image restoration problem.

This paper describes a statistical image reconstruction method for X-ray computed tomography (CT) that is based on a physical model that accounts for the polyenergetic X-ray source spectrum and the measurement nonlinearities caused by energy-dependent attenuation. We assume that the object consists of a given number of nonoverlapping materials, such as soft tissue and bone. The attenuation coefficient of each voxel is the product of its unknown density and a known energy-dependent mass attenuation coefficient. We formulate a penalized-likelihood function for this polyenergetic model and develop an ordered-subsets iterative algorithm for estimating the unknown densities in each voxel. The algorithm monotonically decreases the cost function at each iteration when one subset is used. Applying this method to simulated X-ray CT measurements of objects containing both bone and soft tissue yields images with significantly reduced beam hardening artifacts.

Tomographic image reconstruction using statistical methods can provide more accurate system modeling, statistical models, and physical constraints than the conventional filtered backprojection (FBP) method. Because of the ill posedness of the reconstruction problem, a roughness penalty is often imposed on the solution to control noise. To avoid smoothing of edges, which are important image attributes, various edge-preserving regularization methods have been proposed. Most of these schemes rely on information from local neighborhoods to determine the presence of edges. In this paper, we propose a cost function that incorporates nonlocal boundary information into the regularization method. We use an alternating minimization algorithm with deterministic annealing to minimize the proposed cost function, jointly estimating region boundaries and object pixel values. We apply variational techniques implemented using level-sets methods to update the boundary estimates; then, using the most recent boundary estimate, we minimize a space-variant quadratic cost function to update the image estimate. For the positron emission tomography transmission reconstruction application, we compare the bias-variance tradeoff of this method with that of a "conventional" penalized-likelihood algorithm with local Huber roughness penalty.

The FFT is used widely in signal processing for efficient computation of the Fourier transform (FT) over a set of uniformly spaced frequency locations. However, in many applications, one requires nonuniform sampling in the frequency domain, i.e., a nonuniform FT. Several papers have described fast approximations for the nonuniform FT based on interpolating an oversampled FFT. This paper presents a method for the nonuniform FT that is optimal in a min-max sense. The proposed method minimizes the worst-case approximation error over all signals of unit norm. Unlike many previous methods for the nonuniform FT, the proposed method easily generalizes to multidimensional signals. We are investigating this method as a fast algorithm for computing the Radon transform in 2D iterative tomographic image reconstruction.

In positron emission tomography (PET), positron emission from radiolabeled compounds yields two high energy photons emitted in opposing directions. However, often the photons are not detected due to attenuation within the patient. This attenuation is nonuniform and must be corrected to obtain quantitatively accurate emission images. To measure attenuation effects, one typically acquires a PET transmission scan before or after the injection of radiotracer. In commercially available PET scanners, image reconstruction is performed sequentially in two steps regardless of the reconstruction method: 1. Attenuation correction factor computation (ACF) from transmission scans, 2. Emission image reconstruction using the computed ACFs. This two-step reconstruction scheme does not use all the information in the transmission and emission scans. Postinjection transmission scans contain emission contamination that includes information about emission parameters. Similarly, emission scans contain information about the attenuating medium. To use all the available information, we propose a joint estimation approach that estimates the attenuation map and the emission image simultaneously from these two scans. The penalized-likelihood objective function is nonconvex for this problem. We propose an algorithm based on paraboloidal surrogates that alternates between updating emission and attenuation parameters and is guaranteed to monotonically decrease the objective function.

We introduce a new fast converging parallelizable algorithm for image restoration. This algorithm is based on paraboloidal surrogate functions to simplify the optimization problem and a concavity technique developed by De Pierro (1995) to simultaneously update a set of pixels. To obtain large step sizes which affect the convergence rate, we choose the paraboloidal surrogate functions that have small curvatures. The concavity technique is applied to separate pixels into partitioned sets so that parallel processors can be assigned to each set. The partitioned separable paraboloidal surrogates are maximized by using coordinate ascent (CA) algorithms. Unlike other existing algorithms such EM and CA algorithms, the proposed algorithm not only requires less time per iteration to converge, but is guaranteed to monotonically increase the objective function and intrinsically accommodates nonnegativity constraints as well.

We derive the exact statistical distribution of maximum a posteriori (MAP) estimators having edge-preserving nonGaussian priors. Such estimators have been widely advocated for image restoration and reconstruction problems. Previous investigations of these image recovery methods have been primarily empirical; the distribution we derive enables theoretical analysis. The signal model is linear with Gaussian measurement noise. We assume that the energy function of the prior distribution is chosen to ensure a unimodal posterior distribution (for which convexity of the energy function is sufficient), and that the energy function satisfies a uniform Lipschitz regularity condition. The regularity conditions are sufficiently general to encompass popular priors such as the generalized Gaussian Markov random field prior and the Huber prior, even though those priors are not everywhere twice continuously differentiable.

Traditional space-invariant regularization methods in tomographic image reconstruction using penalized-likelihood estimators produce images with nonuniform spatial resolution properties. The local point spread functions that quantify the smoothing properties of such estimators are space-variant, asymmetric, and object-dependent even for space-invariant imaging systems. We propose a new quadratic regularization scheme for tomographic imaging systems that yields increased spatial uniformity and is motivated by the least-squares fitting of a parameterized local impulse response to a desired global response. We have developed computationally efficient methods for PET systems with shift-invariant geometric responses. We demonstrate the increased spatial uniformity of this new method versus conventional quadratic regularization schemes in simulated PET thorax scans.

In many transmission imaging geometries, the transmitted "beams" of photons overlap on the detector, such that a detector element may record photons that originated in different sources or source locations and thus traversed different paths through the object. Examples include systems based on scanning line sources or on multiple parallel rod sources. The overlap of these beams has been disregarded by both conventional analytical reconstruction methods as well as by previous statistical reconstruction methods. We propose a new algorithm for statistical image reconstruction of attenuation maps that explicitly accounts for overlapping beams in transmission scans. The algorithm is guaranteed to monotonically increase the objective function at each iteration. The availability of this algorithm enables the possibility of deliberately increasing the beam overlap so as to increase count rates. Simulated single photon emission tomography transmission scans based on a multiple line source array demonstrate that the proposed method yields improved resolution/noise tradeoffs relative to "conventional" reconstruction algorithms, both statistical and nonstatistical.

Gradient-based iterative methods often converge slowly for tomographic image reconstruction and image restoration problems, but can be accelerated by suitable preconditioners. Diagonal preconditioners offer some improvement in convergence rate, but do not incorporate the structure of the Hessian matrices in imaging problems. Circulant preconditioners can provide remarkable acceleration for inverse problems that are approximately shift-invariant, i.e., for those with approximately block-Toeplitz or block-circulant Hessians. However, in applications with nonuniform noise variance, such as arises from Poisson statistics in emission tomography and in quantum-limited optical imaging, the Hessian of the weighted least-squares objective function is quite shift-variant, and circulant preconditioners perform poorly. Additional shift-variance is caused by edge-preserving regularization methods based on nonquadratic penalty functions. This paper describes new preconditioners that approximate more accurately the Hessian matrices of shift-variant imaging problems. Compared to diagonal or circulant preconditioning, the new preconditioners lead to significantly faster convergence rates for the unconstrained conjugate-gradient (CG) iteration. We also propose a new efficient method for the line-search step required by CG methods. Applications to positron emission tomography (PET) illustrate the method.

In this paper a minimax methodology is presented for combining information from two imaging modalities having different intrinsic spatial resolutions. The focus application is emission computed tomography (ECT), a low-resolution modality for reconstruction of radionuclide tracer density, when supplemented by high-resolution anatomical boundary information extracted from a magnetic resonance image (MRI) of the same imaging volume. The MRI boundary within the two-dimensional (2-D) slice of interest is parameterized by a closed planar curve. The Cramer-Rao (CR) lower bound is used to analyze estimation errors for different boundary shapes. Under a spatially inhomogeneous Gibbs field model for the tracer density a representation for the minimax MRI-enhanced tracer density estimator is obtained. It is shown that the estimator is asymptotically equivalent to a penalized maximum likelihood (PML) estimator with resolution-selective Gibbs penalty. Quantitative comparisons are presented using the iterative space alternating generalized expectation maximization (SAGE-FM) algorithm to implement the PML estimator with and without minimax weight averaging.

We present a framework for designing fast and monotonic algorithms for transmission tomography penalized-likelihood image reconstruction. The new algorithms are based on paraboloidal surrogate functions for the log likelihood. Due to the form of the log-likelihood function it is possible to find low curvature surrogate functions that guarantee monotonicity. Unlike previous methods, the proposed surrogate functions lead to monotonic algorithms even for the nonconvex log likelihood that arises due to background events, such as scatter and random coincidences. The gradient and the curvature of the likelihood terms are evaluated only once per iteration. Since the problem is simplified at each iteration, the CPU time is less than that of current algorithms which directly minimize the objective, yet the convergence rate is comparable. The simplicity, monotonicity, and speed of the new algorithms are quite attractive. The convergence rates of the algorithms are demonstrated using real and simulated PET transmission scans.

This paper analyzes and compares image reconstruction methods based on practical approximations to the exact log likelihood of randoms-precorrected positron emission tomography (PET) measurements. The methods apply to both emission and transmission tomography, however, in this paper we focus on transmission tomography. The results of experimental PET transmission scans and variance approximations demonstrate that the shifted Poisson (SP) method avoids the systematic bias of the conventional data-weighted least squares (WLS) method and leads to significantly lower variance than conventional statistical methods based on the log likelihood of the ordinary Poisson (OP) model. We develop covariance approximations to analyze the propagation of noise from attenuation maps into emission images via the attenuation correction factors (ACF's). Empirical pixel and region variances from real transmission data agree closely with the analytical predictions. Both the approximations and the empirical results show that the performance differences between the OP model and SP model are even larger, when considering noise propagation from the transmission images into the final emission images, than the differences in the attenuation maps themselves.

We present a framework for designing fast and monotonic algorithms for transmission tomography penalized likelihood image reconstruction. The new algorithms are based on paraboloidal surrogate functions for the log-likelihood. Due to the form of the log-likelihood function, it is possible to find low curvature surrogate functions that guarantee monotonicity. Unlike previous methods, the proposed surrogate functions lead to monotonic algorithms even for the nonconvex log-likelihood that arises due to background events such as scatter and random coincidences. The gradient and the curvature of the likelihood terms are evaluated only once per iteration. Since the problem is simplified, the CPU time per iteration is less than that of current algorithms which directly minimize the objective, yet the convergence rate is comparable. The simplicity, monotonicity and speed of the new algorithms are quite attractive. The convergence rates of the algorithms are demonstrated using real PET transmission scans.

Traditional space-invariant regularization schemes in tomographic image reconstruction using penalized likelihood estimators produce images with nonuniform resolution properties. The local point spread functions that quantify the local smoothing properties of such estimators are not only space-variant and asymmetric, but are also object-dependent even for space-invariant systems. We propose a new regularization scheme for increased spatial uniformity and demonstrate the resolution properties of this new method versus conventional regularization schemes through an investigation of local point spread functions.

We propose a new objective function for the image reconstruction problem, where the image is comprised of piecewise smooth regions separated by sharp boundaries. We use alternating minimization to minimize our objective function. We use the level set technique to minimize with regard to the boundary. The advantage of this new approach is shown through the bias/variance analysis of a hot spot.

The author comments that Zhang (see ibid., vol.5, no.12, p.1662-65, 1996) attempts to establish the convergence of a mean field iteration for an Ising Markov random field (MRF) for large values of the hyperparameter value /spl beta/ is not correct for the function T/sub /spl delta//(u) defined. The convergence of the mean field iteration remains an open question for large /spl beta/.

In this paper, a method is introduced for incorporating perfectly registered MRI boundary information into a penalized likelihood emission reconstruction scheme. The boundary curve is modeled as a periodic spline whose coefficients are estimated from the MRI image. The resulting boundary estimate is mapped to a spatially variant set of Gibbs weights. When incorporated into a quadratic roughness penalty, these weights improve emission reconstruction bias/variance performance by preventing smoothing across the estimated boundary. Finally, we derive a new penalty function that accounts for the uncertainty inherent in the boundary estimates.

Preconditioning methods can accelerate the convergence of gradient-based iterative methods for tomographic image reconstruction and image restoration. Circulant preconditioners have been used extensively for shift-invariant problems. Diagonal preconditioners offer some improvement in convergence rate, but do not incorporate the structure of the Hessian matrices in imaging problems. For inverse problems that are approximately shift-invariant (i.e. approximately block-Toeplitz or block-circulant Hessians), circulant or Fourier-based preconditioners can provide remarkable acceleration. However, in applications with nonuniform noise variance (such as arises from Poisson statistics in emission tomography and in quantum-limited optical imaging), the Hessian of the (penalized) weighted least-squares objective function is quite shift-variant, and the Fourier preconditioner performs poorly. Additional shift-variance is caused by edge-preserving regularization methods based on nonquadratic penalty functions. This paper describes new preconditioners that more accurately approximate the Hessian matrices of shift-variant imaging problems. Compared to diagonal or Fourier preconditioning, the new preconditioners lead to significantly faster convergence rates for the unconstrained conjugate-gradient (CG) iteration. Applications to position emission tomography (PET) illustrate the method.

This paper presents a new class of algorithms for penalized-likelihood reconstruction of attenuation maps from low-count transmission scans. We derive the algorithms by applying to the transmission log-likelihood a version of the convexity technique developed by De Pierro for emission tomography. The new class includes the single-coordinate ascent (SCA) algorithm and Lange's convex algorithm for transmission tomography as special cases. The new grouped-coordinate ascent (GCA) algorithms in the class overcome several limitations associated with previous algorithms. 1) Fewer exponentiations are required than in the transmission maximum likelihood-expectation maximization (ML-EM) algorithm or in the SCA algorithm. 2) The algorithms intrinsically accommodate nonnegativity constraints, unlike many gradient-based methods. 3) The algorithms are easily parallelizable, unlike the SCA algorithm and perhaps line-search algorithms. We show that the GCA algorithms converge faster than the SCA algorithm, even on conventional workstations. An example from a low-count positron emission tomography (PET) transmission scan illustrates the method.

This paper presents a method for incorporating anatomical NMR boundary side information into penalized maximum likelihood (PML) emission image reconstructions. The NMR boundary is parameterized as a periodic spline curve of fixed order and number of knots that is known a priori. Maximum likelihood (ML) estimation of the spline coefficients yields an "extracted" boundary, which is used to define a set of Gibbs weights on the emission image space. These weights, when coupled with a quadratic penalty function, create an edge-preserving penalty that incorporates our prior knowledge effectively. Qualitative analysis demonstrates that our method results in smooth images that do not suffer loss of edge contrast, while quantitative estimates of bias and variance for various values of the smoothing parameter show an improvement over standard quadratically penalized maximum likelihood.

Many estimators in signal processing problems are defined implicitly as the maximum of some objective function. Examples of implicitly defined estimators include maximum likelihood, penalized likelihood, maximum a posteriori, and nonlinear least squares estimation. For such estimators, exact analytical expressions for the mean and variance are usually unavailable. Therefore, investigators usually resort to numerical simulations to examine the properties of the mean and variance of such estimators. This paper describes approximate expressions for the mean and variance of implicitly defined estimators of unconstrained continuous parameters. We derive the approximations using the implicit function theorem, the Taylor expansion, and the chain rule. The expressions are defined solely in terms of the partial derivatives of whatever objective function one uses for estimation. As illustrations, we demonstrate that the approximations work well in two tomographic imaging applications with Poisson statistics. We also describe a "plug-in" approximation that provides a remarkably accurate estimate of variability even from a single noisy Poisson sinogram measurement. The approximations should be useful in a wide range of estimation problems.

This paper examines the spatial resolution properties of penalized-likelihood image reconstruction methods by analyzing the local impulse response. The analysis shows that standard regularization penalties induce space-variant local impulse response functions, even for space-invariant tomographic systems. Paradoxically, for emission image reconstruction, the local resolution is generally poorest in high-count regions. We show that the linearized local impulse response induced by quadratic roughness penalties depends on the object only through its projections. This analysis leads naturally to a modified regularization penalty that yields reconstructed images with nearly uniform resolution. The modified penalty also provides a very practical method for choosing the regularization parameter to obtain a specified resolution in images reconstructed by penalized-likelihood methods.

Many estimators in signal processing problems are defined implicitly as the maximum of an objective function, such as maximum likelihood (ML) and maximum a posteriori (MAP) methods. Exact analytical expressions for the mean and variance of such estimators are usually unavailable, so investigators usually resort to numerical simulations. The paper describes approximate analytical expressions for the mean and variance of implicitly defined estimators. The expressions are defined solely in terms of the partial derivatives of whatever objective function one uses for estimation. The authors demonstrate the utility and accuracy of the approximations in a PET transmission computed tomography application with Poisson statistics. The approximations should be useful in a wide range of estimation problems.

We give estimation error bounds and specify optimal estimators for continuous, closed boundary curves in an NMR image. The boundary is parameterized using periodic B-splines. A Cramer-Rao lower bound on mean-square-estimate error in the presence of system smoothing and Gaussian noise is derived, and the performance of maximum likelihood and penalized maximum likelihood estimators is compared to this bound. Finally, we comment on the usefulness of estimates of the boundary for providing anatomical side information in the reconstruction of functional tomographic images like those of a PET or SPECT system.

Single photon emission computed tomographic images (SPECT) have relatively poor resolution. In an attempt to improve SPECT image quality, many methods have been developed for including anatomic information, extracted from higher resolution, structurally correlated magnetic resonance images (MRI), into the SPECT reconstruction process. These methods provide improved SPECT reconstruction accuracy if the anatomic information is perfectly correlated with the SPECT functional information. However there exist mismatches between MRI anatomical structures and SPECT functional structures due to different imaging mechanisms. It has been reported that if the MR structures are applied into SPECT, the mismatched part will cause artifacts. The paper describes a joint estimation approach which unifies MR information extraction and SPECT reconstruction processes to avoid these artifacts. Both qualitative and quantitative evaluations show that the method improves the SPECT reconstruction where the MR information matches and is robust to mismatched MR information.

Iterative methods for tomographic image reconstruction often converge slowly. Preconditioning methods can often accelerate gradient-based iterations. Previous preconditioning methods for PET reconstruction have used either diagonal or Fourier-based preconditioners. Fourier-based preconditioners are well suited to problems with near-circulant Hessian matrices. However, due to the nonuniform Poisson noise variance in PET, the circulant approximation to the Hessian is suboptimal. This paper shows that a particular combined diagonal/Fourier preconditioner yields a more accurate approximation to the Hessian and gives significantly faster convergence rates than does either preconditioner used alone.

This paper describes rapidly converging algorithms for computing attenuation maps from Poisson transmission measurements using penalized-likelihood objective functions. We demonstrate that an under-relaxed cyclic coordinate-ascent algorithm converges faster than the convex algorithm of Lange (see ibid., vol.4, no.10, p.1430-1438, 1995), which in turn converges faster than the expectation-maximization (EM) algorithm for transmission tomography. To further reduce computation, one could replace the log-likelihood objective with a quadratic approximation. However, we show with simulations and analysis that the quadratic objective function leads to biased estimates for low-count measurements. Therefore we introduce hybrid Poisson/polynomial objective functions that use the exact Poisson log-likelihood for detector measurements with low counts, but use computationally efficient quadratic or cubic approximations for the high-count detector measurements. We demonstrate that the hybrid objective functions reduce computation time without increasing estimation bias.

Most expectation-maximization (EM) type algorithms for penalized maximum-likelihood image reconstruction converge slowly, particularly when one incorporates additive background effects such as scatter, random coincidences, dark current, or cosmic radiation. In addition, regularizing smoothness penalties (or priors) introduce parameter coupling, rendering intractable the M-steps of most EM-type algorithms. This paper presents space-alternating generalized EM (SAGE) algorithms for image reconstruction, which update the parameters sequentially using a sequence of small "hidden" data spaces, rather than simultaneously using one large complete-data space. The sequential update decouples the M-step, so the maximization can typically be performed analytically. We introduce new hidden-data spaces that are less informative than the conventional complete-data space for Poisson data and that yield significant improvements in convergence rate. This acceleration is due to statistical considerations, not numerical overrelaxation methods, so monotonic increases in the objective function are guaranteed. We provide a general global convergence proof for SAGE methods with nonnegativity constraints.

This paper reviews and compares three maximum likelihood algorithms for transmission tomography. One of these algorithms is the EM algorithm, one is based on a convexity argument devised by De Pierro (see IEEE Trans. Med. Imaging, vol.12, p.328-333, 1993) in the context of emission tomography, and one is an ad hoc gradient algorithm. The algorithms enjoy desirable local and global convergence properties and combine gracefully with Bayesian smoothing priors. Preliminary numerical testing of the algorithms on simulated data suggest that the convex algorithm and the ad hoc gradient algorithm are computationally superior to the EM algorithm. This superiority stems from the larger number of exponentiations required by the EM algorithm. The convex and gradient algorithms are well adapted to parallel computing.

The paper summarizes a new Bayesian method for edge-preserving image restoration from noisy measurements. The line-site method of Geman and Geman (1984) forces region boundaries to lie along pixel boundaries, which is unnatural, particularly for 3D data. The present authors augment the intensity process with a binary “mixture site” process, which has one parameter for each pixel indicating the presence of a boundary at some unknown location within that pixel. The method was motivated by the PET and SPECT transmission images with partial volume effects, and is easily extended to 3D data sets.

We apply a uniform Cramer-Rao (CR) bound to study the bias-variance trade-offs in parameter estimation. The uniform CR bound is used to specify achievable and unachievable regions in the bias-variance trade-off plane. The applications considered are: (1) two-dimensional single photon emission computed tomography (SPECT) system, and (2) one dimensional edge localization.>

We give a recursive algorithm to calculate submatrices of the Cramer-Rao (CR) matrix bound on the covariance of any unbiased estimator of a vector parameter /spl theta/_. Our algorithm computes a sequence of lower bounds that converges monotonically to the CR bound with exponential speed of convergence. The recursive algorithm uses an invertible "splitting matrix" to successively approximate the inverse Fisher information matrix. We present a statistical approach to selecting the splitting matrix based on a "complete-data-incomplete-data" formulation similar to that of the well-known EM parameter estimation algorithm. As a concrete illustration we consider image reconstruction from projections for emission computed tomography.>

The authors develop a strategy for joint estimation of physiological parameters and myocardial boundaries using ECT (emission computed tomography). They construct an observation model to relate parameters of interest to the projection data and to account for limited ECT system resolution and measurement noise. The authors then use a maximum likelihood (ML) estimator to jointly estimate all the parameters directly from the projection data without reconstruction of intermediate images. They also simulate myocardial perfusion studies based on a simplified heart model to evaluate the performance of the model-based joint ML estimator and compare this performance to the Cramer-Rao lower bound. Finally, the authors discuss model assumptions and potential uses of the joint estimation strategy.

The authors have previously developed a model-based strategy for joint estimation of myocardial perfusion and boundaries using ECT (emission computed tomography). They have also reported difficulties with boundary estimation in low contrast and low count rate situations. Here they propose using boundary side information (obtainable from high resolution MRI and CT images) or boundary regularization to improve both perfusion and boundary estimation in these situations. To fuse boundary side information into the emission measurements, the authors formulate a joint log-likelihood function to include auxiliary boundary measurements as well as ECT projection measurements. In addition, they introduce registration parameters to align auxiliary boundary measurements with ECT measurements and jointly estimate these parameters with other parameters of interest from the composite measurements. In simulated PET O-15 water myocardial perfusion studies using a simplified model, the authors show that the joint estimation improves perfusion estimation performance and gives boundary alignment accuracy of <0.5 mm even at 0.2 million counts. They implement boundary regularization through formulating a penalized log-likelihood function. They also demonstrate in simulations that simultaneous regularization of the epicardial boundary and myocardial thickness gives comparable perfusion estimation accuracy with the use of boundary side information.

Presents an image reconstruction method for positron-emission tomography (PET) based on a penalized, weighted least-squares (PWLS) objective. For PET measurements that are precorrected for accidental coincidences, the author argues statistically that a least-squares objective function is as appropriate, if not more so, than the popular Poisson likelihood objective. The author proposes a simple data-based method for determining the weights that accounts for attenuation and detector efficiency. A nonnegative successive over-relaxation (+SOR) algorithm converges rapidly to the global minimum of the PWLS objective. Quantitative simulation results demonstrate that the bias/variance tradeoff of the PWLS+SOR method is comparable to the maximum-likelihood expectation-maximization (ML-EM) method (but with fewer iterations), and is improved relative to the conventional filtered backprojection (FBP) method. Qualitative results suggest that the streak artifacts common to the FBP method are nearly eliminated by the PWLS+SOR method, and indicate that the proposed method for weighting the measurements is a significant factor in the improvement over FBP.

An estimation approach is described for three-dimensional reconstruction from line integral projections using incomplete and very noisy data. Generalized cylinders parameterized by stochastic dynamic models are used to represent prior knowledge about the properties of objects of interest in the probed domain. The object models, a statistical measurement model, and the maximum a posteriori probability performance criterion are combined to reformulate the reconstruction problem as a computationally challenging nonlinear estimation problem. For computational feasibility, a suboptimal hierarchical algorithm is described whose individual steps are locally optimal and are combined to satisfy a global optimality criterion. The formulation and algorithm are restricted to objects whose center axis is a single-valued function of a fixed spatial coordinate. Simulation examples demonstrate accurate reconstructions with as few as four views in a 135 degrees sector, at an average signal-to-noise ratio of 3.3.>