Axial accuracy and signal enhancement in acoustic-resolution photoacoustic microscopy by laser jitter effect correction and pulse energy compensation

: In recent years, photoacoustic imaging has found vast applications in biomedical imaging. Photoacoustic imaging has high optical contrast and high ultrasound resolution allowing deep tissue non-invasive imaging beyond the optical diffusion limit. Q-switched lasers are extensively used in photoacoustic imaging due to the availability of high energy and short laser pulses, which are essential for high-resolution photoacoustic imaging. In most cases, this type of light source suffers from pulse peak-power energy variations and timing jitter noise, resulting in uncertainty in the output power and arrival time of the laser pulses. These problems cause intensity degradation and temporal displacement of generated photoacoustic signals which in turn deteriorate the quality of the acquired photoacoustic images. In this study, we used a high-speed data acquisition system in combination with a fast photodetector and a software-based approach to capture laser pulses precisely in order to reduce the effect of timing jitter and normalization of the photoacoustic signals based on pulse peak-powers simultaneously. In the experiments, maximum axial accuracy enhancement of 14 µm was achieved in maximum-amplitude projected images on XZ and YZ planes with ± 13.5 ns laser timing jitter. Furthermore, photoacoustic signal enhancement of 77% was obtained for 75% laser pulses peak-power stability.


Introduction
The photoacoustic (PA) effect was first discovered by Alexander Graham Bell over a century ago; however, this effect has attracted more attention and had applications in recent decades due to the availability of advanced short pulsed lasers, fast data acquisition (DAQ) devices, and computing capabilities. One of the major applications of using the PA effect is photoacoustic imaging (PAI), in which a light absorbing specimen is illuminated by the short laser pulse. The light absorption results in a rise in temperature leading to rapid thermoelastic expansion and consequently generation of PA waves. The detection of generated waves provides the possibility of locating initial PA wave sources, i.e., the spatial distribution of light absorbers. These absorbers could be microvasculature in tissue. This imaging method enjoys the advantages of high optical contrast as well as low acoustic scattering in soft tissue simultaneously and can be used for deep tissue imaging beyond the optical diffusion limit [1][2][3][4]. PAI can be utilized for structural imaging of vascular networks, functional imaging of the brain, and blood flow hemodynamics of the soft tissue in vivo [5,6]. PAI is divided into a number of categories, such as photoacoustic tomography (PAT), photoacoustic microscopy (PAM), and photoacoustic endoscopy (PAE) systems [7]. In PAM, the 3D images of absorbers can be directly captured by scanning the sample in confined optical illumination and acoustic detection conditions [8] without further reconstruction algorithms as used in PAT systems [9][10][11][12][13]. Based on the type of optical and/or acoustic focusing used, PAM systems are divided into two categories: 1) optical-resolution (OR-PAM) and 2) acoustic-resolution (AR-PAM) systems. In OR-PAM, tightly focused laser light is utilized to generate acoustic waves detected by a focused or unfocused high-frequency, broadband ultrasound transducer (UT). Alternatively, in AR-PAM, weakly focused laser pulses induce acoustic waves detected by a focused UT. Lateral and axial resolutions of both systems depend on the degree of optical or acoustic focus and the bandwidth of UT [14][15][16]. AR-PAM has been also used for many biomedical imaging applications [17][18][19][20].
In PAI systems, light sources are often Q-switched lasers with the pulse-width of several nanoseconds. Since the lasing process initially starts by spontaneous emission, Q-switched lasers have an inevitable so-called timing jitter (TJ) noise. TJ is the variation in the output time of laser pulses after the external triggering. It can be caused by several phenomena such as the temperature fluctuations of the gain medium, response time of the electronic devices, and pump power/spectrum [21,22]. This jitter causes temporal displacement of PA signals, thereby degrading the axial accuracy of images and not allowing signals temporal averaging to increase SNR [23]. A TJ of 10 ns axially shifts the signals by approximately 15 µm (in water or soft tissue) and thus a random pattern is obtained from a smooth sample [24]. When constructing the images of maximum-amplitude projection on XZ or YZ planes (MAPxz, MAPyz), the integration of these random uneven patterns increases the sample image thickness (as illustrated in Fig. 1). In Fig. 1, (a) is an example of imaging set-up for a tube-like sample. Dashed lines on the sample indicate the direction of the scan lines to capture B-mode images which are comprised of A-mode signals. (b) and (d) are 2D B-mode and 3D images of the sample in the absence of TJ, while (c) and (e) refer to the same images with TJ. Since the sample is opaque, measured thickness is the envelope width of the generated A-mode signals and not the sample thickness. Fig. 1. Effect of timing jitter on photoacoustic signals, B-mode, and MAP images. The system scans the sample according to the scan lines marked on the sample. Signals affected by TJ undergo a temporal displacement and thus the constructed B-mode images will have an uneven pattern. The combination of these uneven B-mode images results in a thicker 3D or MAP images.
Although the axial resolution of PAM systems utilizing UTs with the bandwidth of 70-80% is around 30 µm, in case of imaging fine structures such as microvasculature networks, the axial accuracy and visibility of the images are degraded due to the TJ of the lasers, since the diameter and axial distance of the structures are comparable to the axial resolution of the system. Moreover, in very fast scanning OR-PAM systems for imaging fine structures where TJ is in the order of microseconds, TJ could also reduce the accuracy of MAP images on the XY plane (MAPxy) and lateral resolution due to the mismatch between the scanning and laser excitation positions at the XY coordinate. This is only observable in the fast scanning axis [25]. Proposed solutions for TJ problem are discussed later in this section.
Another problem with the Q-switched light sources is the laser pulses peak-power (PP) variations. Since the amplitude of generated PA signals is related to the pulse energies based on the equation p 0 = Γµ a F where p 0 is the initial induced pressure (pa), Γ the dimensionless Grüneisen coefficient, µ a the absorption coefficient (cm −1 ), and F the local light fluence(J/cm 2 ). PP variations can randomly change the amplitude of acquired signals. In the case of quantitative and functional imaging, where multiple wavelengths are used to measure oxygen saturation (SO 2 ), PP variations can have a major effect on the precision of images [26]. Therefore, accurate measurement of pulse PP is vital for normalization of the captured signals and thus improving PA signals strength. Two solutions could be used to reduce the effect of PP variations: 1) utilizing highly stable lasers and 2) measuring PP variations and normalization of signals. The first solution requires expensive lasers while the second one needs fast photodetectors (PD) and data acquisition (DAQ) systems.
As regards the TJ problem, two main solutions have been offered: 1) reducing TJ in the lasers and 2) compensating for the TJ effect in PA signals. Khurgin et al. reduced the TJ of a passively Q-switched laser by composite pumping. Two main limitations with this method were low repetition-rate and a high TJ of about 500 ns [21,27]. Using a combined active and passive Q-switched laser, Wang et al. reduced the TJ of a passively Q-switched laser from 58 µs to less than 400 ns. As for an actively Q-switched laser, they could reduce the TJ to below 200 ns at a repetition-rate of 5 kHz [28,29]. In another study, Cole et al. reduced the TJ of a passively Q-switched laser by direct bleaching of saturable absorber from 241 ns to 20 ns at the repetition-rate of 10 Hz [22]. In comparison to other works, they could reduce TJ to a greater extent, yet the repetition-rate was far beyond what is needed in scanning PAM systems. With regard to reducing the effect of laser TJ on PA signals, several methods have been used. One method was to use an external trigger source to synchronize DAQ and an actively Q-switched laser [30][31][32]. This, however, merely compensates for the delay between laser triggering and pulse generation, and thus there still will be TJ noise due to uncertainty about pulse generation after triggering the laser. In some studies, laser pulses were captured by a PD whose output was used to trigger DAQ and in some studies to estimate PP variation [33][34][35]. In this technique, the clock and sampling frequencies of most DAQ systems were not high enough and reduced the timing accuracy of triggering as well as the accuracy with which PP variation is measured for short pulses [36]. In a super-resolution photoacoustic microscopy, Kim et al. used PD and DAQ with an increased sampling-rate of 500 MS/s to further enhance the timing accuracy of triggering and capturing pulse-to-pulse energy variations [37]. Using PA signals generated from transducer surface for temporal alignment of PA signals was another proposed technique to compensate for the laser TJ effect. In an acoustic radiation force impulse (ARFI) imaging technique for measuring cell mechanics, Kang et al. used the PA signals arrival time generated by the ultrasonic transducer surface to align displacement signals of the sample [38]. Despite its strengths in PAM systems, this method is only applicable to thin and transparent samples in the transmission mode and cannot be used in bulky samples or in vivo experiments. In another study, Zhang et al. used a similar technique in a super-resolution photoacoustic tomography system [39]. They utilized DAQ with the sampling-rate of 40 MS/s in the set-up. Since the low sampling-rate of the DAQ does not allow for accurate measurement of TJ and PP variation, this method cannot be applied to high-resolution PAM systems.
In this study, we proposed and experimentally demonstrated direct and simultaneous measurement of laser pulses' PP variation and their temporal location using a fast PD and custom-built high sampling-rate DAQ in an AR-PAM system. By applying timing jitter effect correction (TJEC) and peak-power compensation (PPC) for all signals prior to constructing 3D PA images, a significant improvement was achieved regarding the axial accuracy and intensity of B-mode images as well as MAPxz and MAPyz images. This method is not limited to AR-PAM and can be applied to other PAI systems, including OR-PAM, PAT, and PAE.

Experimental set-up
The block diagram of the AR-PAM system is shown in Fig. 2(a). An actively Q-switched laser (MSL-AO-532-500µJ, CNI Lasers) with the wavelength of 532 nm, pulse-width of 8 ns, and repetition rate of 1-5 kHz was used as the light source. The laser output was coupled to a multi-mode fiber (M29L05 -Ø600 µm, 0.39 NA, Thorlabs) using a coupling lens (focal length = 50 mm). A small portion of laser light was directed toward a fast PD (Silicon PIN Detector, 818-BB-22, Newport) by a beam splitter in order to achieve TJEC and PPC. In addition, we placed a pair of polarizers in front of the PD to control the energy of the incident light and prevent saturation or damage. PD can be placed before or after coupling lens. The output of multi-mode fiber was focused on the sample by an objective lens (focal length = 50 mm) after passing through an opto-acoustic-combiner comprised of two right angle prisms filled with silicon oil [40]. Silicon oil serves as the index matching material for the laser light and as the reflective interface for the PA waves. Generated PA waves are detected by a UT (V214-BB-RM, 50 MHz, Olympus) attached to the right-angle prisms by epoxy resin (G14250, General Purpose, Thorlabs). An acoustic lens (LC4573, f = -10.0 mm, Ø6 mm, UV Fused Silica, Plano-Concave, Thorlabs) was also glued on the bottom of opto-acoustic-combiner using UV-curing adhesive (UV 30-24, LOXEAL Engineering Adhesives) to collect acoustic signals from a focused area. We employed a water tank as the acoustic coupling media with a thin transparent and flexible membrane at its center between the sample and system. Furthermore, a motorized XYZ stage comprised of fast (X direction) and slow-scanning-stages (Y and Z direction) was utilized for XY scanning of the system in a triangular pattern [41]. The fast-scanning-stage scans the imaging head along the X direction whereas the slow-scanning-stages scan both the water tank and the sample along the Y and Z directions. Acquired PA signals along with the laser pulse signals coming from PD are visualized by an oscilloscope (GDS-2304A, 2GS/s, 4Channel, GWINSTEK) and captured by the DAQ.

Data acquisition
A home-made control circuit includes an 8-bit microcontroller (Atmega2560, microchip) with 16-bit timer/counters (T/C) generating two PWM pulses at 1 kHz is used to simultaneously control the laser and DAQ. The duty-cycle of PWM pulses was 20% and 0.1%, respectively. Figure 2(b) illustrates the control circuit and both PWM pulses. The delay of about 1.2 µs was merely used to compensate for the delay between the laser triggering and pulse generation and not for correcting the TJ effect. A signal from the stage is read by the microcontroller to determine the scan state of the stage (High/Low) and is then converted to two signals one of which still indicating the scan state of the stage and the other indicating the scan direction. The DAQ was a 2-channel, 12-bit resolution, and 1.8 GS/s simultaneous-sampling custom-built card (Xilinx Kintex-7 FPGA). The output of UT was amplified by two cascade amplifiers (ZFL-500-LN, 0.1-500MHz, 24dB, Mini-Circuits) along with laser pulse signals from PD were captured with desired durations and digitized by DAQ. Amplifiers were supplied with batteries to eliminate electronic noises from switching power supplies. Also, it was assured that impedance mismatching and external noises in electronic signal transmission lines did not occur. The signals were saved as a single binary file for the entire scan of the sample on a personal computer (CPU: Intel corei9-9900k, Ram: 32 GB). This binary file contains information about the direction of the scan provided by the microcontroller, as well as PD and PA signals for any XY scan location/pixel. PD signals are used for TJEC and PPC of PA signals at the processing stage. Further, the directional information is utilized for rearranging 1D signals to 3D PA data.

Data processing
A C-Sharp (C#) analysis platform was prepared in the Visual Studio (Microsoft Visual Studio 2015) environment. Parallel and multi-thread processing were utilized in the platform to have good performance in loading and processing the large raw data. At first, the raw data is loaded and the PPC is applied to the data based on the PP of laser pulses. Next, TJEC is executed on PA signals by shifting/aligning signals based on the temporal location of laser pulses. Methods of PPC and TJEC are illustrated in Fig. 3 schematically. Figure 3(a, b) and (c, d) are laser pulses and PA signals for two specific locations/pixels, respectively. Based on the temporal location of the laser pulses and their PP, TJEC by realigning signals and PPC by correcting signals amplitudes were executed on signals. In Fig. 3, dashed lines refer to acquired laser pulses and PA signals, and solid lines illustrate corrected signals. Bandpass filtering to eliminate out-of-range frequencies (noise reduction), Hilbert transformation to calculate PA signals envelope, down-sampling to reduce the size of final images, preparation of 2D B-mode images/slices and MAPs on different planes, as well as color-coding were performed in the main processing stage. Further, 16-bit grayscale images were used in the analysis instead of 8-bit mapped images to maintain vertical voltage resolution of the captured PA and PD signals. For a 25-square mm area with a depth of 3 mm and step size of 10 µm, the size of each binary file was approximately 2.6 GB. The whole processing procedure took 30 to 50 seconds depending on the setting.

System performance
In order to examine the system performance and analysis platform, the PA image of a USAF 1951 resolution test target was captured. Figure 4(a, b) illustrates the photograph of the target taken by a stereo microscope and the constructed MAPxy image, respectively. Figure 4(c) shows the intensity profile, edge-spread-function (ESF), and line-spread-function (LSF) for the edge of the square indicated by the red line in Fig. 4(b). The width of LSF, which is the lateral resolution of the system, was measured to be 48 µm. Figure 4(d) shows the simulated photoacoustic shifted signals and indicates the possibility of distinguishing two separate absorbers within an axial distance of 27.4 µm. By using very small samples such as gold nanoparticles, the measured axial resolution of the system can even be better. The obtained results is deemed satisfactory considering the lateral resolution of such imaging systems [8].

Effect of laser timing jitter and peak-power variation on system resolution
The lateral resolution of an AR-PAM system can be calculated theoretically using the equation R L = 0.72 λ / NA , where λ is the acoustic wavelength and NA the numerical aperture of the acoustic lens [14]. Empirically, the lateral resolution of the system can be measured by imaging a small sample or a sharp edge in the system's focus and calculating ESF as well as LSF. The width of the LSF will be the lateral resolution of the system. Based on the equation and also the fact that temporal displacement of the signals doesn't affect the lateral distribution of the signals due to TJ, the lateral resolution will not be disturbed by TJ of the laser pulses. PP variations could make the profile of ESF and LSF noisy but smoothing or fitting ESF by a Sigmoid-like function and LSF by a Gaussian-like function eliminates the effect of PP variations on lateral resolution measurement.
The axial resolution of the AR-PAM system is calculated by the equation R A = 0.88 c / BW , where c stands for the speed of acoustic waves in the medium and BW for the frequency bandwidth of the transducer [14]. In practice, axial resolution is measured by axial FWHM of the B-mode image of a fine sample, such as a nanoparticle which is the FWHM of its PA signal. Based on the axial resolution equation and due to the fact that the temporal/axial displacement of the signals has no effect on the FWHM of the signals or B-mode images, the axial resolution of the system is the same regardless of TJ. PP variations (similar to the case of lateral resolution) could only decrease PA signals amplitude. Since the amplitude of the signal envelope (as a Gaussian-like function) has no considerable effect on the FWHM of the envelope, the axial resolution is not degraded due to PP variations.
TJ affects the PA image axial accuracy when constructing 3D or MAPxz/MAPyz images from B-mode images. First, for small samples such as microvascular networks, structures with a small distance in depth are combined in the images. Second, as a result of this combination and also the axial displacement of signals in B-mode images, the thickness of samples in the images increases. To show this in practice, tiny laser-jet printer toner particles immobilized on a glass slide were imaged. To immobilize particles on the glass slide, a small amount of toner was solved in 70% (v/v) isopropyl alcohol and distilled water solution and then sprayed on the glass slide. It was left to dry at the room temperature and then airflow of approximately 300°C was applied on the surface of the glass slide using a heat gun. After the preparation process was completed, the sample was submerged in the water tank and imaged at an inclined angle to have particles in different positions close to the focus of the system. Figure 5 illustrates the results. Figure 5(a) shows the MAPxy image of the toner sample, in which the presence of particles in different system focus positions is clear. TJ had no effect on the MAPxy image since it was responsible only for the temporal/axial displacement of the signals. Figure 5(b-d) are the MAPxz images of the sample with TJ. Figure 5(e-g) are the corresponding images following TJEC. In these images, the graphs inside dotted boxes show the intensity profiles of raw and corrected images over dotted lines. Dashed graphs correspond to raw images and solid graphs to corrected images. The analogy of the images and graphs with TJ and corrected ones indicates two points. First, mixing is eliminated and the visibility of images is enhanced. Second, the thickness of particles is decreased due to TJEC

Quantitative investigation of impact of laser TJEC and PPC
For a better and quantitative evaluation of the impact of TJEC and PPC, three sets of experiments with different samples including surgical blade, human black hair, and floppy disk were conducted. In all the experiments, an area of approximately 25 square mm with 3 mm in depth was imaged with the laser repetition-rate of 1 kHz. The fast scanning stage speed was 10 mm/s and the scanning and data acquisition lasted 4.2 minutes for the 10 µm step size. Since the validation of the aforementioned corrections in different system conditions is crucial, three experiments in random laser states (warm-up times and energies) were carried out; thus, TJs and laser stability/PP variations were different. Additionally, samples were imaged in different positions of the system focus. Figure 6-8 illustrate the obtained results.  Figure 6 shows the PA imaging results for the surgical blade in different views and modes. In this experiment, the imaging system was slightly out of focus, and maximum TJ and peak-power stability (PPS) of the laser were ±12.3 ns and 75%, respectively. Histograms of the laser TJ and PP are illustrated in Fig. 6(g) and (i). Laser PPS was described using the equation PPS= (100 -∆I / I c * 100), where ∆I is the FWHM of the laser PP histogram and I c its center intensity in Fig. 6(h). Figure 6(a) and (b) are the B-mode images of raw and TJ effect corrected signals. The conversion of random and uneven pattern to a smooth one is clear in these images after TJEC. In Fig. 6(a, b), images in dotted boxes are the edge of corresponding image surfaces extracted by ImageJ and LabVIEW. Surface profiles were generated from edge images and the deviation from mean profiles was calculated. The standard deviations of the profiles were 0.56 and 0.3 for raw and corrected images, respectively. There was an enhancement of 46% in the smoothness of surfaces. Figure 6(c) and (d) show the MAPxz images of raw and TJ effect corrected signals. By plotting the intensity profile for Fig. 6(c) and (d) on the dashed lines, the width of the MAPxz for the surgical blade sample was measured 122 and 113 µm (Fig. 6(h)). Axial accuracy enhancement of 9 µm in the MAPxz image was obtained after TJEC.  and (f) are the MAPxy images of the TJ effect corrected and the TJ effect corrected plus PP compensated signals. As indicated in Fig. 6(j), plotting the intensity profile on the dashed lines in (e) and (f) results in an increase of 45-113% in the intensity of the image. The horizontal colored ribbon in (j) indicates the PP enhancement factor for any pixel on the dashed lines and the vertical color-bar is its scale-bar. This ribbon indicates that the intensity enhancement is not the same for all signals/pixels and that it depends on the PP of the laser pulse in that specific pixel. One-by-one comparisons of the whole pixels of both MAPxy images in Fig. 6(e) and (f) indicated that the mean and maximum enhancement after PPC were 77% and 180%, respectively. The detailed enhancement factors statistics are summarized in Table 1. The results obtained for the human black hair sample are demonstrated Fig. 7. The hairs were fixed on a heavy metal ring prior to imaging in order to be mechanically stable during scanning. Also, imaging was performed in an out-of-focus system configuration; that is, the out-of-focus state of the system in this experiment was adjusted to a greater extent compared to the previous one. Here, the TJ and PPS of the laser were ±11.4 ns and 90%. The standard deviations of surface profiles were 2.56 and 1.81 for raw and corrected images, respectively, and an enhancement of 29% in the smoothness of B-mode images was obtained. Axial accuracy enhancement of 14 µm in the MAPyz image after TJEC and an increase of 10-48% in the intensity of the pixels following PPC on the dashed lines were obtained. For the whole pixels of the MAPxy images, a mean and maximum enhancement of 30% and 64% were yielded. Further information on the enhancement factors after PPC can be found in Table 2.  Figure 8 shows the results of the experiment using the floppy disk sample. Here, imaging was performed in the system focus, and the TJ and PPS of the laser were ±13.5 ns and 95%, respectively. The standard deviation of surface profiles were 0.53 and 0.26 for raw and corrected images, respectively, and an enhancement of 50% in the smoothness of B-mode images was obtained. Axial accuracy enhancement of 14 µm in the MAPxz image and intensity enhancement of 6-14% were obtained. The mean and maximum intensity enhancement for the whole area of MAPxy images were measured 11% and 23%. In this experiment, signal intensity enhancement was smaller than previous ones due to better stability of the laser. Table 3 summarizes additional information about the intensity enhancement due to PPC.

Conclusion
In this study, the influence of laser timing jitter effect correction (TJEC) and pulse peak-power compensation (PPC) in an acoustic-resolution photoacoustic microscopy (AR-PAM) system were investigated. The effect of TJEC and PPC was verified by imaging different samples in different conditions in terms of timing jitters, peak-powers, and systems focus. It was found that our proposed method improved the axial accuracy and mean intensity of acquired photoacoustic B-mode and MAP images by 14 µm and 77%, respectively. Capturing laser pulses with a high temporal resolution, applying TJEC and PPC, and obtaining high processing performance were among the main challenges. To capture laser pulses precisely, a high-speed custom-build DAQ system and photodetector were employed. To have more effective processing performance and to handle big captured data more efficiently, an analysis platform was prepared using C#, which is a fast, fully object-oriented, and statically-typed compiled language. Moreover, parallel and multi-thread processing was used to accelerate processing procedures. Further efforts could be using CUDA supported graphic cards to have an enhanced parallel GPU computing power as well as using processing power capabilities of FPGAs of DAQs which seems to be vital for the clinical application of the photoacoustic imaging systems. The proposed method can be used in multi-wavelength functional and quantitative imaging of living samples where the laser pulse intensities need to be measured accurately. There are no limitations in the application of the method in OR-PAM, PAT, and PAE imaging systems.