Cone-beam and fan-beam image reconstruction algorithms based on spherical and circular harmonics

A cone-beam image reconstruction algorithm using spherical harmonic expansions is proposed. The reconstruction algorithm is in the form of a summation of inner products of two discrete arrays of spherical harmonic expansion coefficients at each cone-beam point of acquisition. This form is different from the common filtered backprojection algorithm and the direct Fourier reconstruction algorithm. There is no re-sampling of the data, and spherical harmonic expansions are used instead of Fourier expansions. As a special case, a new fan-beam image reconstruction algorithm is also derived in terms of a circular harmonic expansion. Computer simulation results for both cone-beam and fan-beam algorithms are presented for circular planar orbit acquisitions. The algorithms give accurate reconstructions; however, the implementation of the cone-beam reconstruction algorithm is computationally intensive. A relatively efficient algorithm is proposed for reconstructing the central slice of the image when a circular scanning orbit is used.     

Cone-beam reconstruction for the two-circles-plus-one-line trajectory

The Kodak Image Station In-Vivo FX has an x-ray module with cone-beam configuration for radiographic imaging but lacks the functionality of tomography. To introduce x-ray tomography into the system, we choose the two-circles-plus-one-line trajectory by mounting one translation motor and one rotation motor. We establish a reconstruction algorithm by applying the M-line reconstruction method. Numerical studies and preliminary physical phantom experiment demonstrate the feasibility of the proposed design and reconstruction algorithm.     

Excitation-resolved fluorescence tomography with simplified spherical harmonics equations

Fluorescence tomography (FT) reconstructs the three-dimensional (3D) fluorescent reporter probe distribution inside biological tissue. These probes target molecules of biological function, e.g. cell surface receptors or enzymes, and emit fluorescence light upon illumination with an external light source. The fluorescence light is detected on the tissue surface and a source reconstruction algorithm based on the simplified spherical harmonics (SP(N)) equations calculates the unknown 3D probe distribution inside tissue. While current FT approaches require multiple external sources at a defined wavelength range, the proposed FT method uses only a white light source with tunable wavelength selection for fluorescence stimulation and further exploits the spectral dependence of tissue absorption for the purpose of 3D tomographic reconstruction. We will show the feasibility of the proposed hyperspectral excitation-resolved fluorescence tomography method with experimental data. In addition, we will demonstrate the performance and limitations of such a method under ideal and controlled conditions by means of a digital mouse model and synthetic measurement data. Moreover, we will address issues regarding the required amount of wavelength intervals for fluorescent source reconstruction. We will explore the impact of assumed spatially uniform and nonuniform optical parameter maps on the accuracy of the fluorescence source reconstruction. Last, we propose a spectral re-scaling method for overcoming the observed limitations in reconstructing accurate source distributions in optically non-uniform tissue when assuming only uniform optical property maps for the source reconstruction process.     

Pulsed EPR imaging: image reconstruction using selective acquisition sequences.

Pulsed nuclear magnetic resonance imaging (NMRI) apparatus has developed very quickly. On the other hand, the development of apparatus for pulsed electron paramagnetic resonance imaging (EPRI) has been very slow. This fact is due to the extremely reduced relaxation times of the paramagnetic probes. EPR linewidths are larger than typical NMR linewidths. These large linewidths are also responsible for a substantial worsening of spatial resolution. Due to the brevity of the electronic relaxation times, not all the acquisition/reconstruction techniques currently used in NMRI (such as spin-echo, gradient-echo, etc) can be applied in pulsed EPRI. In fact, the usable sequences in pulsed EPRI are only acquisitions from projections, where it is possible to use stationary magnetic field gradients. Moreover, the use of high fixed magnetic field gradients induces a short decay time constant T2*. The low T2* value can make it impossible for the analogue to digital conversion system (ADC) to reproduce signal variations during the whole acquisition interval and the resolution can worsen. A new pulsed EPRI acquisition sequence from projections, based on selective reception, is presented that is particularly useful in solving the problems of worsening of spatial resolution associated with the use of an ADC. In order to demonstrate the capabilities of our acquisition method, simulated numerical tests will also be reported.     

Analyser-based phase contrast image reconstruction using geometrical optics

Analyser-based phase contrast imaging can provide radiographs of exceptional contrast at high resolution (<100 microm), whilst quantitative phase and attenuation information can be extracted using just two images when the approximations of geometrical optics are satisfied. Analytical phase retrieval can be performed by fitting the analyser rocking curve with a symmetric Pearson type VII function. The Pearson VII function provided at least a 10% better fit to experimentally measured rocking curves than linear or Gaussian functions. A test phantom, a hollow nylon cylinder, was imaged at 20 keV using a Si(1 1 1) analyser at the ELETTRA synchrotron radiation facility. Our phase retrieval method yielded a more accurate object reconstruction than methods based on a linear fit to the rocking curve. Where reconstructions failed to map expected values, calculations of the Takagi number permitted distinction between the violation of the geometrical optics conditions and the failure of curve fitting procedures. The need for synchronized object/detector translation stages was removed by using a large, divergent beam and imaging the object in segments. Our image acquisition and reconstruction procedure enables quantitative phase retrieval for systems with a divergent source and accounts for imperfections in the analyser.     

Photoacoustic image reconstruction in an attenuating medium using singular-value decomposition

Attenuation effects can be significant in photoacoustic tomography since the generated pressure signals are broadband, and ignoring them may lead to image artifacts and blurring. La Rivie?re et al. [Opt. Lett. 31(6), pp. 781-783, (2006)] had previously derived a method for modeling the attenuation effect and correcting for it in the image reconstruction. This was done by relating the ideal, unattenuated pressure signals to the attenuated pressure signals via an integral operator. We derive an integral operator relating the attenuated pressure signals to the absorbed optical energy for a planar measurement geometry. The matrix operator relating the two quantities is a function of the temporal frequency, attenuation coefficient and the two-dimensional spatial frequency. We perform singular-value decomposition (SVD) of this integral operator to study the problem further. We find that the smallest singular values correspond to wavelet-like eigenvectors in which most of the energy is concentrated at times corresponding to greater depths in tissue. This allows us to characterize the ill-posedness of recovering the absorbed optical energy distribution at different depths in an attenuating medium. This integral equation can be inverted using standard SVD methods, and the initial pressure distribution can be recovered. We conduct simulations and derive an algorithm for image reconstruction using SVD for a planar measurement geometry. We also study the noise and resolution properties of this image-reconstruction method.     

Subsets and overrelaxation in iterative image reconstruction.

A number of iterative image reconstruction algorithms were integrated into one formula characterizing each algorithm by only two parameters: overrelaxation and number of subsets. From the formula it follows that the ordered-subsets iteration (OS-EM) is equivalent to iteration with overrelaxation, where the OS level corresponds to the overrelaxation parameter. Algorithms represented by the formula were studied with respect to speed of convergence and image characteristics. In particular, OS-EM was compared with a single-projection iteration procedure using an optimized sequence of overrelaxation parameters (HOSP) which combines rapid convergence with reduced storage requirements. As a result, OS-EM with a constant number of subsets either needed more iteration steps than HOSP or provoked additional noise, depending on the number of subsets used during iteration. OS-EM can be improved by using decreasing OS levels, imitating the decreasing overrelaxation parameters used for HOSP. The resulting OS-EM may be slightly more rapid than HOSP, due to the increasing number of projections used simultaneously.     

Cone-beam computed tomography with a flat-panel imager: effects of image lag

A system for cone-beam computed tomography (CBCT) has been developed based upon the technology of active matrix flat-panel imagers (FPIs), and the system has demonstrated the potential for fully three-dimensional volumetric imaging with high spatial and contrast resolution. This paper investigates the effects of image lag (arising from charge trapping and release in the FPI pixels) upon CBCT reconstructions. Hypotheses were derived based upon a simple, geometrical/physical model, suggesting that image lag in the projection data results primarily in two artifacts: a spatial blurring artifact in the direction opposite to the direction of rotation (called a "comet") and a line artifact along the direction of the first few projections (called a "streak"). The hypotheses were tested by means of computer simulations and experimental measurements that yielded CBCT images of a simple cylindrical water phantom containing an attenuating rod of varying size and composition. The computer simulations generated projection images based upon analysis of the system geometry and a simple model of the FPI that allowed free adjustment of the image lag. Experimental measurements involved CBCT scans of the phantom under various conditions and modes of acquisition followed by examination of the resulting CBCT axial slices for lag artifacts. Measurements were performed as a function of exposure level, position and contrast of the rod, and for three modes of acquisition designed to isolate and/or minimize the two hypothesized artifacts. The results clearly illustrate the comet and streak artifacts, particularly in relation to high-contrast objects imaged at high exposure levels. The significance of such artifacts under clinical conditions is expected to be small, considering the magnitude of the effect relative to the morphology and composition of typical anatomy. The artifacts may become appreciable, however, in the presence of high-contrast objects, such as marker BBs, dental fillings, and metal prosthetics. A procedural method of reducing lag artifacts is demonstrated.     

Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units

Diffuse optical tomographic image reconstruction uses advanced numerical models that are computationally costly to be implemented in the real time. The graphics processing units (GPUs) offer desktop massive parallelization that can accelerate these computations. An open-source GPU-accelerated linear algebra library package is used to compute the most intensive matrix-matrix calculations and matrix decompositions that are used in solving the system of linear equations. These open-source functions were integrated into the existing frequency-domain diffuse optical image reconstruction algorithms to evaluate the acceleration capability of the GPUs (NVIDIA Tesla C 1060) with increasing reconstruction problem sizes. These studies indicate that single precision computations are sufficient for diffuse optical tomographic image reconstruction. The acceleration per iteration can be up to 40, using GPUs compared to traditional CPUs in case of three-dimensional reconstruction, where the reconstruction problem is more underdetermined, making the GPUs more attractive in the clinical settings. The current limitation of these GPUs in the available onboard memory (4 GB) that restricts the reconstruction of a large set of optical parameters, more than 13,377.     

Improved compressed sensing-based cone-beam CT reconstruction using adaptive prior image constraints

Volumetric cone-beam CT (CBCT) images are acquired repeatedly during a course of radiation therapy and a natural question to ask is whether CBCT images obtained earlier in the process can be utilized as prior knowledge to reduce patient imaging dose in subsequent scans. The purpose of this work is to develop an adaptive prior image constrained compressed sensing (APICCS) method to solve this problem. Reconstructed images using full projections are taken on the first day of radiation therapy treatment and are used as prior images. The subsequent scans are acquired using a protocol of sparse projections. In the proposed APICCS algorithm, the prior images are utilized as an initial guess and are incorporated into the objective function in the compressed sensing (CS)-based iterative reconstruction process. Furthermore, the prior information is employed to detect any possible mismatched regions between the prior and current images for improved reconstruction. For this purpose, the prior images and the reconstructed images are classified into three anatomical regions: air, soft tissue and bone. Mismatched regions are identified by local differences of the corresponding groups in the two classified sets of images. A distance transformation is then introduced to convert the information into an adaptive voxel-dependent relaxation map. In constructing the relaxation map, the matched regions (unchanged anatomy) between the prior and current images are assigned with smaller weight values, which are translated into less influence on the CS iterative reconstruction process. On the other hand, the mismatched regions (changed anatomy) are associated with larger values and the regions are updated more by the new projection data, thus avoiding any possible adverse effects of prior images. The APICCS approach was systematically assessed by using patient data acquired under standard and low-dose protocols for qualitative and quantitative comparisons. The APICCS method provides an effective way for us to enhance the image quality at the matched regions between the prior and current images compared to the existing PICCS algorithm. Compared to the current CBCT imaging protocols, the APICCS algorithm allows an imaging dose reduction of 10-40 times due to the greatly reduced number of projections and lower x-ray tube current level coming from the low-dose protocol.     

Statistical image reconstruction for transmission tomography using relaxed ordered subset algorithms

Statistical reconstruction methods offer possibilities for improving image quality as compared to analytical methods, but current reconstruction times prohibit routine clinical applications in x-ray computed tomography (CT). To reduce reconstruction times, we have applied (under) relaxation to ordered subset algorithms. This enables us to use subsets consisting of only single projection angle, effectively increasing the number of image updates within an entire iteration. A second advantage of applying relaxation is that it can help improve convergence by removing the limit cycle behaviour of ordered subset algorithms, which normally do not converge to an optimal solution but rather a suboptimal limit cycle consisting of as many points as there are subsets. Relaxation suppresses the limit cycle behaviour by decreasing the stepsize for approaching the solution. A simulation study for a 2D mathematical phantom and three different ordered subset algorithms shows that all three algorithms benefit from relaxation: equal noise-to-resolution trade-off can be achieved using fewer iterations than the conventional algorithms, while a lower minimal normalized mean square error (NMSE) clearly indicates a better convergence. Two different schemes for setting the relaxation parameter are studied, and both schemes yield approximately the same minimal NMSE.     

Application of the modified spherical harmonics method to some problems in biomedical optics

It is demonstrated that a semi-analytical method, modified spherical harmonics (MSH), is a fast and rigorous solution for the radiative transport equation in an infinite medium containing an isotropic source. The optical properties of the medium can be determined using an inverse scheme based on fluence data (including phase if a modulated source is used) without the drawbacks inherent in the diffusion approximation. The MSH method can also be used to solve for the fluorescence fluence due to a point excitation source in an infinite medium.     

Optimization of 2D image reconstruction for positron emission mammography using IDL

The Clear-PEM system is a prototype machine for Positron Emission Mammography (PEM) under development within the Portuguese PET-Mammography consortium. We have embedded 2D image reconstruction algorithms implemented in IDL within the prototype's image analysis package. The IDL implementation of these algorithms proved to be accurate and computationally efficient. In this paper, we present the implementation of the MLEM, OSEM and ART 2D iterative image reconstruction algorithms for PEM using IDL. C and IDL implementations are compared using realistic Monte Carlo simulated data. We show that IDL can be used for the easy implementation of image reconstruction algorithms for emission tomography.     

Internet2-based 3D PET image reconstruction using a PC cluster

We describe an approach to fast iterative reconstruction from fully three-dimensional (3D) PET data using a network of PentiumIII PCs configured as a Beowulf cluster. To facilitate the use of this system, we have developed a browser-based interface using Java. The system compresses PET data on the user's machine, sends these data over a network, and instructs the PC cluster to reconstruct the image. The cluster implements a parallelized version of our preconditioned conjugate gradient method for fully 3D MAP image reconstruction. We report on the speed-up factors using the Beowulf approach and the impacts of communication latencies in the local cluster network and the network connection between the user's machine and our PC cluster.     

In vitro 3D reconstruction of long bones using B-scan image processing

An echographic image processing method has been developed, and validated by in vitro experiments, for the 3D reconstruction of the long bones of the newborn. The reconstruction of successive parallel cross-sections is obtained by a 2D reconstruction technique using radial B-scan image processing. The automatic segmentation of all the calculated images allows the extraction of the external contours of the skeleton. After structuring the explored volume using a contour association method, a contour interpolation step is required to solve the anisotropy problem, to obtain a 3D representation with cubic voxel lists. The results are encouraging, and a new mechanical part prototype of the acquisition system is under test for in vivo experiments. The main originality of the paper lies in the combination of different steps to obtain a practical solution to a clinical problem.     

Three-dimensional reconstruction of cranial defect using active contour model and image registration

In neurosurgery, cranial incisions during craniotomy can be recovered by cranioplasty—a surgical operation using cranial implants to repair skull defects. However, surgeons often encounter difficulties when grafting prefabricated cranial plates into defective areas, since a perfect match to the cranial incision is difficult to achieve. Previous studies using mirroring technique, surface interpolation, or deformed template had limitations in skull reconstruction to match the patient’s original appearance. For this study, we utilized low-resolution and high-resolution computed tomography images from the patient to repair skull defects, whilst preserving the original shape. Since the accuracy of skull reconstruction was associated with the partial volume effects in the low-resolution images and the percentage of the skull defect in the high-resolution images, the low-resolution images with intact skull were resampled and thresholded followed by active contour model to suppress partial volume artifacts. The resulting low-resolution images were registered with the high-resolution ones, which exhibited different percentages of cranial defect, to extract the incised cranial part. Finally, mesh smoothing refined the three-dimensional model of the cranial defect. Simulation results indicate that the reconstruction was 93.94% accurate for a 20% skull material removal, and 97.76% accurate for 40% skull material removal. Experimental results demonstrate that the proposed algorithm effectively creates a customized implant, which can readily be used in cranioplasty.     

Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach

The reconstruction of images in photoacoustic tomography is reliant on specifying the speed of sound within the propagation medium. However, for in vivo imaging, this value is not normally accurately known. Here, an autofocus approach for automatically selecting the sound speed is proposed. This is based on maximizing the sharpness of the reconstructed image as quantified by a focus function. Several focus functions are investigated, and their performance is discussed. The method is demonstrated using phantom measurements made in a medium with a known sound speed and in vivo measurements of the vasculature in the flank of an adult mouse.