In vivo deconvolution acoustic-resolution photoacoustic microscopy in three dimensions

Acoustic-resolution photoacoustic microscopy (ARPAM) provides a spatial resolution on the order of tens of micrometers, and is becoming an essential tool for imaging fine structures, such as the subcutaneous microvasculature. High lateral resolution of ARPAM is achieved using high numerical aperture (NA) of acoustic transducer; however, the depth of focus and working distance will be deteriorated correspondingly, thus sacrificing the imaging range and accessible depth. The axial resolution of ARPAM is limited by the transducer’s bandwidth. In this work, we develop deconvolution ARPAM (D-ARPAM) in three dimensions that can improve the lateral resolution by 1.8 and 3.7 times and the axial resolution by 1.7 and 2.7 times, depending on the adopted criteria, using a 20-MHz focused transducer without physically increasing its NA and bandwidth. The resolution enhancement in three dimensions by DARPAM is also demonstrated by in vivo imaging of the microvasculature of a chick embryo. The proposed D-ARPAM has potential for biomedical imaging that simultaneously requires high spatial resolution, extended imaging range, and long accessible depth. ©2016 Optical Society of America OCIS codes: (170.5120) Photoacoustic imaging; (170.3880) Medical and biological imaging; (170.6900) Three-dimensional microscopy. References and links 1. L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics 3(9), 503–509 (2009). 2. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). 3. K. Maslov, H. F. Zhang, S. Hu, and L. V. Wang, “Optical-resolution photoacoustic microscopy for in vivo imaging of single capillaries,” Opt. Lett. 33(9), 929–931 (2008). 4. V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods 7(8), 603–614 (2010). 5. S.-L. Chen, Z. Xie, P. L. Carson, X. Wang, and L. J. Guo, “In vivo flow speed measurement of capillaries by photoacoustic correlation spectroscopy,” Opt. Lett. 36(20), 4017–4019 (2011). 6. M. Heijblom, D. Piras, W. Xia, J. C. G. van Hespen, J. M. Klaase, F. M. van den Engh, T. G. van Leeuwen, W. Steenbergen, and S. Manohar, “Visualizing breast cancer using the Twente photoacoustic mammoscope: What do we learn from twelve new patient measurements?” Opt. Express 20(11), 11582–11597 (2012). 7. L. Xi, S. R. Grobmyer, G. Zhou, W. Qian, L. Yang, and H. Jiang, “Molecular photoacoustic tomography of breast cancer using receptor targeted magnetic iron oxide nanoparticles as contrast agents,” J. Biophotonics 7(6), 401–409 (2014). 8. J.-M. Yang, C. Favazza, R. Chen, J. Yao, X. Cai, K. Maslov, Q. Zhou, K. K. Shung, and L. V. Wang, “Simultaneous functional photoacoustic and ultrasonic endoscopy of internal organs in vivo,” Nat. Med. 18(8), 1297–1302 (2012). 9. P. Wang, T. Ma, M. N. Slipchenko, S. Liang, J. Hui, K. K. Shung, S. Roy, M. Sturek, Q. Zhou, Z. Chen, and J.X. Cheng, “High-speed intravascular photoacoustic imaging of lipid-laden atherosclerotic plaque enabled by a 2kHz barium nitrite raman laser,” Sci. Rep. 4, 6889 (2014). 10. K. Maslov, G. Stoica, and L. V. Wang, “In vivo dark-field reflection-mode photoacoustic microscopy,” Opt. Lett. 30(6), 625–627 (2005). 11. H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechnol. 24(7), 848–851 (2006). 12. K. H. Song, E. W. Stein, J. A. Margenthaler, and L. V. Wang, “Noninvasive photoacoustic identification of sentinel lymph nodes containing methylene blue in vivo in a rat model,” J. Biomed. Opt. 13(5), 054033 (2008). #253158 Received 3 Nov 2015; revised 12 Dec 2015; accepted 3 Jan 2016; published 7 Jan 2016 (C) 2016 OSA 1 Feb 2016 | Vol. 7, No. 2 | DOI:10.1364/BOE.7.000369 | BIOMEDICAL OPTICS EXPRESS 369 13. Z. Xie, S.-L. Chen, T. Ling, L. J. Guo, P. L. Carson, and X. Wang, “Pure optical photoacoustic microscopy,” Opt. Express 19(10), 9027–9034 (2011). 14. J. Y. Kim, C. Lee, K. Park, G. Lim, and C. Kim, “Fast optical-resolution photoacoustic microscopy using a 2axis water-proofing MEMS scanner,” Sci. Rep. 5, 7932 (2015). 15. I. M. Braverman, “The cutaneous microcirculation,” J. Investig. Dermatol. Symp. Proc. 5(1), 3–9 (2000). 16. G. E. Koehl, A. Gaumann, and E. K. Geissler, “Intravital microscopy of tumor angiogenesis and regression in the dorsal skin fold chamber: mechanistic insights and preclinical testing of therapeutic strategies,” Clin. Exp. Metastasis 26(4), 329–344 (2009). 17. S.-L. Chen, J. Burnett, D. Sun, X. Wei, Z. Xie, and X. Wang, “Photoacoustic microscopy: a potential new tool for evaluation of angiogenesis inhibitor,” Biomed. Opt. Express 4(11), 2657–2666 (2013). 18. C. Zhang, K. Maslov, J. Yao, and L. V. Wang, “In vivo photoacoustic microscopy with 7.6-μm axial resolution using a commercial 125-MHz ultrasonic transducer,” J. Biomed. Opt. 17(11), 116016 (2012). 19. M.-L. Li, H. E. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Improved in vivo photoacoustic microscopy based on a virtual-detector concept,” Opt. Lett. 31(4), 474–476 (2006). 20. Z. Deng, X. Yang, H. Gong, and Q. Luo, “Two-dimensional synthetic-aperture focusing technique in photoacoustic microscopy,” J. Appl. Phys. 109(10), 104701 (2011). 21. R. Ma, S. Söntges, S. Shoham, V. Ntziachristos, and D. Razansky, “Fast scanning coaxial optoacoustic microscopy,” Biomed. Opt. Express 3(7), 1724–1731 (2012). 22. M. Omar, D. Soliman, J. Gateau, and V. Ntziachristos, “Ultrawideband reflection-mode optoacoustic mesoscopy,” Opt. Lett. 39(13), 3911–3914 (2014). 23. H. Estrada, J. Turner, M. Kneipp, and D. Razansky, “Real-time optoacoustic brain microscopy with hybrid optical and acoustic resolution,” Laser Phys. Lett. 11(4), 045601 (2014). 24. D. Cai, Z. Li, and S.-L. Chen, “Photoacoustic microscopy by scanning mirror-based synthetic aperture focusing technique,” Chin. Opt. Lett. 13(10), 101101 (2015). 25. A. Olivo, L. Rigon, F. Arfelli, G. Cantatore, R. Longo, R. H. Menk, S. Pani, M. Prest, P. Poropat, G. Tromba, E. Vallazza, and E. Castelli, “Experimental evaluation of a simple algorithm to enhance the spatial resolution in scanned radiographic systems,” Med. Phys. 27(11), 2609–2616 (2000). 26. O. Michailovich and A. Tannenbaum, “Blind deconvolution of medical ultrasound images: a parametric inverse filtering approach,” IEEE Trans. Image Process. 16(12), 3005–3019 (2007). 27. C. Yu, C. Zhang, and L. Xie, “A blind deconvolution approach to ultrasound imaging,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control 59(2), 271–280 (2012). 28. T. S. Ralston, D. L. Marks, F. Kamalabadi, and S. A. Boppart, “Deconvolution methods for mitigation of transverse blurring in optical coherence tomography,” IEEE Trans. Image Process. 14(9), 1254–1264 (2005). 29. Y. Liu, Y. Liang, G. Mu, and X. Zhu, “Deconvolution methods for image deblurring in optical coherence tomography,” J. Opt. Soc. Am. A 26(1), 72–77 (2009). 30. Y. Wang, D. Xing, Y. Zeng, and Q. Chen, “Photoacoustic imaging with deconvolution algorithm,” Phys. Med. Biol. 49(14), 3117–3124 (2004). 31. C. Zhang, C. Li, and L. 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). 32. T. Jetzfellner and V. Ntziachristos, “Performance of blind deconvolution in optoacoustic tomography,” J. Innov. Opt. Health Sci. 4(4), 385–393 (2011). 33. J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express 21(6), 7316–7327 (2013). 34. L. Zhu, L. Li, L. Gao, and L. V. Wang, “Multi-view optical resolution photoacoustic microscopy,” Optica 1(4), 217–222 (2014). 35. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. 62(1), 55–59 (1972). 36. L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. 79(6), 745–754 (1974). 37. D. S. C. Biggs and M. Andrews, “Acceleration of iterative image restoration algorithms,” Appl. Opt. 36(8), 1766–1775 (1997). 38. H. C. Yalcin, A. Shekhar, A. A. Rane, and J. T. Butcher, “An ex-ovo chicken embryo culture system suitable for imaging and microsurgery applications,” J. Vis. Exp. 44(44), 2154 (2010).


Introduction
Photoacoustic (PA) imaging is a rapidly emerging imaging modality that combines the advantages of high optical contrast and low acoustic scattering [1,2].In PA imaging, biological tissues are illuminated by a short laser pulse, which induces an instantaneous temperature rise and thermoelastic expansion, and then leads to the emission of acoustic waves, also called PA waves.The generated PA waves are detected by ultrasonic transducers and used to form images which map the optical absorption distribution in biological tissues.PA imaging has been demonstrated in a wide variety of biomedical applications such as structural and functional imaging of microcirculation [3][4][5], study of breast cancer [6,7], and endoscopic and intravascular imaging [8,9].
To achieve high-resolution PA imaging, two types of PA microscopy (PAM) can be implemented using a wideband ultrasonic transducer.In acoustic-resolution PAM (ARPAM), high lateral resolution is enabled by acoustic focusing [10][11][12].To further enhance the lateral resolution, optical-resolution PAM (ORPAM) employs tight optical focusing to realize resolution up to several micrometers or even a sub-micrometer scale [3,13,14].In both types, high axial resolution is achieved by utilizing the wideband response to detect a short acoustic pulse.
High-resolution imaging of subcutaneous blood vessels is critical for the study of abnormal peripheral microcirculation [15] and tumor angiogenesis tracking [16].Although the high lateral resolution of ORPAM is vital for leading biomedical investigations such as imaging of tumor angiogenesis [17] and a melanoma cell [18], optical focusing is ineffective when the photon propagation changes to the optical diffusive regime (~1 mm below human skin), significantly deteriorating high resolution beyond the soft-depth limit [1].
ARPAM shows promise in deep penetration by taking advantage of weak ultrasonic scattering in the optical diffusive regime.High spatial resolution is achieved using an ultrasonic transducer with a high central frequency, a wide bandwidth, and a high numerical aperture (NA) [10,11,19].High-NA transducers have a wider angular coverage, which can be utilized in synthetic aperture focusing technique (SAFT) to improve the lateral resolution for out-of-focus regions in ARPAM systems [19][20][21][22][23][24].However, high resolution still faces challenges in ARPAM.First, although the detection of PA waves at a high central frequency and a wide bandwidth enables high spatial resolution, high-frequency PA waves (≥20 MHz) attenuate sharply in skin, which results in high resolution only in a limited imaging depth when imaging subcutaneous blood vessels.Second, drawbacks are introduced by increasing the acoustic NA.Since the depth of focus (DOF) shortens quickly with the increased NA, the in-focus imaging range is sacrificed considerably using a high-NA transducer.Although the lateral resolution for the out-of-focus region can be improved by SAFT using high-NA transducers, equivalently extending the DOF, the restored lateral resolution is poorer, especially for the region far from the acoustic focus, than the original lateral resolution in the focus region of the transducer [20,24].Besides, the working distance (WD) of a high-NA transducer is short, which may even physically hinder subcutaneous imaging.Therefore, it is of great significance to apply different approaches for improving the spatial resolution without the need of increasing the central frequency, the bandwidth, and the NA of the ARPAM system.
The deconvolution algorithms for enhancing the spatial resolution and/or the signal-tonoise ratio (SNR) have been widely studied in various imaging modalities such as X-ray imaging [25], ultrasound imaging [26,27], and optical coherence tomography [28,29].Deconvolution approaches have recently been investigated in PA imaging in ORPAM and PA computed tomography systems [18,[30][31][32][33][34].Although ARPAM has been demonstrated as a valuable tool for biomedical imaging, there is no study, to our knowledge, on the application of deconvolution methods in scanning-based ARPAM to improve the image quality.Besides, previous exploration of deconvolution PA imaging is limited to two dimensions, yet threedimensional (3D) deconvolution volumetric imaging has not been extended.In practice, employing deconvolution to boost both lateral and axial resolutions can be of great value for in vivo and clinical study.The simple deconvolution such as Wiener filtering [18,25] is generally highly sensitive to the SNR; on the other hand, iterative methods yielding better deconvolved images typically require more computation time, especially for 3D deconvolution imaging.Hence, a feasible approach should be studied and evaluated.
In this work, for the first time, we have examined the deconvolution in ARPAM (D-ARPAM) in three dimensions (3D), where both lateral and axial resolutions are much improved.Compared with the deconvolution in ORPAM, the approach developed in this work can restore high resolution beyond the soft-depth limit, which has greater potential for imaging of subcutaneous microvasculature.Compared with the resolution enhancement by changing to a higher-frequency and/or a higher-NA transducer in ARPAM, the penetration depth in skin by the D-ARPAM is not to be degraded as a result of unchanged central frequency of the transducer.Besides, since the acoustic NA is kept same, the high DOF and long WD inherent in the original ARPAM system are maintained.Imaging of polymer beads and tungsten wires with a 20-MHz transducer demonstrated that the D-ARPAM is capable of improving the lateral resolution and the axial resolution by more than 1.8 and 1.7 times, respectively.Further in vivo imaging of microvasculature in a chick embryo has shown the promise of D-ARPAM and its potential in revealing subcutaneous microvasculature.

ARPAM system
The schematic of our ARPAM imaging system is shown in Fig. 1.A lamp pumped solid state Q-switched laser (LPS-S-532, Changchun New Industries Optoelectronics Tech.Co., China) was employed.The laser was working at a wavelength of 532 nm with pulse duration of 8 ns and a repetition rate of 20 Hz.The laser beam with a diameter of ~7 mm was coupled into a 1.5 mm-core multimode fiber (FT1500UMT, Thorlabs, NJ) via a convex lens.To facilitate in vivo imaging, an imaging head consisting of a 20-MHz ultrasonic transducer (V317-SU, focal length = 13.5 mm, element diameter = 6 mm, Panametrics NDT, MA) and the fiber was constructed, and the head was mounted on a three-axis motorized stage.The laser beam emitted from the fiber was aligned with the focal region of the transducer to maximize the detection sensitivity.For in vivo imaging, the optical fluence deposited on the biological tissue was ~5 mJ/cm 2 , which is under the ANSI safety limit (20 mJ/cm 2 ).The sample was placed under a water tank.The water tank with a bottom window sealed with plastic wrap was used to couple the ultrasonic wave from the sample to the transducer.The acoustic signals were then pre-amplified using a commercial electrical amplifier (ZFL-500LNB-S+, Mini-Circuits) and then recorded by a digitizer (PCIe-9852, ADLINK Technology, Inc.) at a 200-MS/s sampling rate.Volumetric ARPAM images were obtained with mechanical scanning of the imaging head using the motorized stage, which was in constant motion during data acquisition.The lateral resolution of this system is determined by the acoustic focus of the transducer, while the axial resolution is limited by the detection bandwidth of the transducer.

D-ARPAM in 3D
Generally in the field of image processing a blurred image I can be expressed as where * denotes the convolution operator, O is the original image, and PSF is a non-negative blur kernel of size that is small compared to the image size.The PSF is also called point spread function.Assuming the presence of additive noise in practical cases, the relationship between the original image O and the output image I can be represented as where n denotes the random spatial distribution of noise.It is commonly assumed that O and n are uncorrelated.The goal of deconvolution is to optimally restore the original image O from the PSF of the system and I.
In our deconvolution, we used an iterative algorithm based on the Lucy-Richardson (LR) deconvolution with a known system PSF [35,36].The iterative equation to seek the optimal estimation of O was derived from the maximum-likelihood estimate approach and is represented by where k is the number of iteration, J is the restored estimation of O, and PSF'(x,y,z) = PSF(−x, −y,−z).Usually, the initial guess of the J 0 is set as I to start the iteration.Then, the image evaluation J can be found by simply iterating Eq. ( 3) until some stop criteria are satisfied.Two stop criteria were used in our case to quantify the quality of the deconvolved image J.One is the normalized adjacent mean square error (AMSE) of two successively restored images in adjacent iterations, which is defined as The other is the absolute error ratio (AER), which is defined as AMSE was used to evaluate the evolution of fluctuation in successively restored images, while AER was used to assess the consistency between the restored image J and the original image O.
Performing a 3D convolution of two 3D matrices with L × M × N elements takes longer time than performing L-times two-dimensional (2D) convolutions of two 2D matrices with M × N elements and M × N-times one-dimensional (1D) convolutions of two 1D vectors with L elements.In order to improve the speed of the deconvolution in 3D to improve both the lateral and axial resolutions of the system, we did not directly apply the LR iterative algorithm in Eqs.(3)-( 5), which require 3D convolution operations.Instead, we adopted a two-step method in our 3D deconvolution.Step 1: The 2D LR iterative deconvolution was applied to each 2D section image separately along the depth direction, which improves the lateral resolution of the image.Step 2: The processed image after step 1 was further used in step 2, where the 1D LR iterative deconvolution was applied to individual A-line PA signals acquired at the 2D scanning positions.We tested and compared the images obtained by the proposed two-step method and by the one-step 3D deconvolution.We find the deconvolved images have similar quality (shown later), yet the two-step method can save computation time and ease the requirement of memory size.

Resolution of the original ARPAM system
We first measure the original spatial resolutions of the imaging system and then determine the enhanced resolutions after applying deconvolution algorithm to evaluate the improvement by the deconvolution ARPAM in 3D.
Spatial resolution is an important parameter for ARPAM.In order to experimentally measure the lateral resolution of our ARPAM, a 6-μm carbon fiber was imaged with unidirectional B scan.The sample was embedded and fixed with 1% agar.Figure 2(a) plots the PA signal profile as a function of the lateral distance.The full width half maximum (FWHM) resolution was estimated to be 280 μm.On the other hand, the −6 dB axial resolution is provided by the time-resolved ultrasonic detection, and is inversely proportional to the ultrasonic transducer bandwidth.The ARPAM axial resolution was also estimated from the imaged carbon fiber.By studying Hilbert transform of the detected A-line PA signal and using the criteria of the FWHM, the axial resolution of the system was determined as ~76 μm, as shown in Fig. 2(b).A normalized 1D Gaussian function with a FWHM of the measured axial resolution was used as the axial PSF of our original ARPAM system.

Lateral resolution of D-ARPAM
We tested the LR deconvolution to evaluate the lateral resolution of D-ARPAM by imaging a 50-μm red color dyed polystyrene microsphere (1050KR, Phosphorex Inc., MA). Figure 3(a) is the maximum amplitude projection (MAP) image by projecting the 3D PA image onto the lateral plane along the scanning direction.The corresponding D-ARPAM image by applying LR deconvolution algorithm is shown in Fig. 3(b).The PSF used in this computation was a normalized 2D Gaussian function with a FWHM of 280 μm, which was chosen from the measured lateral resolution of the original ARPAM system.The computation was terminated when the AMSE reached its minimum value, which is 7 × 10 −5 in this case.We compare the two numbers of iterations when AMSE and AER achieve their respective minimum values.We find that the two numbers are similar and thus, using the criterion of AMSE is sufficient.Similar observations of this consistency were reported in [33].The 1D lateral profiles from the original ARPAM and D-ARPAM are shown in Fig. 3(c) for comparison.In D-ARPAM, the lateral resolution was improved to 75 μm, ~3.7 times finer than that of the original ARPAM (~280 μm).In practical cases, the absorber may have closely-arranged patterns such as high-density microvasculature and cannot be resolved due to image blur.To better evaluate the capability of D-ARPAM in differentiation of two close objects, we imaged a cross pattern made of two 25-μm tungsten wires.Figures 3(d) and 3(e) show the 2D MAP images before and after the LR deconvolution, respectively.It can be clearly seen that the image by D-ARPAM are improved to finer lateral size, which is closer to the original size of the tungsten wire.The criteria that the two targets are still distinguishable should be determined to quantify the improvement of D-ARPAM.In our case, the value of c/spk was checked at different separations of the two objects, where spk is the smaller of the two peaks in the 1D PA profile, and c is the contrast defined as the difference of the spk and the valley between the two peaks.When the c/spk is larger than 0.5, the two targets are regarded as distinguishable.Figures 3(f) and 3(g) show the 1D lateral profiles along the lines in Figs.3(d) and 3(e).In this comparison, the D-ARPAM was able to differentiate two objects with a separation of 420 μm while the original ARPAM can only differentiate those with a separation of 780 μm.The results suggest a ~1.8 times improvement in the case of two close objects.

Axial resolution of D-ARPAM
Similarly, the axial resolution of D-ARPAM was also calibrated in two ways.We first imaged a single 25-μm tungsten wire to get an A-line PA signal (depth profile).Since the PSF is a non-negative blur kernel and there is a positive constraint for the deconvolved image in the LR deconvolution algorithm, the blurred image I in Eq. ( 3) used for LR deconvolution should also be non-negative.Thus, the Hilbert transform (envelope detection) were applied to the Aline PA signal to obtain a monopolar 1D image, as shown in Fig. 4(a).Figure 4(b) shows the 1D profile before and after the LR deconvolution.Determined by their FWHMs, the axial resolution of the original ARPAM and D-ARPAM was 78 μm and 29 μm, respectively, indicating a ~2.7 times improvement.Second, we imaged a sample consisting of two 25-μm tungsten wires arranged with a small angle between them, as shown in Fig. 4(c), in order to calibrate the capability of D-ARPAM to differentiate two close objects in axial direction.Figure 4(d) shows the original B-mode image scanned along the lateral direction.Applying 1D LR deconvolution to each A-line signal in Fig. 4(d), the axial size of the two tungsten wires was obviously reduced, as shown in Fig. 4(e).The criteria used in Section 3.2 were also adopted to evaluate the improvement of D-ARPAM.By checking some depth profiles at different separations in the axial direction, the minimum differentiable separation of D-ARPAM was 68 μm while that of the original ARPAM was 113 μm, as shown in Figs.4(f) and 4(g), respectively, plotted along the lines in Figs.4(d) and 4(e).This corresponds to a ~1.7 times improvement.The deconvolved lateral and axial resolutions are summarized in Table 1.

Phantom imaging by 3D D-ARPAM
To demonstrate the performance of 3D D-ARPAM, the phantom composed of five 50-μm tungsten wires with different orientations was imaged.As mentioned in Section 2, the twostep method for 3D deconvolution was applied to the original 3D ARPAM image.Each section along the depth direction corresponds to a 7.5-μm-thick layer considering the sampling rate of 200 MS/s and the sound velocity of 1500 m/s in our case.Figure 5(a) is the 2D MAP images on XY plane before and after the 3D deconvolution.The enhanced lateral resolution by D-ARPAM enables clearer identification of the wire pattern.Figure 5(b) shows the 2D MAP images on YZ plane before and after deconvolution.As can be seen, the axial size of tungsten wires is also reduced.We compared the images obtained by the proposed two-step method and by the one-step 3D deconvolution.Similar quality of the deconvolved images by the two methods can be observed in Figs.5(a) and 5(b).Thus, in the following, we present the 3D D-ARPAM results by the two-step method.Figures 5(c) and 5(d) show the volumetric 3D images by the original ARPAM system and by the D-ARPAM, respectively.A clear improvement of spatial resolution in both lateral and axial directions can be seen in Fig. 5(d) after performing the 3D deconvolution.Besides, the peak SNR is improved by 10 dB (from 52 dB to 62 dB).Although extra iterations in the 3D deconvolution can be performed to further reduce the pattern length scale, more undesired artifacts would also be introduced in the deconvolved image, resulting in a less accurate and less reliable result in the restored image [37].

In vivo imaging by 3D D-ARPAM
To demonstrate in vivo imaging by the 3D D-ARPAM, a chick embryo chorioallantoic membrane (CAM) model was employed.Fertilized chick eggs were purchased from a hatchery.For easy access to the chick embryo CAM, an ex ovo chick embryo culture method was used [38].Briefly, eggs were incubated for 72 hours using a commercial incubator.Then, they were taken from the incubator and the embryos were transferred to a sterile container.To have more developed blood vessels on the CAM, the embryos were returned back to the incubator for a few days before imaging.The sample on embryo developmental day 5 was used for imaging.Figures 6(a) and 6(b) show the 2D MAP images of the chick embryo CAM on XY plane and on XZ plane, respectively, by ARPAM and D-ARPAM.As can be seen, the two images show a good match on vessel patterns as well as the branching points.The 3D D-ARPAM image has much finer length scales in both lateral and axial directions than the original ARPAM images.Figures 6(c

Discussion
In this work, we investigate the deconvolution algorithm in 3D.The ARPAM imaging system such as the imaging head and the laser can be further upgraded to facilitate in vivo and functional imaging [10,11] region is poorer, especially for the region far from the acoustic focus, than the original lateral resolution in the focus region of the transducer [20,24].In contrast, the deconvolution algorithm explored in this work is to improve both the lateral and axial resolutions within the focus region to be better than the original spatial resolutions in the focus region of the transducer.
Currently, the LR deconvolution was applied only within the acoustic focus region of the transducer.This is because the lateral PSF for the acoustic focus region is different from that for the out-of-focus region, and therefore, the measured lateral PSF (in Section 3.1) cannot be used for the out-of-focus region.The extension of the deconvolution algorithm in 3D to the out-of-focus region is worthy of future research.
A 20 MHz transducer was employed to test the 3D deconvolution ARPAM, which is sufficient for proof-of-concept experiments.Currently, the resolutions of the deconvolved images (lateral: 75 μm; axial: 29 μm) are not comparable to those of high-resolution ORPAM systems (several micrometers).It will be interesting to compare in vivo biological images acquired by the proposed 3D D-ARPAM system and the high-resolution ORPAM system with comparable resolutions of the two systems while the former has the potential to show better imaging depth than the latter.To have similar high resolutions as the ORPAM system, the 3D D-ARPAM system can employ a high-frequency focused transducer (e.g., 75 MHz).
From Figs. 3 and 4, we observed that the FWHM of the object in the case of the "One individual object" is similar to the FWHMs of the two objects in the case of the "Two close objects."That is, the criterion of FWHM gives similar improvement for one and more-thanone objects.In the case of "Two close objects", if the two peak values in the 1D PA profile are the same, the improvement by the criterion of the c/spk>0.5 (mentioned in Section 3.2) will be similar to that by the criterion of FWHM.The improvement is reduced as the difference of the two peak values increases.Thus, we have less improvement by the criterion of the c/spk>0.5.Using nominally the same two objects, we are able to experimentally calibrate the reduced improvement by the criterion of the c/spk>0.5.
The 25 μm wires were used to calibrate the spatial resolutions of D-ARPAM because their small size is a good approximation of a point source and we can better approach the limit of resolution improvement by the LR deconvolution.If a comparable size of samples to the resolutions for our original ARPAM system is used, we expect the size of the blurred image will be ~2 times larger than the original sample size considering the convolution of the two Gaussian spatial profiles from the PA emission and the acoustic detection.That is, the improvement is ≤ 2 times, which is not the best performance of the D-ARPAM system.
The SNRs of the deconvolved images are better than the original ones as deconvolution will converge the intensity to its original sites, e.g., the aforementioned 10 dB improvement by the deconvolution in 3D in Fig. 5.However, the SNRs seem to be degraded in some deconvolved images such as Figs.6 and 7.This may be explained as follows.Before deconvolution, there are different values (intensity) for different pixels with the PA signals.For simplicity, considering the values a and b (a > b) for two pixels before deconvolution and their values a' and b' after deconvolution, the reason that we observe seemingly degraded SNRs may be because after deconvolution, the value a' is much enhanced than b', making the b'/a' smaller than the b/a.That is, the pixel with the value b' will be dimmer in the deconvolved image than the pixel with the value b in the original image.

Conclusions
In this work, we have developed 3D D-ARPAM using LR deconvolution.The deconvolution of PA imaging was performed in 3D for the first time, to our knowledge, to restore images in both lateral and axial directions.As demonstrated by imaging bead and wire phantoms, D-ARPAM provides 3.7 times and 2.7 times finer lateral and axial resolution, respectively, when considering one individual object than the original ARPAM.In the case of imaging two close objects, D-ARPAM offers 1.8 times and 1.7 times enhancement in the ability to differentiate the minimum separation of the two objects in lateral and axial directions, respectively.In vivo study using the 3D D-ARPAM was also investigated on the chick embryo CAM.We more improvement is possible provided images with a higher SNR.Quantitative evaluation of the improved resolution with enhanced SNR is of interest and importance for future study.The promising results demonstrated in this work suggest that 3D D-ARPAM enables biomedical PA imaging with finer spatial resolution while keeping the system's original DOF and WD.Testing the 3D D-ARPAM using a transducer with a higher central frequency (e.g., 75MHz) to reveal even finer structure is of great interest for future work.

Fig. 2 .
Fig. 2. Calibration of the spatial resolution of the original ARPAM system.(a) Lateral resolution; (b) Axial resolution.

Fig. 3 .
Fig. 3. Lateral resolution of the D-ARPAM.MAP images in lateral scanning direction of the microsphere by the original ARPAM system (a) and by D-ARPAM (b).(c) The 1D profiles along the lines in (a) and (b).MAP images of the cross pattern made of two tungsten wires obtained by the original ARPAM system (d) and by D-ARPAM (e).The minimum differentiable lateral separations of the two close objects by D-ARPAM (f) and by the original ARPAM system (g), which are plotted along the lines f and g, respectively, in (d) and (e).

Fig. 4 .
Fig. 4. Axial resolution of the D-ARPAM.(a) A-line PA signal and its Hilbert transform.(b) The depth profiles by the original ARPAM system and by D-ARPAM.(c) The phantom with two wires arranged with a small angle between them.B-mode images of the phantom in (c) by the original ARPAM system (d) and by D-ARPAM (e).The minimum differentiable axial separations of the two close objects by D-ARPAM (f) and by the original ARPAM system (g), plotted along the line f and g, respectively, in (d) and (e).

Fig. 5 .
Fig. 5. Images of the phantom containing five 50 μm-tungsten wires.(a) MAP images on XY plane by the original ARPAM system (left), by the two-step D-ARPAM (middle), and by the one-step D-ARPAM (right).(b) MAP images on YZ plane by the original ARPAM system (left), by the two-step D-ARPAM (middle), and by the one-step D-ARPAM (right).The 1 mm scale bar is for both (a) and (b).The 3D rendering images of the phantom by the original ARPAM system (c) and by the D-ARPAM (d).
) and6(d)  show the 3D volumetric PAM images obtained by the original ARPAM system and the 3D D-ARPAM, respectively, showing the ability of D-ARPAM to restore blurred vessels.Another chick embryo CAM was imaged and 3D deconvolution was applied.

Figure 7 (
a) shows the 2D MAP images on XY plane by the original ARPAM system and D-ARPAM.The vessel pattern can be more clearly identified by D-ARPAM than by ARPAM because of resolution enhancement.Besides, the background noise, which may result from PA signals of relatively small vessels, is suppressed.Furthermore, a 1D profile along the dashed line shown in Fig. 7(a) is plotted in Fig. 7(b) to manifest the ability of D-ARPAM to distinguish two close vessels.The two vessels at Distance of ~3.5−4 mm in Fig. 7(b), corresponding to the regions denoted by the red arrows in Fig. 7(a), can be better visualized by D-ARPAM.

Fig. 6 .
Fig. 6.In vivo 3D imaging of the chick embryo CAM.(a) MAP images on XY plane by the original ARPAM system (left) and by D-ARPAM (right).(b) The 2D MAP images on XZ plane by the original ARPAM system (left) and by D-ARPAM (right).The 1 mm scale bar is for both (a) and (b).The 3D rendering images by the original ARPAM system (c) and by D-ARPAM (d).3D animations are available as supplementary videos: by the original ARPAM system (Visualization 1) and by D-ARPAM (Visualization 2).

Fig. 7 .
Fig. 7.In vivo imaging of the chick embryo CAM.(a) MAP images on XY plane by the original ARPAM system (left) and by D-ARPAM (right).(b) 1D profile along the dashed line shown in (a).The two vessels at Distance of ~3.5−4 mm in (b) corresponds to the regions denoted by the red arrows in (a).