2.1.

Principal Component Analysis Model

PCA is a useful technique to reduce the dimensionality of the data by transforming the data to a new data domain called the feature space. The aim of PCA is to find a new basis in a low-dimensional space to represent the major features of the high-dimensional signals by using training samples.

If the number of training images is $M$, all training images can be described as a two-dimensional (2-D) matrix $\mathbf{X}=({\mathbf{x}}_{\mathbf{1}},{\mathbf{x}}_{\mathbf{2}},\cdots ,{\mathbf{x}}_M)$, where each training image is a row vector of the matrix with size $N$ ($N$: pixel number of the image) that can be expressed as follows:

The process using PCA to recover signals through training samples can be summarized in the following five steps:

2.2.

Photoacoustic Computed Tomography Reconstruction via Principal Component Analysis

The PCA-PACT reconstruction process is illustrated by the flowchart shown in Fig. 1(a). This flowchart can be divided into three stages. (1) Sampling. In this stage, two sampling strategies of full sampling and sparse sampling are used to acquire the 3-D data of the entire imaging region. Here, full sampling means that all transducer elements of the ultrasonic array are used to collect a B-scan image; otherwise, it is called sparse sampling (a 48-element ultrasonic array is shown as an example in Fig. 1). In our work, two kinds of sampling modes were explored—one of every two and three frames along the elevational scanning direction was fully sampled, and the illustration of them is shown in Fig. 1(b). The largest number of sparse-sampling frames required to obtain satisfactory reconstructed results was $\sim 67\%$ of all frames in our imaging experiments. (2) Construct the PCA basis. To determine the basis vectors of PCA, the B-scan images reconstructed by the BP method with full sampling are used as training samples to construct the sampling matrix $\mathbf{A}$. Each row of $\mathbf{A}$ is a row-vector difference image between the sample and the average image. Then, the covariance matrix can be obtained by multiplying the transpose of the sampling matrix and itself, and the eigenvalues and eigenvectors of the covariance matrix are computed by eigenvector decomposition. According to PCA theory, the PCA basis vectors are formed by a subset of the eigenvectors, and the subset can be determined by selecting the $k$ eigenvectors corresponding to the $k$ largest eigenvalues ($k$ can be determined by tests in practical experiments). (3) PCA reconstruction. The PCA projection matrix is constructed by using the PCA basis vectors determined in the second stage, with each basis vector as a column of the matrix. Next, the reconstructed B-scan images with sparse sampling are transformed into its feature space by the projection matrix using Eq. (2), and then, a more accurate 2-D photoacoustic image can be recovered by the inverse process of PCA using Eq. (3). Finally, the 3-D photoacoustic images can be constructed by using the reconstructed 2-D B-scans.

Fig. 1

Illustration of 3-D PCA-PACT. (a) Flowchart of the reconstruction process and (b) illustration of two kinds of sampling modes.

3.1.

Imaging System

A linear-array PACT system was used in this work, and its major components included a tunable dye laser pumped by a Q-switched Nd: YLF laser, a linear ultrasonic array with 48 elements (30-MHz center frequency and 70% bandwidth), and a multicore PC equipped with an 8-channel PCI DAQ card. Specifically, 570-nm light was used to illuminate the sample surface; this wavelength corresponds to an isosbestic point where oxy- and deoxyhemoglobin molecules have the same molar optical absorption coefficient. Then, the excited pressure waves were detected by the linear ultrasonic array. One B-scan was obtained from six pulses and stored in the PC by the DAQ, and 3-D data were acquired by mechanical scanning. Finally, the B-scan images were reconstructed with the proper algorithms (such as BP), and 3-D images were generated by using the maximum amplitude projection (MAP) technique. The human hand experiments described here were conducted in compliance with Washington University-approved protocols.

3.2.

In Vivo Human Hand Experiments

Noninvasive photoacoustic imaging of a human hand was performed to validate the developed PCA-PACT method. The 3-D photoacoustic image of the human hand consisted of 166 B-scans, and each reconstructed B-scan image was represented by $128\times 128\text{pixels}$. Figure 2 shows the photoacoustic images with measurements from 48, 16, and 12 transducer elements, reconstructed by different strategies. The image acquired by the traditional BP method with data from all 48 transducer elements is assumed as the control [Fig. 2(a)]. Figures 2(b) and 2(c) show the corresponding image reconstructed by BP with 16 and 12 transducer elements, respectively, and Figs. 2(d)–2(g) are the images reconstructed with PCA using different sampling data. Note that the images in Figs. 2(a)–2(g) are shown as the MAP—the maximum photoacoustic amplitudes projected along the depth direction of the surface of the hand. To further compare the reconstructed images, representative B-scan images corresponding to the cross-sections indicated by dash lines in Figs. 2(a)–2(g) are also shown [Figs. 2(a1)–2(g2)]. Compared with the control, the images reconstructed by BP with data from only 16 or 12 transducer elements became much worse [Figs. 2(b) and 2(c)]—many vascular signals were almost buried in the noise. However, the images reconstructed by PCA were dramatically improved [Figs. 2(d)–2(g)]. When one of every two frames was sparse sampled (i.e., 83 training samples were used) with measurements from only 16 transducer elements, the reconstructed photoacoustic images with PCA [Figs. 2(d), 2(d1), and 2(d2)] had almost the same reconstruction accuracy as the control with less background noise. Furthermore, when two of every three frames were sparse sampled (i.e., 56 training samples were used) with the measurements from only 12 transducer elements, the reconstructed results [Figs. 2(g), 2(g1), and 2(g2)] had almost the same quality as the control with few sparse-sampling artifacts. Thus, PCA-PACT can recover high-quality 3-D photoacoustic images, using fewer measurements than the conventional BP method. In the above experiments, 82 and 55 eigenvectors corresponding to the positive eigenvalues were chosen to construct the PCA bases for 83 and 56 training samples being used, respectively. To further clarify the principal components of PCA, images of 20 principal eigenvectors from 83 training samples are listed in Fig. 6, Appendix.

Fig. 2

Photoacoustic images of the subcutaneous vasculature of a human hand. (a–c) MAP images reconstructed by the BP method with data from 48, 16, and 12 transducer elements; (d, e) MAP images reconstructed by PCA using 83 training samples with 16 and 12 transducer elements; (f, g) MAP images reconstructed by PCA using 56 training samples with 16 and 12 transducer elements; (a1–g2) B-scan images along the vertical dashed lines shown in (a)–(g). $x$ is the (lateral) direction of the transducer array, $y$ is the mechanical scanning direction, and $z$ is the depth direction. The number 48, 16, or 12 indicates that the reconstruction is performed with data from 48, 16, or 12 transducer elements, respectively, and 83 or 56 is the number of training samples used in the PCA; the color scale represents relative optical absorption.

3.3.

Quantitative Analysis

To quantitatively evaluate the reconstructed results, two numerical analyses were performed. First, localized comparisons were performed by plotting photoacoustic amplitudes along chosen lines in B-scan images. Figures 3(a) and 3(b) are plots from B-scan images reconstructed by different strategies with data from 16 transducer elements and were obtained from the horizontal dashed lines shown in Figs. 2(a1) and 2(a2), respectively. Figures 3(c) and 3(d) are the plots from reconstructed images with data from 12 transducer elements along the horizontal dashed lines shown in Figs. 2(a1) and 2(a2), respectively. For a better quantitative comparison, the contrast-to-noise ratios (CNRs) of selected signal peaks were also computed (Fig. 3, insets). In our computations, the range from 4.3 to 5.7 mm along the $x$-axis was used as an estimation of the background. It can be seen that the images recovered by PCA had higher CNRs than those obtained by BP. In particular, when only 12 transducer elements and 56 training samples were used, the photoacoustic images reconstructed by PCA exhibited a higher CNR than the control, which demonstrated the superiority of the PCA method in suppressing the reconstruction artifacts of PACT with fewer measurements.

Fig. 3

Plots of photoacoustic amplitudes (relative optical absorption) along the chosen dashed lines in B-scan images. (a, b) Plots of photoacoustic amplitudes from B-scans reconstructed by different strategies with data from 16 transducer elements along the horizontal dashed lines shown in Figs. 2(a1) and 2(a2), respectively; (c, d) plots of photoacoustic amplitudes from images reconstructed with data from 12 transducer elements along the dashed lines shown in Figs. 2(a1) and 2(a2), respectively. Insets: the CNRs of selected signal peaks; the number 48, 16, or 12 after each reconstruction method indicates the number of transducer elements of the data used in the reconstruction, and 56 or 83 is the number of training samples used in the PCA.

The second quantitative approach for comparison was the relative error (Rerr), which is defined as

Fig. 4

Plots of the relative errors of all 166 frames for in vivo human hand imaging. (a, b) Error plots between the control and images of BP16 and PCA16 using 56 and 83 training samples, respectively; (c, d) error plots between the control and images of BP12 and PCA12 using 56 and 83 training samples, respectively. The number 16 or 12 after each reconstruction method indicates that the reconstruction is performed with data from 16 or 12 transducer elements; 56 or 83 is the number of training samples used in the PCA.

3.4.

Principal Component Analysis Versus Compressed Sensing for Photoacoustic Computed Tomography

CS is currently used for photoacoustic imaging to improve its imaging quality when using fewer measurements. The imaging abilities of the PCA and CS methods were compared to further evaluate the performance of the proposed PCA-PACT method. Figures 5(A)–5(C) are MAP images of the vasculature of the human hand reconstructed by BP with data from 48 transducer elements, CS with data from 16 transducer elements, and PCA with data from 16 transducer elements using 56 training samples, respectively. Figures 5(a)–5(c) are B-scan images of the 40th frame reconstructed by the above three methods, respectively, and Figs. 5(d)–5(f) are the corresponding reconstructed results of the 110th frame. With the use of CS, the image quality improves significantly [Figs. 5(B), 5(b), and 5(e)], although some sparse-sampling artifacts are still notable. Surprisingly, when PCA was used, the higher quality photoacoustic images are recovered, and the sparse-sampling artifacts are effectively suppressed [Figs. 5(C), 5(c), and 5(f)]. To compare the two methods, the CNRs of the B-scan images reconstructed by PCA and CS were also computed and are listed in Table 1 [note: the positions of the peak signals and background region used in the CNR computations are indicated in Figs. 5(a) and 5(d)]. In this case, the CNRs of the images reconstructed by PCA were approximately twice as large as those for BP. Further, the difference images between the results reconstructed by CS and PCA and the control were generated and are shown in Fig. 5. Figures 5(I1) and 5(I2) are the difference images between CS16 and the control of the 40th and 110th frames, respectively, and Figs. 5(I3) and 5(I4) are those between PCA16_56 and the control of the two frames. Compared with Figs. 5(I1) and 5(I2), Figs. 5(I3) and 5(I4) are more similar to the background noise, suggesting that the images reconstructed using the PCA method possessing higher fidelity.

Fig. 5

Photoacoustic images reconstructed by CS and PCA methods. (A–C) MAP images reconstructed by BP with data from 48 transducer elements, CS with data from 16 transducer elements, and PCA with data from 16 transducer elements and 56 training samples; (a, d) the 40th and 110th B-scan images from (A); (b, e) the 40th and 110th B-scan images from (B); (c, f) the 40th and 110th B-scan images from (C); (I1, I2) difference images between BP48 and CS16 of the 40th and 110th frames; (I3, I4) difference images between BP48 and PCA16_56 of the 40th and 110th frames. The number 48 or 16 after each reconstruction method indicates that the reconstruction is performed with data from 48 or 16 transducer elements, respectively; 56 is the number of training samples used in the PCA.

Table 1

Comparisons of the performance of different reconstruction methods.

Moreover, the imaging times were computed for the BP, CS and PCA methods and are listed in Table 1 as a measure of the imaging speed (note: all the statistical computations were done on a PC with an Intel Core G640 CPU of 2.8 GHz). For different reconstruction methods with different measurements and training samples, the corresponding image reconstruction time and DAQ time for one frame are all presented in this table. The following conclusions can be drawn from this data: (1) PCA16_56 has an $\sim 50\%$ shorter DAQ time than that of the methods with data from 48 transducer elements and an $\sim 1.5$-fold longer DAQ time than that of the methods with data from 16 transducer elements; (2) PCA16_56 has a nearly eightfold faster reconstruction time than the CS method with data from 16 transducer elements and $\sim 50\%$ faster reconstruction time than BP with data from 48 transducer elements; (3) For the sum of the DAQ and image reconstruction times, the overall 3-D imaging speed of PCA16_56 is ∼sevenfold faster than the speed of the CS method with data from 16 transducer elements and $\sim 1.8$-fold faster than the speed of the BP method with data from 48 transducer elements.