Computational and Mathematical Methods in Medicine

Volume 2014 (2014), Article ID 329350, 13 pages

http://dx.doi.org/10.1155/2014/329350

Research Article

Sparse-View Ultrasound Diffraction Tomography Using Compressed Sensing with Nonuniform FFT

Shaoyan Hua, Mingyue Ding, and Ming Yuchi

Image Processing and Intelligence Control Key Laboratory of Education Ministry of China, Department of Biomedical Engineering, School of Life Science and Technology, Huazhong University of Science and Technology, Wuhan 430074, China

Received 22 December 2013; Revised 16 March 2014; Accepted 19 March 2014; Published 24 April 2014

Academic Editor: William Crum

Copyright © 2014 Shaoyan Hua et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Accurate reconstruction of the object from sparse-view sampling data is an appealing issue for ultrasound diffraction tomography (UDT). In this paper, we present a reconstruction method based on compressed sensing framework for sparse-view UDT. Due to the piecewise uniform characteristics of anatomy structures, the total variation is introduced into the cost function to find a more faithful sparse representation of the object. The inverse problem of UDT is iteratively resolved by conjugate gradient with nonuniform fast Fourier transform. Simulation results show the effectiveness of the proposed method that the main characteristics of the object can be properly presented with only 16 views. Compared to interpolation and multiband method, the proposed method can provide higher resolution and lower artifacts with the same view number. The robustness to noise and the computation complexity are also discussed.

1. Introduction

In recent years, ultrasound diffraction tomography (UDT) has drawn more and more attention in medical imaging field. Different from traditional B-mode ultrasound technique which displays the strength of the echoes with gray scale to show anatomic structure, UDT infers the distribution of acoustic properties such as refractivity, attenuation, and density. Since these acoustic properties of normal and diseased tissues have different value ranges [1], UDT has the potential for providing functional information of the object. For example, in the breast cancer exam, the malignant tumor, the benign mass, and the normal tissue can be differentiated by UDT [2].

Under the assumption of weak scatting, the Fourier diffraction theory (FDT) [3, 4] is adopted for the image reconstruction of UDT. Firstly, the object is illuminated by plane sound wave from one certain direction and the scattering waves are measured and sampled. Secondly, the spatial Fourier transform is performed on the stored data. Under the Born or Rytov approximation [5], the corresponding frequency values are considered to represent the 2D Fourier transform of the object nonuniform distributed along a semicircle arc. The above processes are implemented oriented at various angles around the object to acquire sufficient spatial frequency information. Finally, the UDT image is reconstructed through inverse spatial Fourier transform. To avoid undersampling, normally more spatial frequency samples are required [6]; thus more views from different directions are needed. However, the scan and process time and the associated cost will increase. Besides that, it will impose rigorous requirements on controlling precision of the imaging system. Moreover, redundant information is introduced by multiple views, which results in the waste of the system resources. Hence, accurate reconstruction from nonuniform distributed frequency samples in sparse-view situation has great practical significance.

Currently, in addition to reconstructing the underlying object through beamforming [7], there are other two main approaches based on FDT to recover the complete information from the sparse-view data: interpolation method and iterative method. The former calculates the frequency values on the rectangular grids by a predetermined interpolation function based on sampled frequency on the circular arc grid [8, 9]. Although the interpolation method is computationally efficient, it is liable to introduce coordinate conversion error, which can severely introduce artifacts and distort the image, especially in sparse-view situation. The latter repeatedly corrects the reconstructed image by minimizing the inconsistency between the sampled and the estimated frequency. Through designing appropriate optimization strategy, this process can maximally reduce the reconstruction error and suppress the artifacts. Various researches have been done on the UDT image reconstruction with the iterative method under sparse-view situation. Bronstein et al. [10] proposed an iterative reconstruction framework based on nonuniform fast Fourier transform (NUFFT) and reduced the number of view with broadband sound wave. LaRoque et al. [11] introduced a kind of diffraction tomography with few-view (sparse-view) and limited-angle data through total variation (TV) minimization algorithm for absorberless media. Tingting Li et al. [12] combined TV regularization with the iterative next-neighbor regridding (INNG) algorithm to suppress artifacts and noise.

Compressed sensing (CS) [13, 14] is built on the sparse nature of real signals in certain transform domain. CS has the ability to reconstruct signals with samples which are much less than those required by the Nyquist criterion. CS brings great innovations in image reconstruction and has been widely used in medical imaging field such as MR [15], CT [16, 17], PET [18], and photoacoustic imaging [19, 20]. In recent years, some progress has been made in applying this emerging theory in ultrasound imaging field to reduce the amount of the data and the complexity of the imaging system. Schiffner et al. [21–23] investigated the performance of CS in solving the inverse scattering problem in pule-echo diagnostic ultrasound imaging under the constraint that the scatterer distribution is sparse. With the same assumption, Wagner et al. [24–26] proposed a method based on finite rate of innovation and Xampling for the reconstruction of the beamformed image from channel RF data. Shen et al. [27, 28] presented a measurement-domain adaptive beamforming approach based on distributed CS to reconstruct an image of sparse targets. Achim et al. [29] introduced a framework based on CS for ultrasonic signals reconstruction under the assumption of RF echoes with -stable distributions. Liebgott et al. [30–32] studied the feasibility of CS for the reconstruction of channel RF data, the quantity of channel RF data is reduced through introducing the wave atoms as a representation basis for prebeamformed RF signal. Quinsac et al. [33–37] and Dobigeon et al. [38] applied CS theory to recover beamformed 2D RF images using conjugate gradient descent or Bayesian approach. Zobly et al. [39, 40] and Richy et al. [41] applied CS for doppler imaging.

In this paper, we propose an iterative reconstruction method based on CS for sparse-view UDT. The above researchers mainly focused on applying CS to the pulse-echo imaging systems. To the best of our knowledge, CS was rarely applied in the investigation of UDT image reconstruction, and this is the main motivation of this work. According to FDT, the sample points are bounded in a circle with the radius of (: wave number) [9] for transmission UDT; this results in a sharp sampling cutoff in spatial frequency space. Thus, the reconstructed object is a low-pass version of the original and the quality of image will be distorted by the Gibbs aliasing. This problem can be further deteriorated in the sparse-view situation due to limited sample data. In this context, besides the general norm constraint, we introduce TV penalty into the cost function, which can reduce the oscillation and preserve edges of the object [42]. Empirical observations [10] showed that the majority of nature images, particularly medical images, demonstrated piecewise continuous behavior; that is, parts of anatomy structures were supposed to show uniform characteristics, which belonged to the class of functions of bounded TV [43]. Furthermore, to improve computational efficiency, a fast and accurate NUFFT is adopted to calculate the forward scatter field.

Through the numerical simulation, the proposed method is verified and compared with interpolation method and iteration method with broadband signal [10]. The quantitative evaluation and the robustness to noise of the method are discussed. Simulation shows that the object can be faithfully reconstructed in sparse-view situation without noticeable loss of image quality and the reconstructed error is reduced.

The rest of the paper is organized as follows. The basic principle of UDT and CS theory are described in Section 2. The proposed method is presented in Section 3. Section 4 describes the simulation experiments for UDT image reconstruction and results. Section 5 provides discussion about noise robustness and computational complexity of the proposed method. Conclusions and future work are summarized in Section 6.

2. Background

2.1. UDT Based on FDT

The classical 2D UDT imaging configuration is shown in Figure 1. The inhomogeneous object with a distribution function , , is surrounded by homogeneous medium such as water; the is refractive index. The reconstruction object of UDT is to infer the unknown through transmitted signals measured by ultrasound transducer.

Figure 1: Scheme of 2D ultrasound diffraction tomography: illuminating the object with plane wave at an angle and the scattered measured along .

Assume the object is illuminated by a monochromatic plane wave with wave number and angular frequency at an angle . The wave is scattered within and at the boundary of the inhomogeneous object. Set up a rotated cartesian coordinate system such that -axis is coincident with the view angle as shown in Figure 1. The relation between and can be derived through coordinate transformation: , .

The sound pressure field satisfies the following wave equation [4]:

where denotes the Laplacian operator. can be modeled as the superposition of the incident wave and the scatter wave : .

Under the assumption of weak scattering, , or first order Born, or Rytov approximation, one can derive the so-called FDT which relates the scattered field measured along the line to the object by Fourier transform:

where is the Fourier transform of received scattered data with respect to :

and , , and is the complex amplitude of the illuminating plane wave. The quantity of in (2) is known or measurable, and the integral contains the underlying object that we want to reconstruct. Detailed analysis reveals that is the Fourier transform of the object along a semicircle of radius and centered at [4], as the arc AOB depicted in Figure 2, and is the unit vector along direction . Mathematically, the image reconstruction by (2) is a typical inverse problem which can be expressed as the following formulation:

where is the sample data in -space, is the diffraction operator, and is the object which we want to reconstruct. For formulation (4), classical interpolation method based on FDT is not applicable in sparse-view situation due to violating Nyquist limitation. With fast and accurate nonuniform fast Fourier transform, we propose a CS framework for UDT reconstruction in sparse-view situation to improve the image quality of the reconstruction.

Figure 2: Fourier diffraction projection theorem: the Fourier transform of the scattered data equals the Fourier transform of the object along semicircle AOB.

2.2. Compressed Sensing: A Short Overview

Compressed sensing is a kind of signal processing technique for efficiently acquiring and reconstructing a signal. CS has been increasingly adopted in a variety of applications by applied mathematicians, computer scientists, and engineers since it was initiated in 2006 [44]. The idea of CS can be expressed by the following linear measurement model:

where is the unknown signal such as an image that we want to reconstruct, is the measured signal, and is one measure matrix, which is decided by the imaging system. In this work, is the Fourier transform operator and is the corresponding Fourier transform of the received scatter field. Assume the unknown signal can be sparsely represented in terms of a known basis:

where is the basis, is the corresponding representation coefficients. Here, the sparse means that the number of nonzero coefficients of is small.

Substituting (6) into (5), we obtain (7) as follows:

where ; since , the matrix is not invertible. Taking advantage of the sparsity of , CS theory shows that if the meets the so-called restricted isometry property (RIP) condition, we can exactly recover with overwhelming probability by solving the following minimization problem [45, 46]:

where norm counts the number of no-zero entry of a vector. Equation (8) seeks the sparsest one among all the possible solutions of .

However, (8) belongs to the class of NP-hard problem, which is difficult to obtain solutions for nearly all real applications. One computationally tractable alternative for (8) is to solve the following problem [13, 45]:

where .

The above CS theory requires to be sparse and exactly, whereas, in most practical situations, the object is approximately sparse or compressible. Here, approximately sparse means that contains a small number of components with magnitudes significantly larger than those of the rest, which are not necessarily zero; compressible means the coefficients of decay exponentially in absolute value. To address these problems, Candès et al. [47] extended (9) to the following form:

where and represents energy bound of error.

Considering that (10) is a convex optimization problem, it can be further recast as the following regularization equation [48]:

where is a regularization parameter.

3. Compressed Sensing for UDT Imaging Reconstruction

In this section, we present an image reconstruction method for UDT in 2D case, while it can be readily extended to 3D. Figure 3 provides the comparison between the interpolation method and the proposed scheme, which is composed of three major steps. Firstly, the object is illuminated from random angles and the number of views can be far below that restricted by the Nyquist limitation. The measured sparse-view data are processed by Fourier transform to obtain the corresponding spatial frequency samples along semicircular arcs oriented at the view angles. Secondly, the inverse problem (4) is formulated within CS framework by building the measure matrix and exploiting the sparsity of the object. Thirdly, the object is reconstructed through the NUFFT and CG algorithm.

(a)

(b)

Figure 3: Schematic of the frequency interpolation method (a) and proposed method (b). In (a), the object is reconstructed through IFFT after frequency interpolation. In (b), the object is reconstructed through CG and NUFFT.

3.1. Sparse-View Data Sampling

As we discussed in Section 2.2, when the sensing matrix satisfies the RIP [45, 46], the NP-hard inverse problem (8) can be transformed to computationally tractable norm minimization problem ((9) or (11)). However, even for moderate dimensional operators , it is computationally impractical to verify the RIP. Fortunately, a few classes of matrices are shown to hold RIP for almost certainly. It is shown in [49, 50] that, when is a Gaussian or partial Fourier, that is, the entries of are randomly selected using a Gaussian pdf or rows of are randomly selected from the rows of Fourier matrix, satisfies RIP.

For UDT scanner, the sample operator coincides with the matrix mentioned above since the samples in spatial frequency domain are obtained through FDT, Which makes it possible to reconstruct the object through minimization with sparse data. Here, we generate the sparse-view data by randomly choosing view angles and the number of view is much less than that required by the conventional interpolation method.

3.2. Inverse Problem Formulation under CS Framework

The sparse-view data sampling in Section 3.1 generates the measured value in problem (9). To formulate the UDT inverse problem under CS framework, the sensing matrix in (9) must be constructed. Here, is dependent on the imaging principle of UDT, while is the basis for the object to be reconstructed. The explicit expression of and for UDT will be derived as below.

Since the object has limited physical dimensions, we assume the object or in cartesian coordinate system has bounded support ; that is, when or . Let , () be the discrete form of the underlying object that we want to reconstruct, where is the sample period for - and -axis. and are the Fourier transform of and , respectively. Assume the frequency response of is band-limited; that is, there exists a cutoff frequency such that , when or . In practice, the loss of resolution by this band limit is negligible; the reconstructed imaging quantity is more influenced by other factors such as the aperture sizes of transmitting and receiving transducers [9]. Based on the Nyquist theorem, if , for and .

Since the receiving array of UDT has a limited number of elements, we denote the measured discrete field , where is the number of elements and is pitch (the distance between the centers of two adjacent elements). is the spatial sampling interval of UDT system. and are the Fourier transform of and , respectively. for all as . The relation between and can be formulated as the following equations:

where , are derived from the coordinate transform and (2); let arc grid points and () be on half arc AO, OB, respectively; then we can get , , . For more details about this transformation please refer to [4].

According to (2) and (12), if sufficiently approximates and , , then we have . Generally, this requirement can be satisfied in the real UDT system because the pitch of normal medical ultrasound probe is among ( is the sound wave length).

Let be the finite set of view angles. In UDT system, is obtained over . Then the corresponding can be calculated offline for all (3). Therefore, the can be reconstructed through the set of observation :

This set of equations can be written in matrix form: , where , , , , , , , or is approximation error. With (13), in inverse problem (9) is obtained.

CS theory utilizes the sparse nature of the object and reconstructs the object by minimizing the corresponding norm in transform domain. According to FDT, the spatial frequency samples are distributed along the arc AOB (Figure 2). As the incident wave revolves around the object, the AOB describes a disk of radius centered at the origin; that is, the reconstructed object is a low-pass version of the original. Besides that, we also exploit the fact that the structural morphology of human soft tissue is expected to demonstrate piecewise continuous behavior. That means the object belongs to the class of bounded TV [43] and the gradient of the object is sparse. In this work, the sparsity of the underlying object is exploited not only through wavelet transform which provides sparse representations for rapidly varied regions, but also by TV which affords sparse transformation for piecewise smooth object. Furthermore, TV constraint can help to suppress Gibbs effect and preserve edges [10].

For a discrete object , TV is defined as

where , denote the forward finite difference operator in horizontal () and vertical () coordinates, respectively. Combining norm with TV constraints, we extend the problem (11) to the following minimization problem:

where , are two positive regularization parameters, is spatial frequency samples, is the selected basis for sparse representation of the object , and is the coefficients of in basis .

3.3. Object Reconstruction

An iterative method based on CG is adopted to solve the inverse problem (15). For UDT system, the receive elements are equally spaced that means the along is equally spaced sampled. However, the measurements of are unequally spaced distributed, since the measurements along the line are projected perpendicularly onto frequency domain of the object along semicircular arc. Thus, the Fourier transform must be computed for every nonuniform frequency points ( or ). Although the result of direct nonuniform discrete Fourier transform (NDFT) is exact, the computation time required by the NDFT restricts its real application. To speed up, a fast NUFFT is employed to approximate NDFT in every iteration of reconstruction.

3.3.1. Conjugate Gradient Method

To solve the inverse problem (15) iteratively by CG, the gradient of the objective function must be computed as:

where represents the conjugate transpose. Since the absolute value function in norm and TV is nonsmooth function, we use approximation techniques to compute the corresponding gradient. For norm, the absolute value function is approximated with a smooth function by using the relation , where is a small positive smooth parameter. For TV, we use the following approximation strategy to avoid a zero denominator: , where is a small positive parameter. Therefore the gradient of norm and TV can be calculated:

3.3.2. NUFFT

NUFFT developed by Fessler and Sutton [51] is adopted. Consider the following 1D NUFFT case:

where is a vector of equally spaced samples of a signal and is a vector of nonuniform distributed frequencies. In matrix notation

where is nonuniform Fourier transform matrix. The NUFFT is implemented by two steps: firstly, project on an oversampled uniform Fourier basis by standard FFT

that is,

where . Secondly, approximate each by interpolating the using uniform samples

that is, . is the th row of interpolation matrix which makes use of neighboring uniform samples of for approximation of each nonuniform sample of . In [51], Fessler and Sutton designed a kind of min-max criterion to choose interpolation coefficients for every :

The analytical solution of (24) is

where denotes Hermitian transpose. Fessler and Sutton have shown that the overall complexity of such method is .

4. Simulation and Results

4.1. Simulation Parameters

In order to evaluate the performance of the proposed method for UDT reconstruction, we have performed a series of numerical experiments for the phantom in Figure 4(a). The phantom consists of ten ellipses which looks like the well-known Shepp-Logan “head phantom” for CT imaging. However, for UDT system, we have modified the gray levels to those used by [8, 9]. The gray levels represent the relative change in refractive index from the background value of 1.0; the maximum and minimum gray intensity are set to 1.0 and 0, respectively. The speed of sound of the background media is 1500 m/s. To evaluate our method, the scattered field was calculated based on FDT under Born approximation. Although the Born approximation imposes limitation on the dimension of the object for real application [4] and cannot distinguish the features of the object spaced less than [52, 53], it can provide a simple and direct method to reconstruct the structure of an object from the measurement of the scattered field. According to FDT, the Fourier transform of the scattered field measured on is proportional to Fourier transform of the object over an arc (2), while the Fourier transform of each ellipse has simple analytical expression; hence we can generate the scattered data through inverse Fourier transform. This procedure not only is fast but also allows the scattered date to be calculated for testing the reconstruction algorithms and experiments parameters such as pitch and number of elements [4, 8–10].

(a)

(b)

(c)

(d)

Figure 4: (a) Original image. (b) Image reconstructed using interpolation method. (c) Image reconstructed using broadband signal. (d) Image reconstructed using CS.

In numerical experiments, the imaging system utilizes a pair of parallel linear array probes [54, 55]. Referring to exiting commercial ultrasound linear transducer and ultrasound tomography system, the frequency of incident wave is set to 1.5 MHz, and the number of elements and pitch of the probes are 128 and , respectively, where is the wavelength of incident plane wave. The distance between the two probes is . The phantom Figure 4(a) is discretized on a Cartesian grid. According to the diffraction limitation [52, 53], the spatial sample step or the resolution of the system is set to . It is also necessary to point out that the numbers of iterations for iterative methods are all set to 8; we did not employ error threshold as the iteration criterion, because we want to compare the iterative algorithms after the same iteration numbers. According to the recommendation of Fessler and Sutton [51], the values of , are set to 6 and 2, respectively. The regularization parameters and are set to 0.01 and 0.001, respectively.

4.2. Results

Figure 4 shows the reconstructed images through different methods from simulated sparse-view data with no added noise. Figure 4(a) is the original phantom, Figure 4(b) is the reconstructed image using bilinear frequency interpolation, Figure 4(c) is the reconstructed image using the method proposed by Bronstein et al. [10] with broadband incident sound wave, and Figure 4(d) is the reconstructed image using the proposed method where the sparse transform basis is Haar wavelet. The number of view is 16, and the number of iterations is 8 for Figures 4(c) and 4(d). Due to sparse-view sampling, the Nyquist-Shannon sampling limitation cannot be satisfied. The reconstructed image of Figure 4(b) is severely blurred and distorted by interpolation error. The small scale features of original phantom cannot be recognized and larger ones are distorted by ring artifacts. Compared with the interpolation method, the two iterative methods can suppress the artifacts and reduce the noises and Gibbs effect as shown in Figures 4(c) and 4(d). In Figure 4(c), we adopted ten different frequency waves used by Bronstein in [10] to illuminate the object in each view. The image shows oversmoothing effect that the edges and details of the features are blurred. This may be caused by resolving the inverse problem of UDT under overdetermined framework; that is, there are more known than unknown .

Figure 4(d) shows the result of the proposed method. Most of the features can be clearly represented. The artifacts and oscillation noises are efficiently reduced. For the homogenous background (black region), there is no visible artifact and grayscale aberration, which accords with the characteristic of the original homogenous medium. For the narrow ring region , the outer and inner boundary can be preserved. The inner background is uniformly displayed except for parts of slightly blurred region. is set to be a low-contrast region (the gray intensity is 0.5, while it is 0.6 for background ), which can be used to evaluate the sensitivity of imaging methods. In Figure 4(d), the region can be clearly distinguished from the surrounding background with a sharp edge definition. Regions and are low and high refractive regions, with the gray intensity 0.3 and 0.75, respectively. The boundary of is preserved and the inner of is clean. The small region can be distinguished from surroundings. It is worth noting that the three small ellipses at region can be identified in Figure 4(d), while they are almost invisible in Figures 4(b) and 4(c). This implies the proposed method can improve the resolution efficiently in sparse-view situation.

Figure 5 shows the magnitude of error in the frequency domains. Since the maximum error of the interpolation method in Figure 5(a) is 4311, which is much bigger than the corresponding ones of the broadband method in Figure 5(c) (249.6) and the proposed method in Figure 5(d) (54.2), we also show the magnitude of error within the interval of for the interpolation method (Figure 5(b)). From Figure 5, we can find that the error of interpolation method around center frequency is significantly larger than the remaining two methods, which accords with Figure 4(b) that the details of the features are seriously distorted comparing to other two methods. In Figure 5(c), the maximum error around center frequency is lower than that of Figure 5(b); however, the error of frequencies other than center frequency area is generally higher than the corresponding ones in Figure 5(b). This means the noises in the reconstructed image through broadband method are higher than the reconstructed image through interpolation method. Figure 5(d) shows that the proposed method can markedly reduce the frequency domain error particularly in the low frequency components. This indicates that the reconstructed image can faithfully represent the feature details while efficiently suppressing the noises for sparse-view sampling data. Table 1 gives the relative mean square error (RMSE) in frequency domains, which is defined as

and are the distribution function of the original phantom and the reconstructed object in frequency domain, respectively. Table 1 also lists the structural similarity (SSIM) index for different methods. Compared with RMSE, the SSIM has proven to be consistent with human eye perception [56]. In this paper, the original object (Figure 4(a)) is the reference image for SSIM. Compared to interpolation method, the proposed method has relatively smaller RMSE and higher SSIM values which coincides with the description above.

Table 1: SSIM and RMSE for different methods.

Interpolation
Broadband signal
CS

SSIM
0.291
0.346
0.820

RMSE in frequency domain
0.765
0.352
0.255

(a)

(b)

(c)

(d)

Figure 5: Magnitude of error in frequency domain. (a) Interpolation method; (b) error value of interpolation method within ; (c) broadband method; (d) CS method.

In (15), the regularization parameters , determine the trade-off between the data consistency and the sparsity of the object. Furthermore, they can be used to adjust the relative weights of the different components in the cost function. Figure 6 shows the images reconstructed under different penalties with 16 views, 128 elements, and pitch of . Figure 6(a) is reconstructed with TV and norm regularization terms the same as Figure 4(d); the regularization parameters are set to 0.01 and 0.001, respectively. Figure 6(b) is reconstructed with solely TV regularization term; that is, , . Compared to Figure 6(a), although the TV penalty can preserve the edges of the object, the inner of the object is blurred. Figure 6(c) is reconstructed with solely norm regularization term; that is, , . The quality of the image, Figure 6(c), for inner region can match Figure 6(a), but the edges of the object are distorted. Figure 6(d) is directly reconstructed from the limited sample data without any penalties; that is, . Compared to Figure 6(a), not only are the qualities of the inner and edge seriously declined in Figure 6(d), but also the background is filled with noise. In practical, one can choose the regularization parameters , empirically based on the noise level and image contrast [57] or adaptively adjust the regularization parameters in the process of reconstruction [58, 59]. The primary target of this work is to accurately reconstruct the object in sparse-view, the selection of the optimal regularization parameters is not the focus of this paper, but it is worth further investigation.

(a)

(b)

(c)

(d)

Figure 6: Images reconstructed with different penalties: (a) , ; (b) , ; (c) , ; (d) .

To fully evaluate the proposed method, we also reconstruct the object with various numbers of views (Figure 7). In Figure 7(b), with 32 views, the smaller scale features (the three small ellipses) can be fully distinguished. Although the quality of reconstructed image for interpolation is also improved visually with the increase of views, the ring artifacts and oscillation noises cannot be eliminated even with 96 views.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 7: The reconstructed image through CS and interpolation methods with different views. (a), (b), (c), and (d) Reconstructed images by CS for 32 views, 48 views, 64 views, and 96 views, respectively; (e), (f), (g), and (h) reconstructed images by interpolation for 32 views, 48 views, 64 views, and 96 views, respectively.

5. Discussion

Not only can the sparse-view sampling scheme save the scan time of UDT but also it can reduce the complexity of the imaging device. In Section 4.2, we show the feasibility of image reconstruction for UDT in sparse-view situation based on CS framework. The two main reasons that CS is efficiently employed can be concluded as follows. Firstly, the acoustic indexes of the object are sparsely represented through orthogonal transform and finite difference transform; the former provides sparse representations for rapidly varied regions through wavelet transform, while the latter affords sparse transformation for piecewise smooth regions by TV. Secondly, the acoustic index coefficients with finite main components in the transform domain can be faithfully recovered through the iterative method based on CG and NUFFT.

To full evaluate our method, we will analyze the noise robustness of our method. Besides that, the computational complexity is also discussed in this part.

5.1. Robustness to Noise

In real situation, in addition to systematic errors such as misaligned transducers, the detected signals contain different kinds of noises such as thermal noise of transducer, electronic noise of amplifier. We added white Gauss noise in the received scatter data to test the robustness to noise of the proposed method. The signal-to-noise ratio (SNR) of simulated noisy data is 20 dB and 10 dB, respectively, where SNR is defined as

The pitch and the number of elements are and 192, respectively.

Figure 8 shows the reconstructed images with none, 20 dB, and 10 dB noise under 16 views and 32 views. The images reconstructed from the signal with an SNR of 20 dB have hardly any difference with the images reconstructed from noise-free signals. The homogenous medium background is not contaminated by noise. The edges of different structures can be distinguished.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: The robustness to noise of the proposed method. (a) 16 views, noise-free; (b) 16 views,  dB; (c) 16 views,  dB; (d) 32 views, noise-free; (e) 32 views,  dB; (f) 32 views,  dB.

However, when the noise increases to 10 dB, the quality of reconstructed image is affected by granulation noise. Although the main features of the object are still visible, the small details are seen to seriously deteriorate, especially for 16 views. How to increase the robustness to noise remains a further research topic.

5.2. Computational Complexity

In the proposed method, the computation time is mainly occupied by NUFFT. For an digital object, the complexity of NUFFT is [51], where is the oversampling constant and is the number of neighbors for interpolation. For normal reconstruction with the iteration number and the linear search times for each iteration, the total complexity of our method is estimated as . For comparison, the theoretical complexity of the frequency domain interpolation requires . Furthermore, the practical view number is generally more than to avoid aliasing. The sample points for every view are at least two times of . The complexity of filtered backpropagation (FBP) is about . For the tested image used in this paper, the frequency interpolation method requires about operations; spatial domain interpolation method (FBP) would require operations and our method requires about operations.

6. Conclusion

UDT is an important image modality and can afford functional exams in application. CS is one of the most exciting advances in signal theory which takes advantage of compressibility of the object to break Nyquist limitation and recover the major component in transform domain. The paper presents one CS framework for UDT image reconstruction. The numerical experiments show that the proposed method can improve image quality. The relative error is smaller than conventional interpolation method and broadband method. Combining norm with TV, not only is the edge of the object preserved, but also the contrast and resolution are improved.

In this paper our effort mainly focuses on integrating CS and UDT to develop practical framework for image reconstruction. Future work will be done on the evaluation of our study with in vivo data. Other important works include how to choose suitable transform basis and how to choose regularization parameters.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is partly supported by the National 973 Project of China (Grant no. 2011CB933103), the Project of the National 12th-Five Year Research Program of China (Grant no. 2012BAI13B02), the National Project of Developing Scientific Instrument of China (Grant no. 2013YQ160551), and the Self Innovation Research Fund of Huazhong University of Science and Technology (Grant no. 2013QN089).

References

J. F. Greenleaf and R. C. Bahn, “Clinical imaging with transmissive ultrasonic computerized tomography,” IEEE Transactions on Biomedical Engineering, vol. 28, no. 2, pp. 177–185, 1981.

P. Huthwaite, F. Simonetti, and N. Duric, “Combining time of flight and diffraction tomography for high resolution breast imaging: initial invivo results (l),” The Journal of the Acoustical Society of America, vol. 132, pp. 1249–1252, 2012.

R. K. Mueller, M. Kaveh, and G. Wade, “Reconstructive tomography and applications to ultrasonics,” Proceedings of the IEEE, vol. 67, no. 4, pp. 567–587, 1979.

A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, Society of Industrial and Applied Mathematics, 2001.

K. Iwata and R. Nagata, “Calculation of refractive index distribution from interferograms using the Born and Rytov's approximation,” Japanese Journal of Applied Physics, vol. 14, supplement 1, pp. 379–384, 1975.

F. Simonetti, L. Huang, and N. Duric, “On the spatial sampling of wave fields with circular ring apertures,” Journal of Applied Physics, vol. 101, no. 8, Article ID 083103, 2007.

F. Simonetti and L. Huang, “From beamforming to diffraction tomography,” Journal of Applied Physics, vol. 103, no. 10, Article ID 103110, 2008.

A. J. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrasonic Imaging, vol. 4, no. 4, pp. 336–350, 1982.

S. X. Pan and A. C. Kak, “A computational study of reconstruction algorithms for diffraction tomography: interpolation versus filtered backpropagation,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 31, no. 5, pp. 1262–1275, 1983.

M. M. Bronstein, A. M. Bronstein, M. Zibulevsky, and H. Azhari, “Reconstruction in diffraction ultrasound tomography using nonuniform FFT,” IEEE Transactions on Medical Imaging, vol. 21, no. 11, pp. 1395–1401, 2002.

S. J. LaRoque, E. Y. Sidky, and X. Pan, “Accurate image reconstruction from few-view and limited-angle data in diffraction tomography,” Journal of the Optical Society of America A: Optics and Image Science, and Vision, vol. 25, no. 7, pp. 1772–1782, 2008.

T. Li, Q. Liu, and J. Luo, “The Iterative Next-neighbor Regridding (INNG) algorithm combined with TV regularization used for reconstruction in diffraction tomography,” in Proceedings of the International Conference on Bioinformatics and Biomedical Technology (ICBBT '10), pp. 47–50, April 2010.

E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.

D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.

M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: the application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2008.

X. Pan and E. Y. Sidky, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Physics in Medicine and Biology, vol. 53, no. 17, pp. 4777–4807, 2008.

D. Wu, L. Li, and L. Zhang, “Feature constrained compressed sensing CT image reconstruction from incomplete data via robust principal component analysis of the database,” Physics in Medicine and Biology, vol. 58, pp. 4047–4070, 2013.

G. Chinn, P. D. Olcott, and C. S. Levin, “Sparse signal recovery methods for multiplexing PET detector readout,” IEEE Transactions on Medical Imaging, vol. 32, pp. 932–942, 2013.

Y. Zhang, Y. Wang, and C. Zhang, “Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction,” Ultrasonics, vol. 52, pp. 1046–1055, 2012.

G. Wang, “Guest editorial compressive sensing for biomedical imaging,” Magnetic Resonance in Medicine, vol. 30, pp. 1013–1016, 2011.

M. Schiffner and G. Schmitz, “Rapid measurement of ultrasound transducer fields in water employingcompressive sensing,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '10), pp. 1849–1852, October 2010.

M. Schiffner and G. Schmitz, “Fast pulse-echo ultrasound imaging employing compressive sensing,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '11), pp. 688–691, 2011.

M. Schiffner, T. Jansen, and G. Schmitz, “Compressed sensing for fast image acquisition in pulse-echo ultrasound,” Biomedical Engineering/Biomedizinische Technik, vol. 57, pp. 192–195, 2012.

N. Wagner, Y. Eldar, A. Feuer, and Z. Friedman, “Compressed beamforming applied to B-mode ultrasound imaging,” in Proceedings of the 9th IEEE International Symposium on Biomedical Imaging (ISBI '12), pp. 1080–1083, 2012.

N. Wagner, Y. Eldar, A. Feuer, and Z. Friedman, “Compressed beamforming with applications to ultrasound imaging,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP '12), pp. 3641–3644, 2012.

N. Wagner, Y. Eldar, and Z. Friedman, “Compressed beamforming in ultrasound imaging,” IEEE Transactions on Signal Processing, vol. 60, pp. 4643–4657, 2012.

M. Shen, Q. Zhang, and J. Yang, “A novel receive beamforming approach of ultrasound signals based on distributed compressed sensing,” in Proceedings of the IEEE International Instrumentation and Measurement Technology Conference (I2MTC '11), pp. 732–736, May 2011.

Q. Zhang, B. Li, and M. Shen, “A measurement-domain adaptive beamforming approach for ultrasound instrument based on distributed compressed sensing: initial development,” Ultrasonics, vol. 53, pp. 255–264, 2013.

A. Achim, B. Buxton, G. Tzagkarakis, and P. Tsakalides, “Compressive sensing for ultrasound RF echoes using α-stable distributions,” in Proceedings of the 32nd Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '10), pp. 4304–4307, Buenos Aires, Argentina, September 2010.

D. Friboulet, H. Liebgott, and R. Prost, “Compressive sensing for raw RF signals reconstruction in ultrasound,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '10), pp. 367–370, October 2010.

H. Liebgott, A. Basarab, D. Kouame, O. Bernard, and D. Friboulet, “Compressive sensing in medical ultrasound,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '12), pp. 1–6, 2012.

H. Liebgott, R. Prost, and D. Friboulet, “Pre-beamformed RF signal reconstruction in medical ultrasound using compressive sensing,” Ultrasonics, vol. 53, pp. 525–533, 2013.

C. Quinsac, A. Basarab, J. Girault, and D. Kouamé, “Compressed sensing of ultrasound images: sampling of spatial and frequency domains,” in Proceedings of the IEEE Workshop on Signal Processing Systems (SiPS '10), pp. 231–236, October 2010.

C. Quinsac, A. Basarab, D. Kouamé, and J. Grégoire, “3D Compressed sensing ultrasound imaging,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '10), pp. 363–366, October 2010.

C. Quinsac, A. Basarab, and D. Kouame, “Frequency domain compressive sampling for ultrasound imaging,” Advances in Acoustics and Vibration, vol. 12, pp. 1–16, 2012.

C. Quinsac, N. Dobigeon, A. Basarab, D. Kouame, and J. Tourneret, “Bayesian compressed sensing in ultrasound imaging,” in Proceedings of the 4th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP '11), pp. 101–104, December 2011.

C. Quinsac, F. de Vieilleville, A. Basarab, and D. Kouame, “Compressed sensing of ultrasound single-orthant analytical signals,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '11), pp. 1419–1422, 2011.

N. Dobigeon, A. Basarab, D. Kouame, and J. Tourneret, “Regularized bayesian compressed sensing in ultrasound imaging,” in Proceedings of the 20th European Signal Processing Conference (EUSIPCO '12), pp. 2600–2604, 2012.

S. M. S. Zobly and Y. M. Kakah, “Compressed sensing: doppler ultrasound signal recovery by using non-uniform sampling & random sampling,” in Proceedings of the 28th National Radio Science Conference (NRSC '11), pp. 1–9, April 2011.

S. Zobly and Y. Kadah, “Orthogonal matching pursuit & compressive sampling matching pursuit for doppler ultrasound signal reconstruction,” in Proceedings of the Cairo International Biomedical Engineering Conference (CIBEC '12), pp. 52–55, 2012.

J. Richy, D. Friboulet, A. Bernard, O. Bernard, and H. Liebgott, “Blood velocity estimation using compressive sensing,” IEEE Transactions on Medical Imaging, vol. 32, pp. 1979–1988, 2013.

L. He, T. C. Chang, S. Osher, T. Fang, and P. Speier, “MR image reconstruction by using the iterative refinement method and nonlinear inverse scale space methods,” Tech. Rep., UCLA CAM Report, 2006.

S. Mallat and E. Polytechique, A Wavelet Tour of Signal Processing: The Sparse Way, Academic Press, 3nd edition, 2008.

“Compressive sensing resources,” Rice University, Houston, Tex, USA, http://dsp.rice.edu/cs

E. J. Candès, “The restricted isometry property and its implications for compressed sensing,” Comptes Rendus Mathematique, vol. 346, no. 9-10, pp. 589–592, 2008.

E. J. Candès and T. Tao, “Decoding by linear programming,” IEEE Transactions on Information Theory, vol. 51, no. 12, pp. 4203–4215, 2005.

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Communications on Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207–1223, 2006.

J. Nocedal and S. Wright, A Wavelet Tour of Signal Processing, Springer, New York, NY, USA, 2nd edition, 2006.

M. Rudelson and R. Vershynin, “Geometric approach to error-correcting codes and reconstruction of signals,” International Mathematics Research Notices, no. 64, pp. 4019–4041, 2005.

E. J. Candès and T. Tao, “Near-optimal signal recovery from random projections: universal encoding strategies?” IEEE Transactions on Information Theory, vol. 52, no. 12, pp. 5406–5425, 2006.

J. A. Fessler and B. P. Sutton, “Nonuniform fast Fourier transforms using min-max interpolation,” IEEE Transactions on Signal Processing, vol. 51, no. 2, pp. 560–574, 2003.

E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Optics Communications, vol. 1, no. 4, pp. 153–156, 1969.

F. Simonetti, “Multiple scattering: the key to unravel the subwavelength world from the far-field pattern of a scattered wave,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, vol. 73, no. 3, Article ID 036619, 2006.

W. Jeong, T. Kim, D. Shin, S. Do, M. Singh, and V. Marmarelis, “Soft tissue differentiation using multiband signatures of high resolution ultrasonic transmission tomography,” IEEE Transactions on Medical Imaging, vol. 24, pp. 399–408, 2005.

G. Zografos, D. Koulocheri, P. Liakou et al., “Novel technology of multimodal ultrasound tomography detects breast lesions,” European Radiology, vol. 23, pp. 673–683, 2013.

Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.

T. F. Chan and P. Mulet, “On the convergence of the lagged diffusivity fixed point method in total variation image restoration,” SIAM Journal on Numerical Analysis, vol. 36, no. 2, pp. 354–367, 1999.

E. T. Hale, W. Yin, and Y. Zhang, “A fixed-point continuation method for l1-regularized minimization with applications to compressed sensing,” Tech. Rep. TR07-07, Rice University, 2007.

S. Ma, W. Yin, Y. Zhang, and A. Chakraborty, “An efficient algorithm for compressed MR imaging using total variation and wavelets,” in Proceedings of the 26th IEEE Conference on Computer Vision and Pattern Recognition (CVPR '08), pp. 135–142, IEEE, June 2008.