Enhancing sparse-view photoacoustic tomography with combined virtually parallel projecting and spatially adaptive filtering

To fully realize the potential of photoacoustic tomography (PAT) in preclinical and clinical applications, rapid measurements and robust reconstructions are needed. Sparse-view measurements have been adopted effectively to accelerate the data acquisition. However, since the reconstruction from the sparse-view sampling data is challenging, both the effective measurement and the appropriate reconstruction should be taken into account. In this study, we present an iterative sparse-view PAT reconstruction scheme, where a concept of virtual parallel-projection matching the measurement condition is introduced to aid the “compressive sensing” in the reconstruction procedure, and meanwhile, the non-local spatially adaptive filtering exploring the a priori information of the mutual similarities in natural images is adopted to recover the unknowns in the transformed sparse domain. Consequently, the reconstructed images with the proposed sparse-view scheme can be evidently improved in comparison to those with the universal back-projection method, for the cases of same sparse views. The proposed approach has been validated by the simulations and ex vivo experiments, which exhibits desirable performances in image fidelity even from a small number of measuring positions.


Introduction
Photoacoustic tomography (PAT) is emerging as a powerful technique for providing deeptissue structural, functional and molecular information in both small animal [1][2][3] and human imaging studies [4][5][6]. In PAT, unfocused ultrasonic transducers, including the single-element or ultrasonic arrays, are usually adopted to receive the photoacoustic (PA) signals emitted from biological tissues. However, a normal ultrasonic array comprising several hundred elements would be relatively expensive, since each element requires its exclusive preamplifier and data acquisition channel. Thus, using one or several single-element transducers for circularly scanning the imaging object is the usual choice in many experimental setups in view of its comparative cost-effectiveness and high measurement flexibility [7,8]. Nevertheless, these schemes normally require hundreds of scanning steps to acquire full-view PAT projections and usually take several minutes. To achieve rapid scanning for increasing demands in practice, sparse-view PAT measurements can be utilized, where the scanning time can be greatly shorten by increasing the interval of the projection angles. Naturally, due to the sparse-view measurements, the conventional algorithms, e.g., the universal back-projection (UBP) [9] or the model-based (also matrix-based) inversion without sparse regularization [10], usually lead to blurred and distorted reconstructions.
In order to guarantee the reconstructed images with high quality in sparse-view PAT, some reconstruction algorithms driven by the concept of compressed sensing (CS) have been studied. At present, the CS-based approaches have been successfully applied to the biomedical imaging field, most notably for Magnetic Resonance Imaging (MRI) and Computed Tomography (CT) [11]. In 2009, the idea of introducing CS theory into PAT was first proposed by Provost and Lesage and was tested in phantom experiments [12]. Thereafter, many efforts about CS-based PAT have been made, and some promising results of phantom, ex vivo or in vivo experiments have also been presented. Most of them import certain a priori information, such as the sparsity of the object representation, to the model-based inverse problem as the regularization item, and then recover images by using the iteration approaches [13][14][15][16][17][18][19]. Commonly, the total variation (TV) as a valid sparsity regularization term utilizing the sparsity of the natural image gradient has been adopted in CS-based PAT reconstructions [13][14][15][16][17]. Moreover, Meng et al. have proposed a reconstruction framework of "compressed sensing with partially known support", which uses a small part of the known nonzero-signals' locations in the transformed sparse domain as a priori information [20]. She has also developed a principal component-analysis-based PAT to reduce the scale of the reconstruction [21]. In addition, Sandbichler et al. have presented a reconstruction using the sparsely transformed integrating-line-detector data followed by the UBP process [22]. Besides, some approaches using patterned excitation light or different scanning schemes have also been proposed [23,24]. The demand of the model-based PAT for the computational resources is large, and in particular, extensive memory space is needed to store the model matrix [25]. Furthermore, it is important but difficult to choose a suitable regularization parameter to balance the fidelity term and the regularization term during the reconstruction. Moreover, the TV-based regularization methods are only the measure of local changes in images and tend to produce over-smoothed image edges and texture details, since these methods favor the piece-wise constant solutions [26]. In addition, when the imaging objects are some small and narrow targets, such as the tiny blood vessels, the TV-based methods may not be effective because the image gradient is not enough sparse in these cases.
Based on the above analysis, to enhance the image quality in sparse-view PAT, both the sparse representation method and the regularization procedure should be optimized. In this study, we present an iterative sparse-view PAT reconstruction scheme that is described as follows. Under a certain measurement condition, the received signals from one measuring position can be viewed as the "virtual parallel-projection profile" of the initial pressure distribution (P 0 image) along the direction orthogonal to the measuring direction (defined as the position vector of the transducer). In this case, similar to the parallel-beam CT imaging, the received PA signals from each measuring position can be transformed to the Fourier domain, and the transformed ones can be considered as the partial Fourier spectrum of the image to be recovered according to the central-slice theorem. In this way, the "compression" process can be achieved. Although a high-fidelity reconstruction can be achieved using a part of the Fourier spectrum, in some sparse-view PAT scenarios, however, the measured spectrum is still insufficient. In principle, it is significantly meaningful to develop a method that is able to effectively exploit the unobserved partial spectrum. For this goal, we propose herein an iterative reconstruction approach. In each iteration step, a spatially adaptive filtering procedure, known as Block-matching and 3D filtering (BM3D) [27], is applied to the image domain to effectively exploit the new features and details of the object. Then the filtered image is transformed back to the Fourier domain. By keeping the part of measured Fourier spectrum unchanged, the next iteration will be executed. The BM3D filtering as a sparse representation process fully considers the a priori information of the existing mutually similar blocks in natural images, and it is a non-local filtering that utilizes more comprehensive feature information of the imaging object compared with the local-based methods. The proposed approach has been validated by the simulation and ex vivo experiments, exhibiting promising performances in imaging fidelity even from a small number of measuring positions.

Virtual parallel-projection
In CT reconstruction, one of the most fundamental concepts is the central-slice theorem. For an imaging object, ( , ) x y F ω ω is its two-dimensional Fourier-transform (2D-FT) result. This theorem states that the values of one-dimensional Fourier transform (1D-FT) of the projection ( , ) P ω α of this object and the values of the profile along a line drawn through the center of ( , ) x y F ω ω are the same. The formula of the central-slice theorem takes the following form where ω is the frequency component in the 1D Fourier domain; α is the angle of one projection; x ω and y ω represent the x-direction and y-direction frequency components in the 2D Fourier domain, respectively. The central-slice theorem provides a direct relationship between the measured signals and the imaging object in the transformed domain where only a part of coefficients are required to reconstruct the image. In 2D-PAT reconstruction, the received signals from one measuring position can be viewed as the projections along different arcs centered on the transducer, as shown in Fig.  1(a). In practice, if the rotation radius and the size of the transducer are appropriately increased, the arc projection path can become linear. This can be regarded as a "virtual parallel-projection" process, as depicted in Fig. 1(b). In this situation, similar to the CT imaging, the measured PA signals can be transformed to the Fourier domain, and the transformed signals can be considered as the partial Fourier spectrum of P 0 . Then, if the remaining part of unknown spectrum can be effectively obtained, the high-quality PAT image is able to be reconstructed.
The determination of the rotation radius is important. The too short scanning distance can't guarantee the "virtual parallel-projection" condition, while the too long scanning distance would lead to a serious decay of the PA signal, especially for the high-frequency components. The consideration for selecting the rotation radius is described as follows, and Fig. 2 is used to help illustrate this process: (1) For a transducer with a large active size (D), the received signals from one sampling position can be regarded as the superposition of the signals from a series of virtual point-transducers arranged along the surface of the transducer, as illustrated in Fig. 2(a). This superposition process leads to a change from the multiple arc back-projection paths to an approximate straight line [29][30][31][32], especially in the central area (with the same size as D) of the reconstruction grid; (2) The PA signals are back projected in the given reconstruction grid from the measuring position. For the column of the grid closest to the measuring position (Column 1 in Fig. 2(a)), if one synthetic projection path is almost completely within Column 1, the other synthetic projection paths further away from the same measuring position are more like straight lines.
As depicted in Fig. 2(b), D 1 is regarded as a virtual-point detector located at the edge of the transducer. E is the edge point of the reconstruction grid with the side length of L. Point F and point G are separately located on the two boundary lines of Column 1, and 1

EG GD
⊥ . The distance between F and G is dl, which represents the minimum resolution of the reconstruction grid. Taking (1) and (2) into account, if one synthetic projection path is almost entirely within Column 1, and it can pass through the upper and lower boundaries of Column 1, at least through the edge points of the boundaries, then the "virtual parallel-projection" condition can be approximately satisfied. It means that if one back-projection path from D 1 can pass through E and F, then the rotation radius (R) calculated in this case is considered to be the minimum rotation radius. R is formulated as follows: where R1 is the distance between D 1 and E, and The derivation of R1 is detailed in the Appendix A.
In our experiments, L is set to 22 mm for the small-animal model, and dl is set slightly less than the axial resolution of the transducer [33]. The central frequency of the selected unfocused transducer is 2.25 MHz (V306, Olympus NDT), thus dl is set to 100 μm accordingly, and D of this transducer is about 12.5 mm. In this case, the scanning radius R is finally selected as 155 mm. Picking three points in one PA signal sequence received from one measuring position and using the above settings, we carry out the simulation to exhibit the projection paths, as shown in Fig. 3(b). For comparison, the projection paths back projected from the point transducer have also been depicted in Fig. 3(a). It can be observed from the results that, using the transducer with a larger active size and the suitable rotation radius can make the arcprojection path become the line-like one. Furthermore, the "virtual parallel-projection" condition is basically satisfied.

Iterative sparse-view PAT reconstruction cooperated with BM3D filtering
In order to exploit the unobserved partial Fourier spectrum of the P 0 image, an iterative sparse-view PAT reconstruction scheme cooperating with BM3D filtering is proposed.

BM3D filtering
The BM3D filtering [27, 34] consists of two main steps: Step 1: Taking the original image ( 0 u ) as the input, an intermediate (i.e., basic estimate) image ( basic u ) is estimated using hard thresholding during the collaborative filtering process: (a) Grouping and collaborative filtering: The input image 0 u is processed block by block. These blocks are called reference blocks. For each reference block, similar blocks to the currently processed one are found using a similarity measure (block matching). A threedimensional (3D) group (array) is built by stacking the matched blocks, and then a collaborative filtering is applied to the grouped blocks. In this step, hard thresholding is applied in shrinking the coefficients in the transform domain. . For each reference block, record the locations of all matched blocks. Use these locations to form two 3D groups (arrays) of image blocks, one from 0 u and the other from basic u . Then, apply a 3D transform on both groups. The collaborative Wiener filtering is performed on the first group ( 0 u -group). Here, the energy spectrum of basic u is used as the true (pilot) energy spectrum in calculating the empirical Wiener coefficients. The estimates of all grouped blocks are obtained by applying the inverse 3D transform on the filtered coefficients and returning to their original locations.
(b) Aggregation: final u is obtained by aggregating all the block-wise estimates using a weighted average. In BM3D filtering process, as the grouped image blocks are similar, the transformation can achieve a very high sparse representation of the original image, and the collaborative filtering is able to reveal the finest details shared by the grouped blocks. Thus, during the PAT reconstruction, the BM3D filtering can help to enhance the useful features of the P 0 image and promote the reconstruction fidelity in the image domain.

Iterative reconstruction process
The unobserved spectrum in Fourier domain can be gradually estimated by an iterative process that is further detailed in the following pseudo codes. At the beginning of the iteration process, we first reconstruct the initial P 0 image from the projection data by using the inverse radon transform (IRT) that has the following form [28]: where ( , ) f x y is the P 0 image to be reconstructed in our work, and ( , ) p s α denotes the projection at the sampling angle α , i.e., the PA signals received from the sampling angle α .
Then the BM3D filtering process is executed (Eq. (4)). After that, the updated P 0 image is transformed to the Fourier domain (Eq. (5)). Keeping the measured Fourier spectrum unchanged (Eq. (6)), the updated 2D Fourier spectrum is transformed back to the image domain (Eq. (7)). Then the next filtering process will be carried out. The proposed method is abbreviated as IRT-BM3D in this work. To describe the implementation process of the proposed IRT-BM3D approach more vividly, the flowchart of the algorithm is also illustrated in Fig. 4.   At each iteration, some pseudo-random noise Ω is only injected in the unobserved spectrum part (Eq. (6)), which helps the BM3D filtering to attenuate the noise robustly and to reveal new features effectively and further to accelerate the convergence. This is similar to the "random search" process. To avoid the injected noise in Fourier domain affecting the global reconstructed image space, the amplitude of Ω is set to very small and exponentially decreasing with the increase of iterations. In the final iteration, the amplitude is set to 0.
According to the central-slice theorem, for the measurements from the sampling angle α   To assess the experimental performance of the proposed imaging system in terms of the spatial resolution, a phantom is built that consists of some absorbing polyethylene microspheres with a diameter of 200 μm. These microspheres are placed approximately in the same plane of the phantom. The full angle measurement with the sampling interval of o 1 is carried out, and then the image is reconstructed by the proposed IRT-BM3D method. Some of the reconstructed microspheres are shown in Fig. 6(a). The profile (along line A-A') of one recovered microsphere is depicted in Fig. 6(b). The reconstructed diameter of the microsphere is estimated as 310 μm, i.e., the average value of the full width at half maximum of three recovered microspheres (I, II, and III).

Simulation experiment
To demonstrate the efficiency and superiority of the IRT-BM3D scheme in sparse-view PAT reconstruction, the simulations are performed in 2D situation. The exact P 0 image is set as the vessel phantom, and two "tumor" targets are added to the phantom as illustrated in Fig. 7(k). , respectively. It can be observed from the images that the IRT-BM3D results are overall superior to the UBP results, especially in #60-and #30view cases. For example, there are lots of artifacts in #60-view image which is reconstructed by the UBP method ( Fig. 7(d)), and some fine structures, e.g. the branched thin vessels, become difficult to be recognized in the low-position tumor target. However, they can be clearly observed in the IRT-BM3D result (Fig. 7(i)). As for the #30-view UBP result shown in Fig. 7(e), even the large structures are blurred. For instance, the whole boundaries of the two tumor targets cannot be extracted clearly, while these structures are reconstructed with high fidelity in the IRT-BM3D image (Fig. 7(j)). Moreover, both of the UBP and the IRT-BM3D schemes obtain the high-quality results in #360-view case, which also demonstrate the accuracy of the proposed IRT-BM3D method. In order to qualitatively analyze the texture details reconstructed by different methods in each sparse-view case, the profiles of all recovered images along the yellow dashed line (A-A') drawn in Fig. 7(k) are showcased in Figs. 7(l)-7(p). These results demonstrate that there exist large fluctuations in the profiles of the UBP results because of the severe artifacts existing in the reconstructed images. In contrast, the profiles of the IRT-BM3D results are much closer to the exact ones, especially in #60-and #30-view cases. For example, the small vessel on the right side of the upper tumor target is indistinguishable in #30-view UBP result, while it is conspicuous in the IRT-BM3D result ( Fig. 7(p)). From the above analysis, the IRT-BM3D scheme is proved to be able to obtain the satisfactory imaging fidelity and structural precision even under the sparser sampling conditions. To quantitatively assess the reconstructed results, the parameter "distance (d)" is calculated to show the difference between the reconstructed image and the exact one. The parameter d is defined as ( ) where u and v are the reconstructed image and the exact phantom image, respectively. The size of the image is M N × . The smaller value of d denotes that the reconstructed image is closer to the exact one. Besides, the peak signal-to-noise ratio (PSNR) of the reconstructed result is also calculated, which is defined as ( ) where max v is the maximum gray value of phantom image, and the normalized max v is set to 1 in this work.  showcase the values of d and the PSNR for each iteration step during the IRT-BM3D reconstruction process in each sparse-view case. To better reconstruct some fine structures and to evaluate the convergence performance of the algorithm, the number of the iterations is set to 300. It can be analyzed from the results that, the convergence can be reached in each case. The smaller d and the larger PSNR in the final results compared with the initial ones demonstrate the proposed method is able to effectively reconstruct the images with high quality. Furthermore, compared with #30-view case, the final values of d are relatively close in #120-, #90-and #60-view cases, so are the PSNR values. These indicate that, to some extent, different number of sampling positions does not entirely affect the quality of the IRT-BM3D results. While with the sampling positions becoming particularly sparse, such as in #30-view case, the quality of the reconstruction is not enough good as the results in other sparse-sampling cases analyzed in our work. Nevertheless, the value of d has a rapid and large decline process in #30-view case.   In practical experiments, the measured PA signals are usually affected by the system noises. To verify the performance of the algorithm under extremely harsh noise conditions, the white Gaussian noises are added to #60-view simulated PA signals to obtain the signal-tonoise ratio (SNR) values of 20, 15, 10 and 5 dB, respectively. The results reconstructed by the UBP and IRT-BM3D methods are illustrated in Fig. 9. The values of d and PSNR of each iteration step during the IRT-BM3D reconstruction in each SNR case have also been displayed in Figs. 10(a) and 10(b). As observed from the results, with the increase of the noise level, the quality of the UBP images deteriorates seriously. For instance, some fine structures, e.g. the branched thin vessels, cannot be recognized clearly in the case of SNR = 5 dB. On the contrary, the IRT-BM3D results show better noise robustness, and all the texture details can be extracted accurately in all the noise cases. In the cases of SNR = 15 dB and SNR = 20 dB, there are no significant reduction in visual quality compared with the non-noise case. Although in the case of SNR = 5 dB, some background noises can be observed in IRT-BM3D image, the fine structures can still be distinguished obviously. The results in Fig. 10 demonstrate that the convergence can be reached in each case, and after the iterations, the PSNR and fidelity of the reconstructed images are dramatically improved. The increases of the PSNR values are 6.8604, 5.2097, 5.7316 and 5.2695 dB, which correspond to the cases of SNR = 5, 10, 15 and 20 dB, respectively. Meanwhile, in all cases, the d values are reduced by at least 0.15.

Ex vivo experiments with biological tissue
To further validate the feasibility of the proposed IRT-BM3D approach for biomedical applications, the ex vivo experiments are conducted on porcine tissues and mouse intestinal tissue, respectively.

Tumor-mimicking tissue imaging
To set up the imaging experiment for the mimic subcutaneous tumor tissue, the biological sample is designed as a sandwich shape, as shown in Figs. 11(k) and 11(l). The top and bottom layers of the sample are composed of porcine tenderloin tissue. Two small pieces of porcine liver tissue with different shapes mimicking the "tumor" targets are sandwiched between two layers. The whole sample with dimensions (length width height) × × of 20 mm 20 mm 12 mm × × is placed into the imaging chamber. Then the remaining space between the sample and chamber is filled with the agar mixture (36°C, A9414, SIGMA). The #360-, #120-, #90-, #60-and #30-view PAT measurements are carried out in the imaging plane close to the interlayer of the ex vivo sample. The measurement settings are the same as those in the simulation experiment.
The images reconstructed by the UBP and the IRT-BM3D algorithms are shown in Figs. 11(a)-11(j), respectively. It can be observed from the images that the UBP and the IRT-BM3D results have the comparable quality in #360-and #120-view cases. Nevertheless, there is a certain degree of blurs in #120-view UBP result. The fidelity of UBP images sharply declines as the sampling positions decrease. There are lots of artifacts in #90-, #60-and #30view UBP results, where the boundaries of "tumor" targets are severely blurred, as indicated by the arrows in Fig. 11(d). In contrast, the IRT-BM3D image has the sharper image edges and clearer target boundaries (Fig. 11(i)). In #30-view case, the two targets can hardly be visually distinguished in the UBP image, while they are able to be extracted in the IRT-BM3D image (Fig. 11(m)). In #360-view case, some texture details can be recovered clearly. For example, some muscle fibers of the porcine tenderloin tissue can be observed. Compared to #360-view result, although these particularly fine structures might not be preserved in #60and #30-view IRT-BM3D images, the reconstructed results are still satisfactory in image fidelity and accuracy of the target structure. Assuming #360-view IRT-BM3D result is the standard image, the PSNR and d values of the final reconstructed UBP and IRT-BM3D results in #120-, #90-, #60-and #30-view cases are calculated and listed in Table 2. As can be seen from the table that compared with the UBP results, the IRT-BM3D images have a significant increase in PSNR values and have smaller d values. The maximum increase of the PSNR value is 11.4036 dB in #120-view case, and the maximum drop of d value is 0.1646 in #30-view case.

Vascular imaging for mouse intestinal tissue
To further demonstrate the applicability of the proposed IRT-BM3D method for imaging more complex biological structures, the mouse intestinal tissue is selected because it has abundant blood vessels with various scales. A 6-week-old male healthy KM mouse is euthanized, and a portion of intestinal tissue is excised for imaging. The tissue is put in the imaging chamber and embedded in turbid agar gel, as shown in Fig. 12(a). The sample is placed in the imaging system < 30 minutes after the animal's death. The experimental procedures in this study are reviewed and approved by the subcommittee on research animal care at Tianjin Medical University Cancer Institute & Hospital. For fine-structure tissue imaging, to balance the measurement cost and the fidelity of the reconstruction, we evaluate the performance of the IRT-BM3D method in #90-view measurement case. The #90-view IRT-BM3D result is illustrated in Fig. 12(b), and the #90-view UBP image is also shown in Fig. 12(c) as comparison.
As seen from the anatomical structure of the tissue in Fig. 12(a), blood vessels of different sizes (about 150 ~500 μm) and shapes are distributed in the sample. The image quality of the #90-view UBP image is not satisfactory and only part of the fine structures can be captured. That's because with the appearance of multiple artifacts in the background, the image contrast is impaired and the visualization of fine details is blocked. This indicates that the UBP algorithm is not suitable for recovering texture details in sparse-view case of sampling interval of o 4 . In contrast, the #90-view IRT-BM3D result presents the better performance in image fidelity and structural precision. The IRT-BM3D result corresponds well to the vascular anatomy, as indicated by the arrows. The vessels with different sizes (arrows 5, 11, 12 and 8) and different shapes (the "Y"-shaped vessels indicated by arrows 2 and 9) can be recovered more clearly in the IRT-BM3D result compared with the UBP image, e.g. the regions in white dotted boxes in Fig. 12(b), and 12(c). The branch and the orientation of the vessels indicated by arrows 1, 2, 3, 8 and 12 in Fig. 12(b) are explicit, while these features cannot be accurately resolved in the #90-view UBP image.

Discussions and conclusions
In order to guarantee the reasonable instrument cost and meet the need of rapid scanning in some cases, such as in dynamic imaging, the sparse-view PAT measurement is adopted. To recover the P 0 image with high quality, the sparse-view PAT reconstruction approaches are widely studied. In this work, a novel iterative IRT-BM3D reconstruction has been proposed. The simulations have been carried out, and the experimental validations have been performed on a self-built PAT imaging system.
For the IRT-BM3D approach, first, we have proposed the "virtual parallel-projection" measurement condition. Under this condition, the arc back-projection paths tend to the approximately parallel lines, as verified and illustrated in Fig. 3. Thus, similar to the parallelbeam CT reconstruction, following the central-slice theorem, the measured signals can be transformed to the 2D Fourier domain where the P 0 image can be sparsely represented, and the transformed signals can be directly acted as the partial Fourier spectrum of the P 0 image. Furthermore, the unobserved partial spectrum can be gradually exploited by iteratively executing the BM3D filtering on the image which is obtained by transforming the original "measured" and the "newly explored" Fourier spectrum to the image domain. This process utilizes another sparse representation assumption of the P 0 image: The natural images consist of a few templates, and these templates can be vessels and any other tissues that have the similar appearance in image. The BM3D filtering process is able to enhance the similarity between blocks that can be described by one template, which leads to a sparse representation with a few block templates for exploiting new features and improving the imaging quality effectively.
As the commonly used sparse-view PAT reconstruction scheme, the model-based inversions with the TV-based regularization item (MB-TV) are widely investigated. To compare the performances of the proposed IRT-BM3D method and the common MB-TV scheme, the simulation experiments of the vessel phantom in #60-and #30-view sampling cases are performed. The images reconstructed by the UBP, TV-based gradient descent (TV-GD) [13] and IRT-BM3D algorithms are shown in Fig. 13, respectively. The same reconstructions have also been carried out on the tumor-mimicking tissue sample, as exhibited in Fig. 14. For the TV-GD reconstruction, the regularization parameter is set to 0.2 [13], and the criterion for the iteration termination is that the parameter d stops decreasing for two successive updating stages. The parameters for the IRT-BM3D reconstruction are the same as the previous part. Table 3 and Table 4 list the PSNR and d values of the final reconstructed TV-GD and IRT-BM3D results of the simulation and ex vivo experiment, respectively.   It can be seen from the results that both the TV-GD and the IRT-BM3D images are conspicuously superior to the UBP results, especially in #30-view case. While the IRT-BM3D images exhibit better performances than TV-GD results in visual observation and quantitative assessment. For example, compared with the #60-view IRT-BM3D result (Fig. 13(c)), some details, e.g. the branched thin vessels in low-position tumor target, are slightly blurred in the TV-GD image ( Fig. 13(b)). Moreover, there are more artifacts in the #30-view TV-GD result (Fig. 13(e)). In general, according to the imaging results, the TV-GD method tends to produce over-smoothed edges and texture details, especially illustrated in Fig. 14(e). This observation can be explained as the fact that the classical TV-based algorithms are prone to perceive images to be local piecewise-constant, while this violates most of the natural images. In contrast, the BM3D method, as a non-local filtering, is more appropriate for reconstruction because the selection of image blocks is not limited to the neighborhoods. Furthermore, the selection of the regularization parameter for balancing two parts of the objective function in the MB-based approach has a significant effect on imaging results, which should be deeply studied. Fortunately, this issue can be avoided in the IRT-BM3D reconstruction.
In this paper, the numerical simulations and ex vivo experiments have verified the effectiveness of the proposed IRT-BM3D algorithm. The reconstructed results having been discussed qualitatively and quantitatively show that the IRT-BM3D approach has the better performance than UBP, in terms of PSNR and imaging fidelity, especially in #30-view case which is a relatively sparse sampling case. Simulation results have also exhibited the high robustness and good convergence performance of the proposed method. The ex vivo experiments have further demonstrated that the IRT-BM3D method is competent for practical sparse-view PAT applications. The IRT-BM3D reconstructions exhibit better performances in effectively reducing the under-sampling artifacts and clearly revealing target boundaries and some texture details. For instance, in #90-view vascular imaging experiment, even the branched thin vessels can be reconstructed accurately, which is crucial for biomedical imaging diagnosis, while this fine structure is indistinguishable in the UBP result. Moreover, it should be noted that, if the imaging objects contain some fine structures, such as in vessel phantom, the quality of the #30-view IRT-BM3D reconstruction is worse than that in #120-, #90-and #60-view cases. However, in tumor-mimicking tissue imaging experiment, where the imaging target has a larger size, the reconstructed results have the similar quality in #90-, #60-and #30-view cases in terms of PSNR, d and visual observation. This indicates that selecting the suitable sparse measuring scheme according to the imaging demand is also important for obtaining good imaging results in sparse-view PAT.
In the setup with only a single transducer, the overall time of data acquisition is generally reduced with increase in the sampling interval (equivalently, the decrease in the sampling positions). Nevertheless, this reduction is not significant in the case of the sparse-view scanning, since the total scanning path of the transducer is always a full circle around the target, which needs much longer time to go through than the time for the sparse sampling. For our imaging system, the total measurement durations in #360-, #120-, #90-, #60-and #30view cases are 56.88s, 35.95s, 34.66s, 32.15s and 25.76s, respectively. In realistic scenarios, if a few more transducers are used, such as only six transducers, and these transducers are fixed on the rotation stage every o 60 , then to achieve #60-view measurement, the rotation stage only needs to rotate nine times at an interval of o 6 . In this case, the sampling time can be reduced to 4.8225s.
In the current work, the imaging system and the scanning geometry are designed for in vivo small animal imaging. Nevertheless, if the imaging target has a larger size, such as the breast, this scanning geometry impractically requires a larger rotation radius to satisfy the "virtual parallel-projection" condition. In this case, instead of increasing the rotation radius, increasing the active size of the transducer, e.g., using an integrating line detector [30], can be the better alternative, as aforementioned.