Efficient discrete cosine transform model–based algorithm for photoacoustic image reconstruction

Yan Zhang; Yuanyuan Wang; Chen Zhang

doi:10.1117/1.JBO.18.6.066008

3 June 2013 Efficient discrete cosine transform model–based algorithm for photoacoustic image reconstruction

Yan Zhang, Yuanyuan Wang, Chen Zhang

Author Affiliations +

Journal of Biomedical Optics, Vol. 18, Issue 6, 066008 (June 2013). https://doi.org/10.1117/1.JBO.18.6.066008

Abstract

The model-based algorithm is an effective reconstruction method for photoacoustic imaging (PAI). Compared with the analytical reconstruction algorithms, the model-based algorithm is able to provide a more accurate and high-resolution reconstructed image. However, the relatively heavy computational complexity and huge memory storage requirement often impose restrictions on its applications. We incorporate the discrete cosine transform (DCT) in PAI reconstruction and establish a new photoacoustic model. With this new model, an efficient algorithm is proposed for PAI reconstruction. Relatively significant DCT coefficients of the measured signals are used to reconstruct the image. As a result, the calculation can be saved. The theoretical computation complexity of the proposed algorithm is figured out and it is proved that the proposed method is efficient in calculation. The proposed algorithm is also verified through the numerical simulations and in vitro experiments. Compared with former developed model-based methods, the proposed algorithm is able to provide an equivalent reconstruction with the cost of much less time. From the theoretical analysis and the experiment results, it would be concluded that the model-based PAI reconstruction can be accelerated by using the proposed algorithm, so that the practical applicability of PAI may be enhanced.

1. Introduction

Photoacoustic imaging (PAI) is a novel noninvasive biomedical imaging technique.¹^–⁶ With its advantages of high contrast, high resolution, and functional imaging ability,⁴ it has great potential in many clinical applications such as structural and functional imaging,⁷^–⁹ tumor detection,¹⁰^,¹¹ ocular imaging,¹² vascular imaging,¹³^,¹⁴ and molecular imaging of biomarkers.¹⁵

The physical basis of PAI is the photoacoustic effect, and it refers to the generation of an ultrasound wave by the electromagnetic energy absorption.² Essentially, it is an energy transformation from the electromagnetic to the acoustic. In the process of PAI, the laser is used to irradiate the biomedical tissue to excite the ultrasound wave, and then the ultrasound signals will be received by the transducers. From the knowledge of measured signals, the photoacoustic image could be reconstructed, and this issue is regarded as an inverse problem. This problem in ultrasonic reflectivity imaging is first investigated for the spherical geometry,¹⁶ and then discussed in detail for plane, cylindrical, and spherical geometries.¹⁷ The exact solution and inversion formula are also obtained in these geometries.¹⁷ Actually, they provided a base for the later research of the photoacoustic image reconstruction. Not only that, more mathematics and formulas for the photoacoustic image reconstruction are deduced.¹⁸^,¹⁹ Based on these exact solutions and reconstruction mathematical formulas, many reconstruction algorithms are developed as follows. The exact reconstruction algorithms in the time and frequency domains have been developed in various geometries²⁰^–²³; Xu et al. proposed a filtered back-projection method²⁴^,²⁵; Zhang et al. proposed a deconvolution algorithm and adopted it into the fast PAI application.²⁶^,²⁷ The inverse spherical Radon transform is another classical reconstruction for PAI, and its fast solution method is proposed in Refs. 28 and 29. These analytical reconstruction algorithms are convenient to implement and have been widely used for PAI reconstruction.³

Besides these analytical algorithms, the model-based algorithm is another type of photoacoustic image reconstruction method³⁰ which has been quickly developed in recent years. Compared with the analytical algorithms, the model-based algorithms are able to provide a more accurate and high-quality reconstruction image.³⁰ Unlike solving the photoacoustic mathematic equations analytically, a model is built up to express the relationship between the measured ultrasound signals and the optical absorption map. The optimization methods such as the iterative method are used to inversely calculate the optical absorption image.³¹ Along with the model-based algorithm, many methods in other aspects can be involved in PAI reconstruction. The compressed sensing (CS)³² was adopted into PAI reconstruction,³³^,³⁴ and a nonlinear conjugate gradient descent algorithm was used to implement the CS-based photoacoustic imaging.³³ It is shown from the results that the image can be reconstructed with fewer measurements and the data acquisition can be accelerated. Additionally, the wavelet transform is also applied into PAI reconstruction and a new reconstruction strategy was proposed.³⁵ By recovering the image with a sparse representation in the wavelet domain, the performance of the reconstruction can possibly be improved. Some parameters in the aspect of image processing are also employed to assist the image reconstruction. As the total variation (TV) coefficient utilized in company with the gradient descent method in PAI reconstruction, the TV-GD method³⁶ is proposed and proved to be useful for PAI reconstruction in the sparse-view sampling condition. From the above discussions, the model-based algorithm is widely accepted for PAI reconstruction for its high-quality performance. However, the relatively heavy computational complexity and huge memory storage requirement are difficulties in implementing these algorithms. Though in some of the developed algorithms, the number of sampling points can be reduced, the model-based reconstruction still faces the problem of heavy computation.

In this paper, we make use of the discrete cosine transform (DCT) and establish a new model for PAI image reconstruction. In this new model, the relationship between the DCT coefficients of the measured signals and the optical absorption image is presented. We calculate the DCT coefficients both from the measured signals and from an assumed optical absorption image using the proposed model. The difference of these two DCT coefficients is minimized to accomplish the image reconstruction. For efficiency, the relatively significant DCT coefficients are used in the reconstruction. In this way, the calculation can be remarkably diminished. Under an appropriate threshold, the amount of the DCT coefficients for reconstruction can be reduced. However, the remaining part still carries the sufficient information for the image reconstruction. Thus it is noted that this modification would not affect the quality of the reconstruction, but the calculation would be reduced. The theoretical calculation complexity is figured out, and the essential advantage of the proposed algorithm is revealed. Through numerical simulations and in vitro experiments, this new algorithm is verified and compared with other developed algorithms. It is demonstrated in the experimental results that the reconstruction can be accelerated effectively and the quality of the reconstruction can be kept at a high level as well.

The paper is organized as follows. Section 2 introduces the DCT-based photoacoustic model and the procedures of the DCT model-based image reconstruction. Section 3 presents the calculation of theoretical algorithm complexity. The numerical simulations are demonstrated in Sec. 4, where the reconstruction results are discussed. Section 5 introduces the experimental system, and the in vitro experiments on this system are shown. Finally, Sec. 6 is the conclusion.

2. Method

In this paper, two-dimensional (2-D) PAI is concerned in our simulations and experiments. According to the physical principle of the photoacoustic effect, assuming the imaging tissue is spatially uniformly illuminated by the laser, the relationship between the photoacoustic signals and the laser energy deposition can be derived as⁴

Eq. (1)

\nabla^{2} p (r, t) - \frac{1}{c^{2}} \frac{\partial^{2} p (r, t)}{\partial t^{2}} = - \frac{β}{C_{p}} A (r) \frac{\partial I (t)}{\partial t},

where

p (r, t)

is the acoustic pressure at the position

r

and at the time

t

,

c

is the sound speed,

C_{p}

is the specific heat,

β

is the isobaric expansion coefficient,

I (t)

is the temporal profile of the laser pulse, and

A (r)

is the spatial distribution of the laser absorption.

Equation (1) can be solved by using Green’s function, and the photoacoustic signal is deduced as⁴

Eq. (2)

p (r_{0}, t) = \frac{β}{4 π C_{p}} \frac{\partial}{\partial t} ∯_{| r^{'} - r_{0} | = c t} \frac{A (r^{'})}{t} d^{2} r^{'},

where

r_{0}

is the position of the ultrasound transducer.

In PAI experiments, a scanning transducer is used to receive the photoacoustic signals at different positions around the imaging sample, and the image reconstruction is regarded as an inverse problem. In the model-based image reconstruction, a model is typically established to connect the measured signals with the optical deposition image. The signal can be calculated based on the optical deposition in the model, and then the image can be reconstructed by minimizing the difference between calculated signals and the real measured signals. In this way, many optimization methods can be used, and various model-based reconstruction algorithms have been developed.

In is paper, involving the discrete cosine transform (DCT), we build up a novel model for the image reconstruction.

For a succinct expression, a new variable $g$ is defined as

Eq. (3)

g (r_{0}, t) = \frac{4 π C_{p} t}{β} \cdot \int_{0}^{t} p (r_{0}, t) d t .

Then Eq. (2) can be transformed to

Eq. (4)

g (r_{0}, t) = ∯_{| r^{'} - r_{0} | = c t} A (r^{'}) d r^{'} .

The reconstructed optical deposition image and the photoacoustic signals are processed discretely. The scanning trajectory is discretized into $M$ points so that the photoacoustic signals in $M$ points are used for the reconstruction. The image $A$ is discretized into a matrix in the size of $N_{x} \times N_{y}$ , and then $A$ is reshaped into a column vector $A^{'}$ in the size of $N_{x} \times N_{y} \times 1$ for calculation. Assuming that the measured vector $g$ is in the size of $T \times 1$ , the relationship between the $g$ and the $A^{'}$ can be expressed in the form of matrix multiplication as

Eq. (5)

g (i) = W (i) \cdot A^{'}, i = 1, 2, 3, \dots, M,

where

M

is the total number of the detection points and

i

means the

i

’th detection point.

In Eq. (5), the size of the matrix $W$ is $T \times N_{x} N_{y}$ , and it can be calculated based on the scanning geometry in PAI. At the $i$ ’th sampling point, the matrix $W (i)$ can be calculated as

Eq. (6)

\begin{array}{l} W_{j k} (i) = {\begin{cases} 1 - | \frac{t_{j}}{Δ t} - \frac{| r_{k} - r_{i} |}{c Δ t} | & when | \frac{t_{j}}{Δ t} - \frac{| r_{k} - r_{i} |}{c Δ t} | < 1 \\ 0 & else \end{cases}, & i = 1, 2, 3, \dots, M, \end{array}

where

c

is the sound speed,

Δ t

is the temporal discretization interval,

t_{j}

is the

j

’th temporal measurement at the

i

’th detection point,

r_{k}

is the location vector of the

k

’th point in

A^{'}

, and

r_{i}

is the location vector of the

i

’th detection point. In the matrix

W

, the locations of the points in the image are recorded, and the elements are valued between 0 and 1, which denotes the contribution to the signal at a certain time. Note that no smoothness will be involved in the reconstructed image because of this discretization.

In our model, the variable $g$ is processed with the DCT and the obtained DCT coefficients are used for reconstructing the image. The arithmetic formula of the DCT is

Eq. (7)

G_{k} = w_{k} \sum_{t = 1}^{T} g_{k} \cos \frac{π (2 t - 1) (k - 1)}{2 T},

where

Eq. (8)

w_{k} = {\begin{matrix} 1 / \sqrt{T} & k = 1 \\ \sqrt{2 / T} & 2 \leq k \leq T \end{matrix} .

The DCT calculation can also be processed in the form of matrix multiplication as

Eq. (9)

G (i) = D g (i),

where

D

is the DCT matrix, and its size is

T \times T

. According to Eqs. (7) and (8), the matrix

D

can be calculated as

Eq. (10)

D_{st} = {\begin{matrix} \frac{1}{\sqrt{T}} \cos [\frac{π (2 t - 1) (s - 1)}{2 T}] & s = 1 \\ \sqrt{\frac{2}{T}} \cos [\frac{π (2 t - 1) (s - 1)}{2 T}] & 2 \leq s \leq T \end{matrix},

where

s

and

t

are the row and column indexes in the DCT matrix

D

, respectively.

In this way, the relationship between the DCT coefficients and the photoacoustic image can be expressed as

Eq. (11)

G (i) = D W (i) \cdot A^{'}, i = 1, 2, 3, \dots, M .

Here, a threshold is set and the significant DCT coefficients with a relatively big value will be picked out. When an appropriate threshold is set, these selected coefficients contain nearly most of the information in the measured signals, and they are sufficient for the image reconstruction. After the selection of the coefficient, the amount of the data is supposed to be reduced, and this procedure is helpful to save the calculation.

Supposing the threshold is TH, whose value will be discussed in the following section, the selection index set can be obtained by

Eq. (12)

{k | k = 1, 2, 3, \dots, T | G_{k} | > TH} .

The DCT coefficient $G$ is reduced to $G^{TH}$ after the selection, and it can be expressed as

Eq. (13)

{G^{TH} | k = 1, 2, 3, \dots, T | G_{k} | > TH} .

Assuming $F (i) = D W (i)$ , the measurement matrix $F$ is modified correspondingly by removing the rows whose indices are not in the set [Eq. (12)]. Then the reduced matrix relative to the threshold, which is denoted by $F^{TH}$ , can be obtained.

The model is established to connect the reduced DCT coefficients with the optical absorption image. Based on the knowledge of the reduced DCT coefficient $G^{TH}$ and the corresponding reduced matrix $F^{TH}$ , the photoacoustic image is reconstructed by iterations. According to the gradient descent method, the iteration can be processed as

Eq. (14)

Δ A^{'} = - \frac{{[F^{TH} (i)]}^{T}}{‖ {[F^{TH} (i)]}^{T} ‖} [F^{TH} (i) \cdot A^{'} - G^{TH} (i)] i = 1, 2, 3, \dots, M,

where

{(*)}^{T}

means the matrix transposition.

According to Eq. (14), the reconstruction vector $A^{'}$ is updated at $M$ different sampling points. Then the iteration is repeated unless the exiting condition is met. Generally, the exiting criterion is set so that the error is under a certain level or the iteration number is more than a prespecified number.

Finally in this section, the proposed DCT model-based algorithm is summarized in four steps as follows:

1. Initialization: Initialize the reconstructed image to be a zero matrix and reshape it into a column vector. Note that the initial matrix can be also set as the analytically reconstructed image, and the remaining artifacts should be taken under consideration. Then apply Eqs. (3) and (9) to calculate $G$ and apply Eqs. (6) and (10) to calculate $F (i)$ .
2. Reduction: Set the threshold TH, and calculate the selection index set by Eq. (12). According to the selection index set, remove the relatively insignificant elements in the DCT coefficients vectors and the corresponding rows in matrix $F (i)$ .
3. Iteration: For $i = 1, 2, 3, \dots, M$ , implement the iteration at each sampling point as shown in Eq. (14).
4. Update or termination: Terminate the iteration when the exiting criterion is met, or update the vector $A^{'}$ and return to step 3 to resume the iteration. Usually the exiting criterion is set so that the error is under a certain level or the iteration number is more than a prespecified number.

3. Theoretical Calculation Complexity

The main motivation for involving the DCT into PAI reconstruction is its effectiveness in saving the computation. In this section, we focus on the calculation complexity in this model-based reconstruction. The theoretical complexity is calculated according to the amount of multiplication in the matrix processing.

In calculating the measurement matrixes, the complexity is determined by the Eqs. (6), (9), and (10). Assuming the reconstructed image is in the size of $N_{x} \times N_{y}$ , the number of sampling angles is $M$ , and the length of the discretized signal is $T$ . In calculating $W$ and the DCT matrix $D$ , by implementing the same procedures on every element in the matrix, the relative complexities are $O (N_{x} N_{y} M T)$ and $O (T^{2})$ . In the multiplication of $W$ by $D$ , it requires $O (N_{x} N_{y} M T^{2})$ operations. From the above discussions, in the model establishment, the calculation of the matrixes leads to a calculation complexity of $O [(N_{x} N_{y} M T) + (T^{2}) + (N_{x} N_{y} M T^{2})]$ . Compared with the iterative reconstruction (IR) algorithm,³¹ the second and the third terms are extra. In other words, the extra calculation is necessary to establish the DCT model for the reconstruction. However, this extra calculation is much less than the saved calculation in the reconstruction as per the following discussion.

In the reconstruction, the naive length of the DCT coefficients is $T$ . By setting a threshold TH, the length of DCT vector can be reduced to $α T$ . Consequently, the amount of the data to process is reduced to $α M T$ . In the first step, the calculation of vector $g$ and the DCT processing require $2 O (T^{3})$ operations. Then in one iteration step, the complexity is $O [N_{x} N_{y} {(α T)}^{2} + {(N_{x} N_{y})}^{2} α M T]$ . Supposing the iteration is repeated $N_{iter}$ times, the resulting calculation complexity is $O [2 T^{3} + N_{iter} N_{x} N_{y} M {(α T)}^{2} + N_{iter} {(N_{x} N_{y})}^{2} α MT]$ . In contrast, with the unreduced data, the calculation complexity of the developed IR algorithm will be $O [2 T^{3} + N_{iter} N_{x} N_{y} M {(T)}^{2} + N_{iter} {(N_{x} N_{y})}^{2} M T]$ .

The first two terms are much smaller than the last one, so the entire complexity is mainly determined by the last term. For instance, if a proper threshold TH is set, the resulting $α$ value will be 0.15. Under a consistent reconstruction condition, when the other parameters are the same, the calculation is determined by the reduction ratio $α$ . As the third term in the complexity shows, the calculation in the reconstruction can be reduced to about 15% of the calculation of the IR algorithm.

In this section, the theoretical complexity of calculation is figured out. From the above theoretical analysis, it may be proved that the proposed DCT model-based algorithm is theoretically advantagious in saving calculations.

4. Simulation

Numerical simulations are used to verify the proposed DCT model-based PAI reconstruction. It is necessary to point out that the simulation and image reconstruction are all performed in 2-D. The image reconstructed by the proposed DCT-based algorithm and the former developed ones are respectively presented in this section. Additionally, the qualitative and quantitative analyses are presented to demonstrate the advantages of the proposed algorithm.

The simulation experiments are performed using Matlab (version 7.8.0.0347) on a PC with 2.0 GHz Quad-CPU and 4.0 GB memory.

In this simulation, the given optical absorption distribution image is shown in Fig. 1. The size of the reconstruction area is $60 \times 60 {mm}^{2}$ , and the radius of the scanning circle is 60 mm. The sound speed is consistently $1500 m / s$ in the simulation.

Fig. 1

The given optical absorption distribution image.

Concerning the proposed algorithm, the quality of the reconstructed image, the selection of the threshold, the few-view sampling condition, and the robustness will be discussed in the following sections.

4.1.

Selection of the Threshold

As mentioned before, a threshold is utilized to establish the model in the proposed method. This threshold is expected to be assigned an appropriate value, which is useful to make an effective compression to the data and preserve the necessary information for the reconstruction.

In this simulation, Eq. (2) is employed to calculate the photoacoustic signals based on the optical absorption map. We calculate the signals at 45 different positions, which are uniformly distributed around the imaging sample. By assigning the thresholds 0.001, 0.01, 0.02, 0.05, different models are built up and the images are reconstructed respectively. For a consistent reconstruction condition, the numbers of the iterations are all set to 20. The size of the reconstructed image is $150 \times 150 pixels$ . The simulated results are shown in Fig. 2.

Fig. 2

The simulation results of the DCT model-based reconstruction algorithm with different thresholds: (a) 0.001, (b) 0.01, (c) 0.02, and (d) 0.05.

In Fig. 2, it is clearly shown that the reconstructed images with a threshold of 0.001 and 0.01 are highly consistent with the original optical deposition image. In the proposed model, a smaller threshold means that more DCT coefficients will be preserved for the reconstruction, and it leads to a better reconstruction. As the threshold increases, the quality of the reconstruction declines. When the threshold increases to 0.05, the reconstructed image is severely blurred, and the optical absorption lines are much thicker than those in the original image. Note that the reconstructed image with a high threshold is not refined with more iteration. The constraint of the algorithm is the threshold rather than the iteration numbers. The high-quality image cannot be reconstructed from the over-reduced sampling data because too much information for the reconstruction is removed.

Besides the quality of reconstruction, we also take the efficiency into consideration to determine the value of the threshold. With different thresholds, we recorded the amounts of the data after the reduction and the times cost in the reconstructions. Additionally, with Fig. 1 as a standard, we calculate the peak signal-to-noise ratios (PSNRs) of four reconstructed images. All of these contents are listed in Table 1.

Table 1

Comparison between the DCT model-based reconstructions with different thresholds.

Threshold	Amount of data in reconstruction	Time cost (s)	PSNR (dB)
0.001	2713	480.6	26.5
0.01	981	189.5	24.0
0.02	636	130.1	23.0
0.05	356	82.0	20.0

In Table 1, the threshold is absolute and determined manually and it is not the factor $α$ . As shown in Table 1, with a smaller-valued threshold, the quality of the reconstructed image is higher, however, the relative amount of the data in reconstruction is larger and more time is needed as a result. It is a compromise of the reconstruction quality and the calculation saving. Considering the visual effect in Fig. 2(a) and 2(b), the reconstructed images are quite similar, and the image details are both refined. In Fig. 2(c), the cross point is slightly coarse with a threshold of 0.02, and the reduction of the amount of data and the time cost in reconstruction are no longer remarkable. Hence, 0.01 seems to be an appropriate threshold in this method. As the threshold is determined, we compare the proposed algorithm with the other developed algorithms, the FBP and the TV-GD methods. The comparison is shown in Fig. 3.

Fig. 3

(a) The original optical absorption image and the reconstructed images with (b) the FBP method, (c) the TV-GD method, and (d) the DCT model-based algorithm.

As in Ref. 24, when the number of sampling views is not sufficient, the FBP method suffers from the problem of artifacts and would not be able to yield an accurate image reconstruction. It is also clearly shown in Fig. 3(a) that the FBP reconstructed image is severely ruined by the remaining artifacts, and the optical absorbers are ambiguous. It could be deduced that the FBP method is not applicative in this sparse-sampling view condition. In contrast, much better reconstructions are obtained by using these two model-based algorithms. Three lines are clearly shown and the contrasts are much higher. The superiority of the model-based PAI reconstruction is revealed.

In order to show the details of reconstructed images clearly, we draw a column of pixels (the 90th column) from each image reconstructed by the FBP method, the TV-GD method, and the DCT model-based algorithm, respectively. The comparisons of pixel profiles are displayed in Fig. 4.

Fig. 4

The grayscale profiles of the original image and the reconstructed images by the FBP method, the TV-GD method, and the DCT model-based algorithm, respectively.

In Fig. 4, the solid, the dashed, the dotted, and the star lines represent the pixel profiles of the original, the FBP, the TV-GD, and the DCT model-based images, respectively.

It is shown in Fig. 4 that the profile of the FBP image significantly differs from the original profile. Because of the severe artifacts in the FBP image, there exist huge vibrations in the pixel profile. The profiles of the TV-GD and the DCT model-based images match with the original one much more precisely. Though the differences still exist, these two iterative methods are proved to possess equivalent reconstruction qualities.

The main advantage of this proposed algorithm is the efficiency. To reconstruct Fig. 3(c), it takes up to 682.4 s to accomplish the reconstruction with the TV-GD algorithm. It is much more than the 189.5 s cost in reconstructing Fig. 3(d) with the DCT-based algorithm. The involved DCT is effective in reducing the data amount.

As reported in Ref. 33, it takes 10 to15 h to reconstruct an image in size of $256 \times 256 pixels$ with the Matlab platform on a PC. For a consistent simulation condition, the proposed algorithm is utilized to reconstruct an image in the same size with the equivalent platform. As a result, it takes only 294 s to accomplish the image reconstruction. It is rational that the proposed DCT model-based algorithm may be more efficient in the computation.

Compared with the wavelet-based algorithm, the DCT-based algorithm is implemented by the iteration instead of the pseudo-inverse matrix processing. The procedures are simpler structures and the storage can be saved. Besides, the algorithm is convenient to adjust according to the scanning geometry. Based on the signals’ sparseness in the wavelet and the DCT domains, the employment of sparseness property can assist the algorithms to obtain high-quality reconstructions.

From the above discussions, it could be concluded that the proposed DCT model-based algorithm is an efficient and practical reconstruction algorithm.

4.2.

Few-View Sampling Condition

In this part, we focus on the obliged number of sampling points for the DCT model-based image reconstruction. It is noted that reducing the sampling points is meaningful in reducing PAI system cost and accelerating the data acquisition. Setting an appropriate sampling number is practically useful in PAI applications. In order to compare the performance of the proposed algorithm with different amounts of sampling views, five different situations are simulated. In these five cases, the sampling angles all cover 360 deg around the imaging object, the angular step size is uniform and set to 4 deg, 6 deg, 8 deg, 12 deg, and 18 deg, respectively. To avoid the noise infection, the simulated signal is noise free. The images are reconstructed based on these simulated data with different numbers of sampling points, and the results are shown in Fig. 5.

Fig. 5

(a) The original optical absorption image and the DCT model-based reconstructed images with simulated data in (b) 20, (c) 30, (d) 45, (e) 60, and (f) 90 sampling points.

It is shown in Fig. 5(d)–5(f) that the reconstructed images have no essential differences when the number of sampling angles is more than 45. To reconstruct an imaging sample with an uncomplicated structure as shown in Fig. 1, the time cost in processing the extra sampled signal is not necessary, and the extra calculation time does not lead to a better reconstruction. However, when the sampling angles are fewer, not enough information will be provided for the image reconstruction. In Fig. 5(b), the reconstructed lines are severely curved, and the reconstructed image is distorted. The quality of the reconstruction is downgraded as a result. In Fig. 5(c), the image is still slightly affeced by the noise, and the noise suppression is not as effective as those shown in Fig. 5(d) and 5(e).

The photoacoustic data with less than 45 sampling points may not be able to provide sufficient information for the reconstruction. The reconstructed image is easily affected and degraded as a result. The data with more sampling points are not necessary, and the reconstruction may not be apparently improved by involving more signals.

In the simulation above presented, the imaging sample is simply structured. To demonstrate the availability of the proposed algorithm in reconstructing a complex structured imaging object, we employ a real MRI image as the optical absorption image and implement the numerical simulation. The TV-GD method and the proposed DCT model-based method are both used for the image reconstruction. Note that the number of sampling points is also set as 45 and the reconstruction results are shown in Fig. 6.

Fig. 6

(a) The original optical absorption image and the reconstructed images from 45-view simulated data by (b) the TV-GD method, and (c) the DCT model-based method.

It is shown that the DCT model-based reconstructed image highly agrees with the original image. The three dark areas are reconstructed with consistent boundaries, and the high-valued skull is also accurately reconstructed. Besides, the proposed algorithm has an equivalent reconstruction quality with the TV-GD method, while the computation is nearly one third of that.

From the above discussions, the appropriate number of sampling points can be set as 45 in the simulation. The image reconstructed with 45 sampled data remains at quite a high-quality level, and the DCT model-based algorithm may be proved to be practical.

4.3.

Robustness

In the progress of the signal acquisition, the detected photoacoustic signals are usually affected by the thermal noises from the ultrasound transducer and the system electronics. The reconstruction algorithm is unavoidably interrupted by these noises, and an algorithm with great robustness is expected to maintain stability in the noisy circumstance. In real applications, the robustness is also an important consideration factor besides the reconstruction quality and the calculation complexity.

To investigate the robustness of the proposed algorithm, white Gaussian noises of different powers are added to the simulated signals. The signal-to-noise ratio (SNR) of the noisy signal is 10, 5, 3, and 0 dB, respectively. With these noisy signals, the DCT model-based algorithm is used to reconstruct the images, and the results are shown in Fig. 7.

Fig. 7

The DCT model-based reconstruction results from noisy data with the SNR of (a) 0 dB, (b) 3 dB, (c) 5 dB, and (d) 10 dB.

It is shown in Figs. 7(b)–7(d) that the images reconstructed from signals with the SNR greater than 3 dB have hardly any difference with the image reconstructed from noise-free signals as shown in Fig. 3(d). The noise suppression is effective, and the details in the reconstructed images are well preserved.

When the noise gets stronger, for example the SNR deteriorates to 0 dB as shown in Fig. 7(a), the image is ruined by the noise, and much impulse noise appears in the background of the reconstructed image. Even the line is segmented by the noises and the reconstructed image is not consistent with the original optical deposition image.

In this section, we implement the reconstructions from signals with different powers of the noise. In this way, the robustness of the proposed algorithm is investigated. As a result, the algorithm is applicative unless the SNR of the signal is worse than 3 dB. Based on the signal with the SNR better than 3 dB, the algorithm is stable and it can avoid the noise influence.

5. In Vitro Experiments

To validate the proposed DCT model-based PAI reconstruction algorithm, the in vitro experiments are conducted based on PAI system illustrated in Fig. 8.

Fig. 8

The scheme of PAI system.

In this system, an $Nd ∶ YAG$ laser generator (Continum, Surelite I) is employed. The wavelength of the emitting laser is 532 nm, the temporal duration of the laser pulse is 5 to 7 ns, the repetition rate is 10 Hz, and the laser energy density in the illumination area is about $6.47 mJ \cdot {cm}^{- 2}$ , which is lower than the laser radiation safety standard of $20 mJ \cdot {cm}^{- 2}$ (Ref. 2). An unfocused water-immersion ultrasound transducer (Panametric, V383-SU) is controlled by the stepping motor (GCD-0301M), and it is used to receive the photoacoustic signals. The diameter of the active element is 9.525 mm, the central frequency is 3.5 MHz, and the bandwidth is 1.12 MHz. The analog signals are received by the pulse receiver (Panametric, 5900PR) and sampled into digital signals by the oscilloscope (Agilent, 54622D). The sampling frequency is set to 16.67 MHz. Then the sampled data is transported to the PC through the general purpose interface bus (GPIB).

In this experiment, the imaging sample is a cylinder made of gelatin in which two black rubber ropes are embedded as the optical absorbers. The radius of the gelatin cylinder is 25 mm and the length of the embedded rubber ropes are 20 and 12 mm. The photograph of the gelatin sample is shown in Fig. 9.

Fig. 9

The photograph of the imaging sample used in the in vitro experiment.

In the experiment, the scanning radius is 38 mm. The angular step of sampling is set as 8 deg, and the signals at 45 positions can be received. The transducer is placed in the same plane as the optical absorbers, and it tends to receive in-plane only. As the cross-sectional image is mainly determined by the measured data in the certain plane, a set of circular measurement data in the same plane is sufficient for the image reconstruction.²⁴ Then the images are respectively reconstructed by the FBP, the IR, the TV-GD, and the DCT model-based algorithms. The reconstruction results are shown in Fig. 10.

Fig. 10

The reconstructed images by (a) the FBP, (b) the IR, (c) the TV-GD, and (d) the DCT model-based algorithms in in vitro experiments.

It is shown in Fig. 10 that the reconstructions of iterative algorithms are much better than the FBP reconstruction. In the FBP image, the optical absorbers are seriously deteriorated because of the severe artifacts and the noises. By contrast, in the image reconstructed by the proposed DCT model-based algorithm, the optical absorbers are clearly shown. The edges of the objects are distinct, and they are quite easy to recognize. In the aspect of noise suppression in the background, the proposed algorithm is still compared to the other methods. Fewer artifacts and noises may contribute to better application in clinical analysis. The high performance of the DCT model-based reconstruction algorithm is revealed in the in vitro experiments.

Similar to the simulation discussions, we recorded the times cost in the reconstructions with three iterative algorithms as well. The times are 119.98, 120.07, and 42.34 s, respectively, by using the IR, the TV-GD, and the DCT model-based algorithms. The reconstruction quality is highly acceptable, while the time cost is less than the IR and the TV-GD methods. From the in vitro experimental results, it is proved that the proposed algorithm is practical in application and efficient in computation.

6. Conclusion

In this paper, we have adopted the DCT and established a new model for PAI. In this new model, an appropriate threshold is determined, and an acceleration strategy is formulated. After implementing the DCT to the measured signals, the minor coefficients whose absolute value is less than the threshold are removed. The measurement matrixes calculated according to the scanning geometry can also be reduced correspondingly. As the relationship between the reduced DCT coefficients and the optical deposition image deduced, the model is built up. Based on the proposed model, the iterative method can be used to accomplish the image reconstruction. The simulations and the in vitro experiments are carried out and the results are presented. It is revealed that the quality of the reconstructed images is equivalent to the former developed model-based algorithms. Additionally, the calculation complexity is figured out, and the advantage of the efficiency is demonstrated. From the times cost comparison among the DCT model-based algorithm and the former developed algorithms, it is proved that the proposed algorithm is more efficient in calculation. In conclusion, by employing the DCT model-based algorithm, PAI reconstruction can be effectively accelerated, and it may be useful to enhance the practical applicability of PAI reconstruction.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 61271071 and 11228411), the National Key Technology R&D Program of China (No. 2012BAI13B02), and Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20110071110017).

References

1.

L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron., 14 (1), 171 –179 (2008). http://dx.doi.org/10.1109/JSTQE.2007.913398 IJSQEN 1077-260X Google Scholar

2.

M. XuL. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). http://dx.doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

3.

C. LiL. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol., 54 (19), R59 –R97 (2009). http://dx.doi.org/10.1088/0031-9155/54/19/R01 PHMBA7 0031-9155 Google Scholar

4.

L. V. Wang, “Prospects of photoacoustic tomography,” Med. Phys., 35 (12), 5758 –5767 (2008). http://dx.doi.org/10.1118/1.3013698 MPHYA6 0094-2405 Google Scholar

5.

R. A. KrugerD. R. ReineckeG. A. Kruger, “Thermoacoustic computed tomography—technical considerations,” Med. Phys., 26 (9), 1832 –1837 (1999). http://dx.doi.org/10.1118/1.598688 MPHYA6 0094-2405 Google Scholar

6.

R. A. Krugeret al., “Photoacoustic ultrasound (PAUS)-reconstruction tomography,” Med. Phys., 22 (10), 1605 –1609 (1995). http://dx.doi.org/10.1118/1.597429 MPHYA6 0094-2405 Google Scholar

7.

X. D. Wanget al., “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol., 21 (7), 803 –806 (2003). http://dx.doi.org/10.1038/nbt839 NABIF9 1087-0156 Google Scholar

8.

H. F. Zhanget al., “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechnol., 24 (7), 848 –851 (2006). http://dx.doi.org/10.1038/nbt1220 NABIF9 1087-0156 Google Scholar

9.

L. Liet al., “Photoacoustic imaging of lacZ gene expression in vivo,” J. Biomed. Opt., 12 (2), 020504 (2007). http://dx.doi.org/10.1117/1.2717531 JBOPFO 1083-3668 Google Scholar

10.

B. Guoet al., “Multifrequency microwave-induced thermal acoustic imaging for breast cancer detection,” IEEE Trans. Biomed. Eng., 54 (11), 2000 –2010 (2007). http://dx.doi.org/10.1109/TBME.2007.895108 ITUCER 0885-3010 Google Scholar

11.

M. Pramaniket al., “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys., 35 (6), 2218 –2223 (2008). http://dx.doi.org/10.1118/1.2911157 MPHYA6 0094-2405 Google Scholar

12.

A. de la Zerdaet al., “Photoacoustic ocular imaging,” Opt. Lett., 35 (3), 270 –272 (2010). http://dx.doi.org/10.1364/OL.35.000270 OPLEDP 0146-9592 Google Scholar

13.

E. Z. Zhanget al., “In vivo high-resolution 3D photoacoustic imaging of superficial vascular anatomy,” Phys. Med. Biol., 54 (4), 1035 –1046 (2009). http://dx.doi.org/10.1088/0031-9155/54/4/014 PHMBA7 0031-9155 Google Scholar

14.

J. J. Niederhauseret al., “Combined ultrasound and optoacoustic system for real-time high-contrast vascular imaging in vivo,” IEEE Trans. Med. Imag., 24 (4), 436 –440 (2005). http://dx.doi.org/10.1109/TMI.2004.843199 ITMID4 0278-0062 Google Scholar

15.

A. De La Zerdaet al., “Carbon nanotubes as photoacoustic molecular imaging agents in living mice,” Nat. Nanotechnol., 3 (9), 557 –562 (2008). http://dx.doi.org/10.1038/nnano.2008.231 1748-3387 Google Scholar

16.

S. J. NortonM. Linzer, “Ultrasonic reflectivity imaging in three dimensions: reconstruction with spherical transducer arrays,” Ultrasonic Imag., 1 (3), 210 –231 (1979). http://dx.doi.org/10.1016/0161-7346(79)90017-8 ULIMD4 0161-7346 Google Scholar

17.

S. J. NortonM. Linzer, “Ultrasonic reflectivity imaging in three dimensions: exact inverse scattering solutions for plane, cylindrical, and spherical apertures,” IEEE Trans. Biomed. Eng., BME-28 (2), 202 –220 (1981). http://dx.doi.org/10.1109/TBME.1981.324791 IEBEAX 0018-9294 Google Scholar

18.

D. FinchS. K. PatchRakesh, “Determining a function from its mean values over a family of spheres,” SIAM J. Math. Anal., 35 (5), 1213 –1240 (2004). http://dx.doi.org/10.1137/S0036141002417814 SJMAAH 0036-1410 Google Scholar

19.

D. FinchM. HaltmeierRakesh, “Inversion of spherical means and the wave equation in even dimensions,” SIAM J. Math. Anal., 68 (2), 392 –412 (2007). http://dx.doi.org/10.1137/070682137 SJMAAH 0036-1410 Google Scholar

20.

Y. XuD. Z. FengL. H. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography—I: planar geometry,” IEEE Trans. Med. Imag., 21 (7), 823 –828 (2002). http://dx.doi.org/10.1109/TMI.2002.801172 ITMID4 0278-0062 Google Scholar

21.

Y. XuM. H. XuL. H. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography—II: cylindrical geometry,” IEEE Trans. Med. Imag., 21 (7), 829 –833 (2002). http://dx.doi.org/10.1109/TMI.2002.801171 ITMID4 0278-0062 Google Scholar

22.

M. XuL. V. Wang, “Time-domain reconstruction for thermoacoustic tomography in a spherical geometry,” IEEE Trans. Med. Imag., 21 (7), 814 –822 (2002). http://dx.doi.org/10.1109/TMI.2002.801176 ITMID4 0278-0062 Google Scholar

23.

M. H. XuY. XuL. H. V. Wang, “Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries,” IEEE Trans. Biomed. Eng., 50 (9), 1086 –1099 (2003). http://dx.doi.org/10.1109/TBME.2003.816081 IEBEAX 0018-9294 Google Scholar

24.

M. XuL. V. Wang, “Pulsed-microwave-induced thermoacoustic tomography: filtered back-projection in a circular measurement configuration,” Med. Phys., 29 (8), 1661 –1669 (2002). http://dx.doi.org/10.1118/1.1493778 MPHYA6 0094-2405 Google Scholar

25.

M. XuL. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 (1), 016706 (2005). http://dx.doi.org/10.1103/PhysRevE.71.016706 PLEEE8 1063-651X Google Scholar

26.

C. ZhangY. Wang, “Deconvolution reconstruction of full-view and limited-view photoacoustic tomography: a simulation study,” J. Opt. Soc. Am. A, 25 (10), 2436 –2443 (2008). http://dx.doi.org/10.1364/JOSAA.25.002436 JOAOD6 0740-3232 Google Scholar

27.

C. ZhangC. LiL. V. Wang, “Fast and robust deconvolution-based image reconstruction for photoacoustic tomography in circular geometry experimental validation,” IEEE Photonics J., 2 (1), 57 –66 (2010). http://dx.doi.org/10.1109/JPHOT.2010.2042801 1943-0655 Google Scholar

28.

L. A. Kunyansky, “Explicit inversion formulae for the spherical mean Radon transform,” Inverse Probl., 23 (1), 373 –383 (2007). http://dx.doi.org/10.1088/0266-5611/23/1/021 INPEEY 0266-5611 Google Scholar

29.

L. A. Kunyansky, “A series solution and a fast algorithm for the inversion of the spherical mean Radon transform,” Inverse Probl., 23 (6), S11 –S20 (2007). http://dx.doi.org/10.1088/0266-5611/23/6/S02 INPEEY 0266-5611 Google Scholar

30.

J. Zhanget al., “Effects of different imaging models on least-squares image reconstruction accuracy in photoacoustic tomography,” IEEE Trans. Med. Imag., 28 (11), 1781 –1790 (2009). http://dx.doi.org/10.1109/TMI.2009.2024082 ITMID4 0278-0062 Google Scholar

31.

G. Paltaufet al., “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am., 112 (4), 1536 –1544 (2002). http://dx.doi.org/10.1121/1.1501898 JASMAN 0001-4966 Google Scholar

32.

E. CandesJ. RombergT. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, 52 (2), 489 –509 (2006). http://dx.doi.org/10.1109/TIT.2005.862083 IETTAW 0018-9448 Google Scholar

33.

J. ProvostF. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imag., 28 (4), 585 –594 (2009). http://dx.doi.org/10.1109/TMI.2008.2007825 ITMID4 0278-0062 Google Scholar

34.

Z. Guoet al., “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt., 15 (2), 021311 (2010). http://dx.doi.org/10.1117/1.3381187 JBOPFO 1083-3668 Google Scholar

35.

A. Rosenthalet al., “Efficient framework for model-based tomographic image reconstruction using wavelet packets,” IEEE Trans. Med. Imag., 31 (7), 1346 –1357 (2012). http://dx.doi.org/10.1109/TMI.2012.2187917 ITMID4 0278-0062 Google Scholar

36.

Y. ZhangY. WangC. Zhang, “Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction,” Ultrasonics, 52 (8), 1046 –1055 (2012). http://dx.doi.org/10.1016/j.ultras.2012.08.012 ULTRA3 0041-624X Google Scholar

Citation Download Citation

Yan Zhang, Yuanyuan Wang, and Chen Zhang "Efficient discrete cosine transform model–based algorithm for photoacoustic image reconstruction," Journal of Biomedical Optics 18(6), 066008 (3 June 2013). https://doi.org/10.1117/1.JBO.18.6.066008

Published: 3 June 2013

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 13 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Reconstruction algorithms

Model-based design

Algorithm development

Image restoration

Photoacoustic spectroscopy

Signal to noise ratio

Absorption

1.

Introduction

2.

Method

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

3.

Theoretical Calculation Complexity

4.

Simulation

Fig. 1

4.1.

Selection of the Threshold

Fig. 2

Table 1

Fig. 3

Fig. 4

4.2.

Few-View Sampling Condition

Fig. 5

Fig. 6

4.3.

Robustness

Fig. 7

5.

In Vitro Experiments

Fig. 8

Fig. 9

Fig. 10

6.

Conclusion

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years