We show an iterative reconstruction framework for diffraction ultrasound tomography.
Abstract
We show an iterative reconstruction framework for diffraction ultrasound tomography. The use of broadband illumination allows significantly reduce the number of projections compared to straight ray tomography. The proposed algorithm makes use of forward non-uniform fast Fourier transform (NUFFT) for iterative Fourier inversion. Incorporation of total variation regularization allows reduce noise and Gibbs phenomena whilst preserving the edges.
Introduction
Ultrasound tomography with diffracting sources is an important type of acoustic imaging. Since the used wavelengths are comparable to the object feature dimensions, wave phenomena such as diffraction become significant. Consequently, the straight ray tomography theory is no more applicable.
The analog of the Fourier Slice Theorem used in straight ray tomography is the Fourier Diffraction Theorem. Using this theorem, image reconstruction in diffraction tomography can be considered as a problem of signal reconstruction from non-uniform frequency samples.
Reconstruction methods used beforehand addressed the problem as straightforward approximation of the inverse non-uniform Fourier transform (NUFT) and involved frequency interpolation, which is liable to introduce significant inaccuracies. Recently, fast and accurate approximations of the forward NUFFT were introduced. Inverse NUFFT can be achieved iteratively in this framework. We adopt this approach for iterative reconstruction in diffraction tomography, combining it with total variation regularization in order to suppress noise preserving the sharpness of edges.
The problem
Tomography is a method of obtaining an image of a slice of the object from projections. A projection is obtained by illuminating the object from a certain angle and measuring the traversing radiation. In diffraction ultrasound tomography, the object is illuminated with a plane acoustic wave and the forward scattered field is measured on a line of detectors. As in straight ray tomography, changing the orientation of the incident plane waves, it is possible to acquire projections at different angles. Unlike conventional tomography, incident wave frequency can also be changed.
When the wavelength of the incident wave is of the same order as the typical feature size, wave phenomena cannot be neglected. Particularly, diffraction has a significant effect on the interaction of the radiation with matter, and the model of line integrals and the Radon transform is not applicable in diffraction tomography.
Reconstruction is possible by the solution of the inverse scattering problem. Unlike straight ray tomography, the preferable strategy is the solution of the inverse problem in frequency domain. Under the assumption of weakly scattering inhomogeneities in the object, the Fourier Diffraction Theorem relates the Fourier transform of the measured scattered field projection with the Fourier transform of the object:
The Fourier Diffraction Theorem:
Given a projection 
of the forward scattered field obtained by illuminating an object f with a plane wave as shown in Figure 1(a), the following equation holds:
Figure 1 – (a) acquisition of a single projection in ultrasound diffraction tomography, (b) illustration of the Fourier Diffraction Theorem
Since wave phenomena obey the superposition principle, illuminating the object with a wave consisting of a set of frequencies (referred to as broadband illumination), rather than a monochromatic wave, will produce samples along a set of semi-circular arcs with different radii. By taking advantage of this fact, one can achieve sufficient image quality with fewer projections.
Two main reconstruction strategies in diffraction tomography are interpolation in frequency domain (analogous to the direct Fourier inversion in straight ray tomography) and interpolation in space domain (analogous to the filtered backprojection), usually termed as backpropagation. However, unlike straight ray tomography, interpolation in space domain is computationally extensive, and thus the majority of efficient algorithms are based on frequency domain interpolation.A common way for image reconstruction in frequency domain is the gridding algorithm. The non-uniform data is interpolated to a uniform Cartesian grid using, for example, polynomial interpolation. Afterwards, the inverse Fourier transform is efficiently computed using FFT. However, this approach is liable to introduce inaccuracies and is sensitive to the configuration of the sample points.
The solution
The heart of iterative image reconstruction from non-uniform frequency samples is the forward non-uniform fast Fourier transform. To define the NUFFT problem, we first consider a one-dimensional case.
Let 
be a vector of non-uniformly distributed frequencies and 
be a vector of samples of a signal. The non-uniform Fourier transform is defined by the formula
In matrix notation
where 
is a full column rank matrix containing discrete exponent functions in its rows
Fast approximation 
of the NUFT operator can be achieved by projecting the signal f on some oversampled uniform Fourier basis 
using standard FFT, with consequent efficient interpolation:
where 
denotes the interpolation operator, which makes use of p neighboring uniform samples for approximation of each non-uniform sample and q is the oversampling factor. Selecting the interpolation coefficients is a problem per se, and there are many ways of doing it. Fessler and Sutton proposed obtaining such interpolation coefficients that minimize the maximum approximation error at a given point of the non-uniform grid over all signals with unit Euclidean norm. This approach can be formulated as a min-max problem:
where 
is the non-zero part of the k-th row of the interpolation matrix 
and 
is a part of the overcomplete DFT basis 
, containing p nearest neighbors of the non-uniform basis element 
.
This problem has an analytic solution:
One of the advantages of the min-max approach is the fact that it allows efficient computation of the adjoint operator, crucial for iterative optimization algorithms.
Straightforward solution of the inverse problem of is a computationally extensive operation. It is given by the Moore-Penrose pseudoinverse:
Alternatively, the inverse problem can be reformulated as an optimization problem:
It is possible to add penalty on some kind of signal irregularity to the object function:
where 
is a parameter controlling the influence of the penalty. Such a penalty approach is usually referred to as Tikhonov regularization. This problem can be solved iteratively using various optimization techniques, which require efficient computation of the objective function and its gradient.
Empirical observations show that the majority of images that occur in nature, and particularly in medical imaging applications, belong to the class of functions of bounded total variation. The penalty term for total variation can be used in problem above. For discrete image, the smoothed total variation is given by
Total variation regularization removes small oscillations (resulting from noise and Gibbs phenomena), without significantly affecting the edges.
Simulation results
The Shepp-Logan analytic head phantom was used for reconstruction example. Eight broadband projections (containing 8 frequencies each) were simulated. The standard gridding algorithm involving frequency interpolation was compared to iterative reconstruction. Non-uniform frequency samples were interpolated on a 64×64 uniform Cartesian grid and then inverted by IFFT.
Figure 2 depicts the reconstruction results. It can be seen that unlike iterative reconstruction, the use of gridding introduces significant artifacts. Interpolation in frequency domain appears to be highly influenced by the regularity of frequency sampling.
Figure 2 – (a) original band-unlimited phantom, (b) reconstruction by gridding, (c) iterative least-squares reconstruction, (d) iterative reconstruction with total variation penalty
When introducing the total variation penalty term, small artifacts of oscillatory nature are first to disappear; however, too strong penalty is liable to affect small features of the image. The influence of total variation regularization is especially significant in presence of noise. Since ultrasound images usually suffer from low SNR, this aspect is of high importance. Figure 3 shows images reconstructed from data contaminated by additive Gaussian noise with different variance. One can observe that the total variation penalty improves the image quality, whilst preserving the edges.
Figure 3 – Iterative reconstruction in presence of noise (a, b: SNR = 20dB, c, d: SNR = 10dB) with total variation penalty (b, d) and without penalty (a, c)
Conclusions
The presented method is capable of taking advantage of broadband illumination and requires fewer projections for image reconstruction. It appeared significantly more accurate compared to the common gridding method. Total variation regularization further improves the image quality, suppressing noise and Gibbs phenomena, thus allowing obtain sufficient reconstruction quality from noisy data.
Possible developments of the iterative approach can be generalization of the regularizer for other classes of signals, incorporation of different NUFFT implementations and non-smooth optimization techniques for efficient minimization of functions with non-smooth penalty terms.
Acknowledgment
We are especially grateful to our project supervisor Dr. Michael Zibulevsky for his help and guidance throughout this work, without which it could not be possible and to the VISL engineer Johanan Erez for his patience and assistance. In addition, we would like to thank: Prof. Jeffrey Fessler (The University of Michigan) who attracted our attention to his novel approach of NUFFT computation and kindly provided his code, Dr. Haim Azhari (Department of Biomedical Engineering, Technion) and Prof. Yehoshua Y. Zeevi for their valuable comments.
We are grateful to the Ollendorf Minerva Center Fund for supporting this project and to all the people not mentioned here, who helped us throughout the work.