Sparse photoacoustic microscopy based on low-rank matrix approximation

Ting Liu; Mingjian Sun; Naizhang Feng; Minghua Wang; Deying Chen; Yi Shen

doi:10.3788/COL201614.091701

Abstract

As a high-resulotion biological imaging technology, photoacoustic microscopy (PAM) is difficult to use in real-time imaging due to the long data acquisition time. Herein, a fast data acquisition and image recovery method named sparse PAM based on a low-rank matrix approximation is proposed. Specifically, the process to recover the final image from incomplete data is formulated into a low-rank matrix completion framework, and the “Go Decomposition” algorithm is utilized to solve the problem. Finally, both simulated and real PAM experiments are conducted to verify the performance of the proposed method and demonstrate clinical potential for many biological diseases.

Photoacoustic microscopy (PAM) is probably the fastest-growing area of biomedical imaging technology due to its tremendous potential for use in various biomedical studies and research applications^[1]. It is a non-ionizing and non-invasive hybrid imaging technique that combines strong optical absorption contrast with high ultrasonic spatial resolution^[2,3]. Using intrinsic contrasts such as hemoglobin or melanin, PAM can provide both morphological and functional information^[4]. Moreover, by using molecularly targeted exogenous contrast agents, it can also provide molecular information. It has been proven that the spatial resolution of PAM images is adequate for many biological applications^[5,6]. However, to resolve smaller structures, such as capillaries in a microcirculation system, a much higher spatial resolution is required^[7,8].

Among the existing PAM technologies, the spatial resolution almost solely depends on the ultrasonic parameters, including the center frequency, the bandwidth of the transducer, and the numerical aperture of the objective. Thus, in order to improve the spatial resolution, the hardware requirements of the system will always be increased correspondingly^[9,10]. Many efforts have been made to improve the spatial resolution without increasing the cost of the system, resulting in some advanced data acquisition schemes and novel signal processing methods^[11,12].

Generally, a point-by-point scanning scheme is implemented in the conventional PAM system. More sampling points are employed to acquire a higher-resolution image, which is confirmed to be more time consuming. Recently, many investigations have proved that the photoacoustic image has a spatial sparsity property in a certain domain by an appropriate sparse transform^[13–15]. Additional, the property has been thoughtfully explored, and the final image is reconstructed from far fewer measurements^[12]. In this Letter, we make use of a simple uniform random sampling scheme to accelerate the PAM data acquisition. The final PAM image is recovered directly from the incomplete data. Specifically, what we keep focus on are the relationships among the sampling density and the rank of the image with the corresponding image recovery method.

Sign up for Chinese Optics Letters TOC. Get the latest issue of Chinese Optics Letters delivered right to you！Sign up now

Based on the random sampling scheme, many missing values will occur in the sampling data. Fortunately, matrix completion has been shown recently to be effective in estimating the missing values in a matrix from a small fraction of known entries^[16,17]. At present, it has been utilized for many applications, including image processing, noise reduction, information retrieval, computer vision, machine learning, principal component analysis, and so on. To address this ill-posed problem, optimization methods for matrix completion often assume that the recoved matrix has a low-rank structure; then, they use this as a constraint to resolve the minimization problem.

The low-rank matrix approximation method is a fast way to solve the low-rank matrix completion. Given a matrix $A \in R^{m \times n}$ , finding its low-rank approximate matrix $B$ is a critical task in many technical fields. It minimizes ${‖ A - B ‖}_{F}^{2}$ with rank matrix $B$ limited to an integer $k$ . The problem of low-rank matrix approximation can be easily solved with the truncated singular-Value decomposition (SVD)^[18,19]. Herein, the “Go Decomposition” GoDec algorithm, as a simple low-rank matrix approximation method, is utilized to obtain the final image without loss of resolution, which can simultaneously reduce the measurement noise. In this method, the bilateral random projections (BRP)^[20] of matrix $X \in R^{m \times n}$ are used to replace the SVD to significantly reduce the time cost^[21]. The GoDec algorithm has been proven to be robust to the measurement noise. Moreover, the algorithm achieves high precision but requires less storage and memory. It also can be available for low-rank and sparse decomposition and computer vision applications.

There are two important parameters in GoDec: these are the rank of the image and the total sampling number. In much of the literature, both of the parameters are defined by users’ experience^[22,23]. Herein, two experience formulations have been proposed by many of experiments. Recently, a sparse PAM based on a low-rank matrix approximation (SLRA-PAM) system was investigated for biomedical applications. Fast acquisition has been achieved, and the PAM image was recovered by low-rank matrix approximation method.

As we have mentioned above, the proposed system is no different from the conventional PAM system for the hardware requirements. However, the data acquisition process is a little different. The schematic of the data acquisition process of the proposed system is demonstrated in Fig. 1. For example, the field of view is $m \times n$ pixels for a conventional point-by-point scanning scheme, and $X$ is the measurement data of the PAM image in this case. For the sparse sampling mode, a random sampling mask $M$ , which is a 0,1 matrix, is generated by the computer with a user-defined sampling density. Based on the sampling mask $M$ , the X, Y linear stage is controlled by a controller to achieve sparse optical scanning. Thus, a compressive measurement PAM data $Y_{M}$ , could be acquired under this scheme, which is defined as $Y_{M} = {\begin{matrix} X_{i, j} & if M_{i, j} = 1 \\ 0 & otherwise \end{matrix} .$ (1)

Figure 1.Schematic of the data acquisition process of the proposed system.

Obviously, many values are missing in $Y_{M}$ . Then, what we need to do is recover the final PAM image $X \in R^{m \times n}$ from $Y_{M}$ . It can be regarded as a matrix completion problem.

The sparsity of the photoacoustic image has been verified in many experiments^{[14,15,24–27]}. More recently, the low-rank property has been incorporated into the field of photoacoustic imaging^[12]. To our knowledge, photoacoustic imaging is typically used in tumor and vasculature imaging. Thus, four typical PAM images are demonstrated in Figs. 2(a)–2(d). The sizes of these images are $502 \times 502$ , $94 \times 94$ , $93 \times 93$ , and $176 \times 176$ pixels, respectively. Here, SVD is performed to verify the low-rank characteristcs. For a PAM image $X$ , ${USV}^{T}$ is the SVD of $X$ , where $S$ is a dialog matrix with all the singular values $σ_{1} \geq σ_{2} \geq \dots$ . Here, we define ${Sum}^{i}$ as the sum of the first $i$ -th singular values, which is ${Sum}^{i} = \sum_{j = 1}^{i} σ_{j} .$ (2)

Figure 2.Low-rank property analysis. (a)–(d) PAM images with tumor and blood vessels. (e)–(h) Singular values of PAM image (a)–(d).

The ratio of ${Sum}^{i}$ with the sum of all the singular values are plotted in Figs. 2(e)–2(h). From Figs. 2(e)–2(h), the ratio increases dramatically with the increasing number of the singular values. More specifically, the ratio of the first 125, 35, 31, and 42 will reach 0.9. Based on the above analysis, the low-rank property of the photoacoustic image could be proven by statistics. Therefore, the low-rank property can be utilized as a constraint to the matrix completion problem, which formed a low-rank matrix completion problem. Thus, the formulation of the SLRA-PAM system is defined as $\min rank (X) s.t. X_{i, j} = Y_{M_{i, j}}, if M_{i, j} = 1 .$ (3)

The matrix completion problem can be viewed as a matrix recovery problem, where one has to recover the missing entries of a matrix with limited number of known entries. In order to solve the low-rank matrix completion problem, Zhou and Tao have proposed an equivalent formulation of Eq. (3) in Ref. [21], which is $\min {‖ Y_{M} - X - Z ‖}_{F}^{2} s.t. rank (X) \leq r supp (Z) = \bar{Ω},$ (4)where $Z$ is an estimate of $- P_{\bar{Ω}} (X)$ , and $Ω$ is the sampling index set of the matrix $X$ . $\bar{Ω}$ is the non-sampling index set of the matrix $X$ . $P_{\bar{Ω}} (X)$ is a linear operator that keeps the entries in $\bar{Ω}$ unchanged. Therefore, the GoDec algorithm is used to solve the low-rank matrix completion problem.

For GoDec algorithm, the optimization problem [Eq. (4)] can be solved instead of solving the following two subproblems alternatively: $X_{t} = \underset{rank (L) \leq r}{\arg \min} {‖ Y - X - Z_{t - 1} ‖}_{F}^{2};$ (5) $Z_{t} = \underset{card (S) \leq k}{\arg \min} {‖ Y - X_{t} - Z ‖}_{F}^{2},$ (6)where $card (S)$ is the cardinality of a set, which is a measure of the “number of elements of the set.” A BRP based low-rank approximation is used to solve the subproblem [Eq. (5)], which can efficiently accelerate the calculation speed. Assume that $\tilde{X}$ is a fast rank- $r$ BRP approximation of $X$ , and its dimension is $m \times n$ . We have $Y_{1} = \tilde{X} A_{1}$ and $Y_{2} = {\tilde{X}}^{T} A_{2}$ , where both $A_{1} \in R^{n \times r}$ and $A_{2} \in R^{m \times r}$ are random matrixes. Therefore, $\tilde{X}$ could be computed as $\tilde{X} = Y_{1} {(A_{2}^{T} Y_{1})}^{- 1} Y_{2}^{T} .$ (7)

In order to improve the approximate precision of $\tilde{X}$ , the obtained right random projection $Y_{1}$ is used to build the left projection matrix $A_{2}$ .

Table Infomation Is Not Enable

When solving the optimization problem [Eq. (3)] by GoDec, $r$ and $k$ should be determined first. Here, $r$ is the rank of recovered image and $k$ is the total sampling number of random sampling. Via extensive numerical tests, different pairs of $r$ and $k$ are utilized. In addition, structural similarity (SSIM) and peak signal-to-noise ratio (PSNR) are used to evaluate the performance of the image recovery method. These are calculated as follows: $SSIM (X, Y) = \frac{(2 μ_{X} μ_{Y} + C_{1}) (2 σ_{X Y} + C_{2})}{(μ_{X}^{2} + μ_{Y}^{2} + C_{1}) (σ_{X}^{2} + σ_{Y}^{2} + C_{2})},$ (9) $PSNR = 10 * \log_{10} \frac{M N}{\sum_{x = 1}^{M} \sum_{y = 1}^{N} {[X - Y]}^{2}},$ (10)where $X$ and $Y$ represent the reference image and the recovered image, respectively. $μ_{X}$ and $μ_{Y}$ stand for the average values of images $X$ and $Y$ , while $σ_{X}$ and $σ_{Y}$ denote the variances. $σ_{X Y}$ is the covariance between $X$ and $Y$ . All the results are summarized in Fig. 3.

Figure 3.Parameter selection. (a) The relationship with $r$ and $n$ . (b) The relationship with $k$ , $r$ , and $n$ .

Finally, the empirical formulas we acquired are shown as follows: $r = c_{1} {(n^{2})}^{c_{2}},$ (11) $k = c_{3} {(r \log (n))}^{c_{4}},$ (12)where $n$ is the dimension of matrix $X$ , and $c_{1} = 0.2276$ , $c_{2} = 0.5129$ , $c_{3} = 6.5686$ , and $c_{4} = 1.6761$ are four constants.

To evaluate the efficiency of the proposed SLRA-PAM system, both simulated and real PAM experiments are presented. In detail, the simulated experiment is performed to give a quantitative analysis, while the real PAM experiment is utilized to show the effectiveness of our method in the real case.

A simulated image of size $138 \times 138$ pixels is selected and presented in Fig. 4(a), which could be utilized as a reference image to perform the quantitative analysis. The test image has similar properties to the real PAM mage of the tumor. With further analysis, the image has a low-rank property, as shown in Fig. 4(b). The low-rank property has been verified through SVD analysis, and it is clearly seen that the first 20 singular values could nearly represent the whole image.

Figure 4.Recovery results in the simulated experiment. (a) Full sampling image. (b) Singular values of PAM image matrix. (c) Random sampling image. (d) Recovered image by optspace algorithm. (e) Recovered image by GoDec algorithm.

Based on the low-rank analysis above, it is again proven that the PAM image could be recovered by the low-rank matrix completion method directly used on a random sampling image. As analyzed above, the two parameters used in the GoDec algorithm, the rank $r$ and total sampling points $k$ , are calculated by Eqs. (11) and (12), $r = c_{1} {(n^{2})}^{c_{2}} \approx 35 k = c_{3} {(r \log (n))}^{c_{4}} = 9396 .$ Therefore, a random sampling image with a sampling density 0.4934 ( $k = 9396$ ) is demonstrated in Fig. 4(c). The corresponding recovery image by optspace and GoDec are displayed in Figs. 4(d) and 4(e), respectively. It can be clearly observed that the image could be perfectly recovered. The result of our method is better than optspace.

For further quantitative analysis, different values for $r$ and $k$ are utilized to conduct more experiments to verify our proposed method. Firstly, the total sampling number $k$ is fixed as $0.5 n^{2}$ , while different values for $r$ are employed. Additional, the SSIM is still used to show the performance. Secondly, $r$ is fixed at 26, different values for $k$ are utilized. All the results are summarized in Table 1.

k=0.5*n2=9522	r	10	15	20	25	30	35	40	50	60	70
k=0.5*n2=9522	SSIM	0.8674	0.9274	0.9386	0.9492	0.9471	0.9179	0.8415	0.7051	0.5792	0.5472
r=26	k	0.1*n2	0.2*n2	0.3*n2	0.4*n2	0.5*n2	0.6*n2	0.7*n2	0.8*n2	0.9*n2	1.0*n2
r=26	SSIM	0.1000	0.4713	0.7152	0.8693	0.9399	0.9638	0.9718	0.9741	0.9785	0.9849

Table 1. Numerical Results with Different r and k

View all Tables

As shown in Table 1, when the total sampling number is fixed at 9522, the best result is acquired with $r = 25$ . It is not far from the empirical results we acquired. Further, when $r$ is fixed at 26, better results will be acquired with larger sampling densities. Moreover, from the SSIM values, we think that when the sampling density is 0.5, the recovered image has relatively good resolution. The PSNR between the recovered image and the reference image is 42.5795 dB.

As we have described in detail, the proposed system could achieve fast data acquisition without dramatically decreasing the resolution. The system we used here is a traditional PAM system. The only difference here is the sampling scheme. We applied the SLRA-PAM system to image the micro-vasculature in the ear of a mouse in vivo. Mouse ears were chosen as the initial organ to test our system because the mouse ears have a well-developed vasculature and have been widely used to study tumors and other micro-vascular diseases.

With the conventional system, in order to give a good description of the field of view, $2000 \times 2000$ sampling points are needed. The full sampling image is shown in Fig. 5(a). The capillaries can be clearly seen in the full sampling image. For our proposed SLRA-PAM system, however, a random sampling scheme is used instead. During the experiment, the total sampling points used are 1927266, which means the sampling density is 0.48 and the data acquisition time is dramatically reduced. The incomplete PAM data is shown in Fig. 5(b). Further, the final PAM is recovered by the optspace algorithm and the GoDec algorithm, as demonstrated in Figs. 5(c) and 5(d). Here, the value of $r$ is 554. The SLRA-PAM images of the micro-vasculature agree with the conventional PAM image, and the small vessels could still be observed. It can be seen that the proposed method could achieve exciting visual results.

Figure 5.Results of real PAM system. (a) Full sampling image. (b) Random sampling image. (c) Recovered image by optspace algorithm. (d) Recovered image by GoDec algorithm.

For further analysis, the full sampling image is utilized as a reference image. The SSIM and PSNR values between the full sampling image and the recovered image by the GoDec algorithm are 0.9172 and 39.4108 dB, respectively. For comparision, the values for the optspace algorithm are 0.9085 and 38.5907 dB. Therefore, the proposed method could achieve fast data acquisition without a dramatic loss of resolution.

In this Letter, an SLRA-PAM system is proposed to accelerate the PAM data acquisition process and obtain a relatively high-resolution PAM image. Firstly, based on the low-rank analysis of typical PAM images, the SLRA-PAM is formulated into a low-rank matrix completion framework. Additionally, the GoDec algorithm, a kind of low-rank matrix approximate method, is employed to solve the optimization method effectively and exactly. In particular, two important parameters are discussed in detail to give a guide to using the image recovery algorithm. Finally, both the numerical and PAM experiments have validated that the proposed fast PAM system could acquire a relatively high-resolution image.

References

[1] L. V. Wang. Nat. Photon., 3, 503(2009).

[2] C. Kim, C. Favazza, L. V. Wang. Chem. Rev., 110, 2756(2010).

[3] C. Kim, T. N. Erpelding, K. Maslov, L. Jankovic, W. J. Akers, L. Song, S. Achilefu, J. A. Margenthaler, M. D. Pashley, L. V. Wang. J. Biomed. Opt., 15, 046010(2010).

[4] T. Liu, M. Sun, N. Feng, Z. Wu, Y. Shen. Chin. Opt. Lett., 13, 091701(2015).

[5] M. Xu, L. V. Wang. Rev. Sci. Instrum., 77, 041101(2006).

[6] H. F. Zhang, K. Maslov, M. Sivaramakrishnan. Appl. Phys. Lett., 90, 053901(2007).

[7] K. Maslov, H. F. Zhang, S. Hu, L. V. Wang. Opt. Lett., 33, 929(2008).

[8] G. He, B. Li, S. Yang. Front. Optoelectron., 8, 122(2015).

[9] Y. Yuan, S. Yang, D. Xing. Appl. Phys. Lett., 100, 023702(2012).

[10] J. Liang, L. Gao, C. Li, L. V. Wang. Opt. Lett., 39, 430(2014).

[11] Z. Wu, M. Sun, Q. Wang, T. Liu, N. Feng, J. Liu, Y. Shen. Chin. Opt. Lett., 12, 121701(2014).