Truncated correlation photoacoustic coherence tomography: An axial resolution enhancement imaging modality.

In this report we present a novel photoacoustic (PA) modality using pulsed chirp excitation at a fixed wavelength and spectral analysis based on frequency-domain (FD) processing. We introduce Truncated Correlation Photoacoustic Coherence Tomography (TC-PACT), a PA axial resolution enhancement methodology, with an application to closely stacked subsurface multilayers of plexiglass coated with a thin layer of graphite acting as surface absorber. The experimental results demonstrated that the SNR and the axial resolution were higher in TC-PACT than in conventional FD photoacoustics. A ns pulsed laser was also used as PA excitation source with axial resolution results found to be commensurate with TC-PACT. However, the separation distance between adjacent absorbers when the latter approached the axial resolution limit of the transducer was only resolved in TC-PACT compared to pulsed laser PA owing to the lower uncertainty in the former modality involved in determining the location of the absorbers.

Introduction

Biomedical photoacoustic imaging (PAI) has been a growing diagnostic modality over the recent two decades, a modality using light-to-thermoelastic-ultrasound-energy conversion to image biomaterial properties embedded in PA transients [1]. Multi-wavelength PAI depth probing ranges from near surface (also accessible to optical imaging) to a few centimeters below, consistent with ultrasound imaging capabilities [2]. In frequency-domain (FD) photoacoustics, a laser emitting at a specific wavelength is modulated in the ultrasound (US) MHz frequency range to detect PA transients originated from a region of interest (ROI) [3,4]. Pulsed PAI imaging systems use wavelength specific pulsed ns lasers to probe from near surface to deeply seated ROI [5]. In those PAI systems, pulse trains repeated at a fixed pulse repetition rate (PRR) are used to irradiate tissue structures and the corresponding PA transients are detected in the time domain (TD) with application specific transducers [6]. PA transients are then reconstructed into planar (2D) or spatial (3D) output images using a variety of image processing techniques [[7], [8], [9]]. In-vivo imaging in small animal studies, in applications in clinical medicine, and in in-vitro applications in biological research have been carried out by probing down to a few centimeters below the surface with high spatial and lateral resolution [[10], [11], [12]]. Even though the transducer specifications and the laser beam size are mainly responsible for the determination of spatial and lateral resolution of PA images, axial resolution is another critical parameter that has not been paid adequate attention to, especially in the case of closely stacked absorbing layers.

PA transients originating in biomedical systems are proportional to the absorbed optical energy of the laser. Amplitude () of a PA transient detected by an ultrasonic transducer is approximated by . Here, Γ is the Grüneisen proportionality constant, is the optical absorption coefficient and F is the light fluence. The Grüneisen proportionality constant, Γ, is defined as β / , where β is the isobaric volume expansion coefficient, c is the speed of sound and is the specific heat [[13], [14], [15]]. Analyzing PA signals from local sources without spatial distortion when they are close to the axial resolution of the PAI system, especially in deep subsurface locations, is a challenging problem. It has been shown [16] that even with the much wider pulsed ns laser bandwidth than modulated CW laser PA, there are limits due to the bipolar shape of the PA transients. Enhancing the axial resolution of PA instrumentation and/or signal generation and processing will further improve clinical detection capabilities of early abnormalities in tissues.

In the present study, a FD PA system using an intensity modulated continuous wave (CW) laser is used to introduce a new signal generation and processing modality, truncated-correlation photoacoustic coherence tomography (TC-PACT), which was found to enhance the signal-to-noise ratio (SNR) of closely stacked absorbing layers and improve axial resolution compared to multi-frequency chirps, to a level similar to pulsed laser PA. The new PA concept is rooted in Truncated Correlation Photothermal Coherence Tomography (TC-PCT), a photothermal imaging modality which has been demonstrated to improve imaging depth and blurring limitations imposed by diffusion waves to overcome subsurface visualization of 3D structures with optimized depth penetration and high axial resolution [[17], [18], [19], [20]]. A combined theoretical and experimental study of PA signals including TC-PACT, conventional FD PA and pulsed laser PA, was undertaken and used to evaluate and compare some critical parameters such as SNR, contrast and axial resolution.

Theoretical foundations of TC-PACT

The TC-PACT modality is a hybrid between pulsed and FD PA signal generation and processing. Excitation signals used in FD PA with sinusoidal and square-wave linear frequency modulated (LFM) chirps, and with variable or fixed duty cycles, have been discussed elsewhere [21,22]. Fig. 1 shows a thin layer of an optically absorbing medium underneath a scattering medium. For LFM signal generation, laser intensity is modulated with an excitation waveform of fixed pulse width t. In Eq. (1), the Heaviside step function, H, is used to denote a rectangular pulsed signal with fixed pulse width

Fig. 1

The 1D model used to generate simulated PA signals. I- is the light fluence. μ is the absorption coefficient. μ - is the scattering coefficient. cs (ρs) and ca (ρa) are the speeds of sound (densities) of the scattering and absorbing medium, respectively. Z, L and l are length range in the scattering and the absorbing medium.

The excitation chirped pulse E(t) is described aswhere n is the number of pulses in the LFM excitation signal, A is the amplitude, and δ(t) is the Dirac impulse function. The starting time, t, of each pulse is given bywherewhere the index i indicates the pulse number in the modulation frequency train and changes from 1 for the starting modulation frequency (f1) to n for the final modulation frequency (f). T is the chirp time duration and B is the frequency bandwidth. The total number of pulses p in a pulsed chirp is given by

Using the Fourier transform of Eq. (1)

The excitation chirp spectrum can be described by the Fourier transform of Eq. (2) in the frequency domain

In the coded excitation technique to provide the depth detectivity, one method is to calculate the cross-correlation function, R(t), of the detected and the reference signals. In the ideal case, the detected signal can be identical to the transmitted signal, if the impulse response of the system is unity. In this case, one needs to calculate the autocorrelation function. Employing the Fourier transformation facilitates the calculation of cross-correlation (or autocorrelation), as well the convolution with the system transfer function. From the Fourier transform of the foregoing expressions we obtain

It should be clarified that the amplitude of the reference excitation waveform used in the cross-correlation is assumed to be unity. In experiments, the excitation amplitude A is controlled by the laser intensity while the reference signal is recorded in the signal processing software. The cross-correlation maximum (peak) of the excitation signal R is thus given by

The SNR of the cross-correlation peak, R, for the pulsed chirp depends on the duration of the pulses that constitute the chirp. In addition to the fact that a rectangular LFM chirp provides better SNR than a sinusoidal chirp, the axial resolution is expected to be further improved by selecting a fixed pulse width to minimize the effects by filtering out the higher harmonic components in R This hypothesis was tested and found to be true using simulated and measured PA transients by processing with LFM and TC-PACT waveforms in the frequency range from 300 kHz to 2 MHz and an applied band-pass filter frequency range from 300 kHz to 6 MHz to process the signal harmonics.

In order to generate transient PA signal simulations, the excitation signal was convolved with the PA response of the target. The back-propagating FD pressure signal from a semi-infinitive absorbing medium is given by [21]where Γ is the Grüneisen parameter, μ is the effective optical attenuation coefficient, μ is the absorption coefficient of the medium, L is the thickness of the scattering medium (between the transducer and absorbing medium), cs (ρs) and ca (ρa) are the speeds of sound (densities) of the scattering and absorbing medium, respectively, is the angular frequency and k = ω/cs is the angular wavenumber of the ultrasonic PA wave back propagating through the scattering medium.

The voltage generated across the transducer is proportional to :where are the sensitivity and transfer function, respectively, of the transducer. The Fourier transform of the cross-correlation (CC) signal then becomes

The PA transient response from a thin layer of thickness l can now be calculated from the inverse Fourier transform of , which is defined as:where,

andHere, k is the wavenumber in the absorbing medium and is the frequency spectrum of the laser intensity. For pulsed chirp excitation . Eq. (13) can be employed for pulsed PA where the PA excitation source is I(t) = δ(t). Considering a very thin layer where the light distribution can be approximated as constant along the thickness l, Eq. (13) can be inverted:Eq. (14) was independently derived by Diebold [23]. In what follows, Eq. (13) will be used to simulate pulsed laser PA transients as well as FD PA responses to validate the experimental results.

TC-PACT signal processing algorithm

Fig. 2(a) shows a schematic flowchart of the TC-PACT reconstruction algorithm [17,18]. It is shown that the excitation pulsed chirp waveform is processed through a delay-incremented signal unit (DIC) with user defined slice width (d) based on the required axial resolution for the system. The range of the slice width (d) can be adjusted up to d ( = 1/ f, where f is the sampling frequency of PA, excitation signal). Output of the DIC is the in-phase (IP) component of the excitation signal. The quadrature (Q) component is computed through the output of a binary exclusive-OR (EX-OR) logic operator based on the input taken from the IP component through frequency doubling of the IP component. Fig. 2(b) shows a portion of IP and Q components of the excitation signals and PA transients of the raw data from the experimental results. Fig. 2(c) shows a portion of cross correlation amplitude computed with the TC-PACT signal processing algorithm from the IP and Q of the excitation signals with PA transients. In the TC-PACT modality, IP and Q components are processed in a manner similar to the enhanced truncated-correlation photothermal coherence tomography (eTC-PCT) processing method [20] to generate the CC amplitude spectrum for each slice width (d) through time gating (d ; i = 1,2,…,n). Furthermore, the IP and Q components of the CC amplitude can be used to highlight the advantage of TC-PACT axial resolution especially for probing closely stacked local absorption peaks. Conventionally, FD-PAI algorithms are processed with the Hilbert transformation (HT) to compute the Q component of the excitation signal and the PA responses throughout the full chirp repetition period. Q from EX-OR gives a single pulse in the TC-PACT modality but Q from HT in the FD-PA modality generates double pulses in positive and negative directions. Hence the CC amplitude magnitude of the Q component in FD-PA generates a broader peak due to these 2 pulses. This is the limiting factor to improving the axial resolution of closely stacked absorbing layers. Thus, a Hilbert transformation is not the best signal processing option to identify localized features embedded in the peak of the TC-PACT CC amplitude. Distortions of the PA signal could appear as a result of such broadening effects and valuable localized spectroscopic information of sample properties might be lost or hidden in the broadened CC amplitude. This potential problem becomes critical with the TC-PACT algorithm which targets the narrowing of the primary CC peak also with the help of better SNR of the IP and Q components of the excitation signal so as to allow otherwise hidden local peaks in the CC amplitude and phase to emerge. Ultimately, the TC-PACT algorithm can be used to reconstruct 3D images using the CC amplitude and phase computed for each slice. In this work, the slice width d = 1/f was selected to obtain the best achievable axial resolution in the TC-PACT CC amplitude which may compromise the high SNR of local absorption peaks in the vicinity of the primary absorption peak in a ROI, the primary feature of the TC-PACT modality. The axial resolution of the peaks can be further improved with the local magnitude of IP and Q channel SNR incorporated in the CC amplitude with FWHM = N/f and integer N < 10. However, (N/ f) may be compromised by the limitations of the much lower axial resolution of the transducer (3.5 MHz) for deep lying ROI. As expected, high TC-PACT axial resolution accompanied by acceptable SNR can be achieved with high frequency transducers with some limitations in the probe depth due to ultrasonic attenuation at high frequencies. Regardless of the choice of transducer, TC-PACT has demonstrated substantially improved FD PA axial resolution when compared to the conventional photoacoustic radar. Therefore, integrating the TC-PACT algorithm to any FD-PAI system with pulsed chirp excitation has the potential to yield higher axial resolution with adequate SNR.

Fig. 2

(a). TC-PACT signal processing flowchart showing how In-Phase (IP) and Quadrature (Q) components are constructed from the excitation LFM chirp signal in order to generate depth resolved cross-correlation (CC) amplitude and phase. (b) and (c). Fig. 2(b) shows a portion of IP and Q components of the excitation signals and PA transients of the raw data from the experimental results. Fig. 2(c) shows a portion of cross correlation amplitude computed with TC-PACT signal processing algorithm from the IP and Q of the excitation signals with PA transients of the raw data from the experimental results.

TC-PACT and pulsed PA experimental system and methods

The TC-PACT concept to improve axial resolution and SNR of closely stacked local absorbers was explored with a FD-PA system. The investigated geometry consisted of two surface-absorbing graphite coated plexiglass (PMMA) plates with dimensions 40 mm x 20 mm x 1 mm. One side of each plate was coated with a thin layer of light absorbing graphite particles. The plate surfaces were pre-conditioned to accommodate the graphite layer by rubbing with commercial sandpaper (mesh 1500). The distance, ds, between the plates was variable, 100 μm < ds < 1200 μm. This plate geometry produced merging PA transients, and TC-PACT was used to explore the maximum achievable axial resolution and SNR of the separation-scanned absorption peaks of the absorbers. The best axial resolution as a first order approximation is given by λ/2 [24], where λ is the wavelength of sound and can be estimated using the speed of sound (c) and the resonance frequency of the transducer, f : λ = c/2f = 0.5 x [1480 (m/s)/3.5(MHz)] = 0.211 mm. The CC amplitude of the local absorbers was calculated to estimate the axial resolution from the experimental PA transients with resonance frequency f of the transducer at 3.5 MHz. The experimental results were subsequently compared with results processed with simulated PA transients.

The plexiglass system and transducer setup were immersed in a water tank to maximize the PA signal. Fig. 3(a) shows the TC-PACT experimental setup. The optical source was an infrared CW laser (Laser Light Solutions, NJ, USA) with built-in laser driver operating at 810 nm and an integrated thermoelectric cooling (TEC) system. The laser irradiation was delivered by an optical fiber to the sample. MATLAB software synthesized excitation chirp waveforms were loaded onto the function generator (Model 33522B, Agilent, CA, USA) and PA transients were captured with a data acquisition card (NI-DAQ, PXIe-5442) using a US transducer, model V382, Olympus Panametrics, USA, 3.5 MHz, connected to a low-noise pre-amplifier (Model 5676, Panametrics). Automated instrument control, data acquisition and signal processing were performed in LabVIEW software environment.

Fig. 3

(a) The TC-PACT experimental setup. A NIR laser beam (810 nm) modulated with a chirp excitation waveform, outfitted with an optical fiber collimator, was used to illuminate two graphite coated Plexiglass (PMMA) plates. By changing the distance, d, between the plates, PA transients at each separation distance were captured with a transducer to explore and compare the axial resolution of the graphite absorption peaks using FD-PA and TC-PACT signal processing algorithms. (b) The PA experimental setup using a pulsed Nd:YAG laser.

Initially, the graphite coated surfaces of the plexiglass plates in water were held in close contact with each other while the modulated laser beam was directed normal to the semitransparent absorbing graphite layers. The upper plate was attached to a linear translational stage (Melles-Griot, USA), while the lower plate was held at a fixed position with a vertical fixture mounted on the optical table. Then, the upper plate was gradually moved away at 20-μm steps until the distance between the plates (ds) reached 1200 μm. PA transients at 40 different d locations were measured. Each location corresponded to a scan ID number and the data were acquired with the DAQ card at 100 MHz sampling frequency (f Due to the very weak absorption of bare plexiglass in the NIR spectral range (810 nm), it was confirmed that the measured PA transients originated in the strongly absorbing graphite coated on the plexiglass surfaces. In order to minimize the effect of reflected ultrasonic waves from the various boundaries, the US transducer was focused on the coated plate surface.

Fig. 3(b) shows the alternative PA experimental setup employing a 5-ns pulsed Nd:YAG laser (Continuum, Santa-Clara, CA, 1064 nm, 400 mJ energy/pulse). The transient response of the same sample configuration as in Fig. 3(a) was measured with a goal to compare the pulsed laser and TC-PACT axial resolution originating from both plates. The Nd:YAG laser irradiation was delivered to the target using high damage threshold guide mirrors. In order to reduce the ablation of the graphite coating, the laser pulses were attenuated using passage through an absorbing cell containing an aqueous black ink solution (not shown in Fig. 3(b)). The strong pulsed PA transients were recorded by a digital phosphor oscilloscope (TEKTRONIX, DPO 7104C) without the use of a preamplifier. Each pulsed PA transient data string was generated by a single shot of the pulsed laser. The oscilloscope was externally triggered by the pulsed laser power supply system, so transient recording was synchronized with the trigger signal. Once again, no PA transients from the plexiglass surfaces without graphite coating were detected.

Results and data analysis

The sequence of connected data analysis is detailed in Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8. In Fig. 4, the absolute magnitude of the CC amplitude (black envelope) with the corresponding IP (blue trace) and Q (red trace) components was analyzed. Fig. 4(a) shows that as long as the distance between the absorbers is much larger than the axial resolution limit of the transducer, the conventional FD-PA CC amplitude signal is good enough to accurately compute the distance between the absorbers. Fig. 4(b) shows the deviation between the set and computed values when the absorbers are stacked closer together and close to the axial resolution limit of the transducer. This is the main focus of data analysis. From this stage on, the IP and Q components of the local absorber CC amplitude became the focus of the study instead of the absolute CC amplitude. Fig. 5(a), shows how the same data used in Fig. 4(b) but processed with the TC-PACT modality improved the accuracy of distance measurement between adjacent absorbers. Fig. 5(b) shows how the distance between the absorbers can be measured up to the axial resolution limit of the transducer. Next, Fig. 4, Fig. 5 (simulated results) were validated with experimental results from conventional FD-PA and TC-PACT shown in Fig. 6, Fig. 7. Fig. 6(a) and (b) show how the IP and Q components of the CC amplitude of the local absorbers A and B at 7 different locations change with distance between the absorbers in the traditional FD-PA modality. Fig. 7(a) and (b) show how the IP and Q components of the CC amplitude of the local absorbers A and B at the same 7 different locations change in TC-PACT with distance between adjacent absorbers. In what follows, a comparison between the results of Fig. 6, Fig. 7 shows how the sensitivity (axial resolution) of the IP and Q components of the CC amplitude improved using TC-PACT.

Fig. 4

The CC amplitude processed with a conventional FD-PA algorithm by using simulated PA transients of two closely stacked absorbers A and B are shown with a separation (a) d = 800 μm and (b) 400 μm from each other. In order to compute the SNR of the local peaks in the ROI, the absolute magnitude of the CC amplitude (black envelope) with the corresponding IP (blue trace) and Q (red trace) components were analyzed. The CC amplitude, IP and Q components started to broaden with decreased SNR and deviations from the expected separation distance between the absorbers increased due to the frequency response limitations of the transducer (3.5 MHz). Figures (a) and (b) show the separation distance distortion at d800 and d400 was shifted to 937 μm and 851 μm, respectively, in simulated PA transients processed with a regular FD-PA algorithm.

Fig. 5

The TC-PACT output associated with the IP and Q components of the CC amplitude with maximum achievable slice width d =1/f, where f = 3.5 MHz, after processing with simulated PA transients generated by two closely stacked absorbers at separation distance (a) d = 400 μm and (b) 200 μm. The IP component of the CC amplitude of the moving local peak at point A with increasing SNR can be tracked as it approaches the fixed absorber peak at point B. The TC-PACT outputs in (a) and (b) show that the SNR of the local absorber at the point A improved even though the deviation from the expected separation distance between the absorbers increased at d = 200 due to the axial resolution limits of the transducer (3.5 MHz).

Fig. 6

Experimental results of the IP and Q components of the CC amplitude of the FD-PA modality with varying Scan ID indices. Scan ID index (i) from 1 to 30 represents the selected range of CC amplitude results processed from PA transients when the separation distance between the adjacent plates was measured in gradual steps as described in the experimental procedures. Location of Scan ID 1 denotes the largest distance between the plates (d (Scan ID 1) > d (Scan ID 5) >….> d (Scan ID 25) > d (Scan ID 30)). (a) CC amplitude components processed with the FD-PA IP component shows that broadening changes due to the moving peak A with constant SNR; (b) Q component broadening between A (1), and A (30) on the delay time axis.

Fig. 7

Experimental results of the IP and Q components of the CC amplitude of the TC-PACT modality with varying Scan ID indexes are shown here. (a) The IP component of the CC amplitude shows that the SNR of the peak doubles between A (1), and A (30) in delay time axis. (b) The SNR of the Q component of the CC amplitude peaks not only doubled but was also immune to excessive broadening between A (1), and A (30) on the delay time axis.

Fig. 8

(a), (b). Experimental results for the CC amplitude, IP and Q components processed with FD-PA and TC-PACT, respectively. All color scales are shown in arbitrary units. Variable distances between the absorbers are indicated by the delay time (μs) vs. scan ID index (Odd for IP and even for Q in (a) and (b), respectively. (c): CC amplitude of the PA transients experimentally measured with a 3.5 MHz transducer with TC-PACT as a function of acoustic delay time between the absorbers. (d): The 2nd derivative of the CC amplitude in (c) as a function of Scan ID index. Axial resolution from the 2nd derivative of the CC amplitude bands for the moving absorber A and the fixed absorber B was further improved with the location of the corresponding IP and Q components of the CC amplitude in A and B bands as shown in Fig. 8(b). The vertical color bars in (a) and (b) represent IP and Q channels of FD-PA and TC-PACT and have the same y axis arbitrary units as shown in Fig. 6, Fig. 7. The vertical color bars in (c) represent the same arbitrary units as the absolute values of local IP, Q components shown in (b). The vertical color bar in (d) represents arbitrary units of the second derivative of (c).

Simulation results using conventional FD-PA

Fig. 4(a) shows the CC amplitude processed with a regular FD-PA algorithm using simulated PA transients of two absorbers with 800 μm separation distance between them. For a realistic study of the SNR of the various CC peaks, 15 dB noise was added to the detected signal before calculating the CC. Also, the transducer transfer function was simulated as a bell shaped (cos2) distribution in the Fourier domain [25]. Then, the CC amplitude and its IP and Q components were analyzed in Fig. 4(a) and (b). These figures further show that the calculated separation distances at d = 800 μm and 400 μm between the peaks A, B (actual absorber locations) measured from the local CC amplitude envelope peaks were distorted at 937 μm and 851 μm, respectively. Fig. 4(b) shows that the onset of merging CC peaks as absorber A moves closer to absorber B further compromises the accuracy of separation distance measurement and the ability to identify the presence of the two adjacent absorbers. Probing the location of these local absorbers experimentally would be even more challenging due to the added background noise of the instrumentation.

Simulation results using TC-PACT

Fig. 5(a) shows the simulated TC-PACT output processed with maximum achievable slice width, d (=1/f; f100 MHz), from two closely stacked absorbers at separation distance d = 400 μm from each other. The corresponding calculated distance in TC-PACT was 358 μm. Compared to the FD-PA results, Fig. 4(b), which identified the local absorber separation at the location of d = 400 μm as being 851 μm, the TC-PACT results processed with the same simulated raw data exhibit a much higher localization accuracy. This difference highlights the severity of absorption peak broadening in the FD-PA modality due to the closely stacked absorber layers and its negative effects on separation measurements. The TC-PACT output in Fig. 5(b) shows that the measurement deviation from the actual separation distance between absorbers A and B at d = 200 increased drastically due to the axial resolution limits of the transducer (3.5 MHz). These results clearly indicate the strong functional dependence between transducer response bandwidth and axial resolution, regardless of the PA signal generation and detection waveforms. In order to explore these axial resolution limiting effects, the behavior of the IP and Q components of the CC amplitude was investigated. These components are highly correlated and undergo subtle changes as functions of the relative distance between closely stacked absorbers. The results of this investigation are shown in Fig. 6, Fig. 7, Fig. 8(a) and (b).

In summary, the simulated PA transient detection using the TC-PACT modality at the minimum achievable slice width, d (=1/f) has been shown to have the potential to probe closely stacked planar absorbers with high axial (depth) resolution and high SNR. The minimum detected distance between adjacent absorption peaks was identified to be 200 μm, limited by the 3.5 MHz ultrasonic transducer.

Experimental results with TC-PACT

Fig. 6 shows the experimental results of the IP and Q magnitude of the FD-PA CC amplitude from varying Scan ID indices. Scan ID index (i) from 1 to 30 represents the selected range of CC amplitude results processed from PA transients when the separation distance between the adjacent plates was measured in gradual steps as described in the experimental procedures. Location of Scan ID 1 denotes the largest distance between the plates (d (Scan ID 1) > d (Scan ID 5) >….> d (Scan ID 25) > d (Scan ID 30)).

Fig. 6(a) shows the IP component of the FD-PA CC amplitude as denoted in the moving peak A with constant signal magnitude. Fig. 6(a) further shows that the FD-PA amplitude peak of the CC IP magnitude broadens with inter-absorber distance. The Scan ID index 30 (Red) to Scan ID index 1 (Blue) delay time corresponds to a time span between 3.5 and 4.75 μs. Even though the broadening with distance is expected, the CC amplitude of the IP magnitude does not change measurably. This shows the relative insensitivity of the IP component to the identification of local adjacent absorbers separated by small distances. Fig. 6(b) shows the FD-PA Q magnitude component of the CC amplitude. Similar to the IP component in Fig. 6(a), the CC amplitude peak of the Q component broadens from the Scan ID index 30 (Red) to Scan ID Index 1 (Blue), that is, from 3.6 μs to 4.4 μs in delay time. Unlike the fixed IP signal magnitude component of the CC amplitude shown in Fig. 6(a), the FD-PA Q component of the CC amplitude improves with decreasing distance between the absorbers.

Fig. 7 shows the experimental results of the TC-PACT IP and Q magnitude of the CC amplitude from varying Scan ID indices. Fig. 7(a) shows the IP component with similarly variable scan ID indexes. The IP magnitude doubles when the distance between the plates changes from Scan ID index 30 (Red; nearest) to Scan ID Index 1 (Blue; farthest) in the delay time range of 3.5–4.75 μs signal time span. This result is in juxtaposition to Fig. 6(a), in which the IP component of the FD-PA CC amplitude does not change with distance. Fig. 7(b) shows the TC-PACT Q component of the CC amplitude under the same conditions as in Fig. 7(a). Not only does the Q signal magnitude double, but also it exhibits reduced broadening: Q components are highly sensitive to the separation distance between the adjacent absorbers and can be used to identify closely stacked absorbers.

Fig. 8(a) is the collection of all 40 distances between the absorbers by combining IP and Q components of the conventional FD-PA CC amplitude. Fig. 6(a) and (b) show only 7 of those distances and explain how the sensitivity of those IP and Q components of the CC amplitude changes with shorter distance between absorbers including the limiting distance close to the axial resolution limit of transducer. The color scale of Fig. 8(a) reflects the same range for all 40 distances of IP and Q magnitudes of the CC amplitude shown in Fig. 6(a, b). The color scale of Fig. 8(a) represents the 1.3–13.0 range with 1 order of magnitude changes to match the 1 order of magnitude changes in Fig. 8(b)).

Specifically, in Fig. 8(a) the distance d between the adjacent absorbers was measured starting at the contact point, d < 100 μm, and measured at 40 locations with an incremental step size of 20 μm. The axial resolution of the 3.5 MHz transducer was 200–300 μm [24]. The Scan ID indices of the CC amplitude IP and Q components in Fig. 8(a) are denoted with consecutive odd index (2i-1) and even index (2i), respectively, where index i changes from 1 to 40. Because of it, scan ID indexes doubled the 40 distance between the absorbers. All IP and Q results of 40 Scan IDs are plotted in detail to identify the location of local absorbers. The CC Amplitude of each Scan ID contains IP and Q components, so that the x-axis in Fig. 8(a) comprises 80 CC amplitude results (40 × 2) with consecutive odd and even Scan ID indices corresponding to CC amplitude IP and Q. Fig. 8(a) further shows that in FD-PA, when the distance d decreases in the Scan ID 40–80 range (d < 500 μm), the CC amplitude IP component of the fixed absorber B (red color band in the 4.4–4.6 μs delay time range) is dominant compared to the Q component (broader yellow color peak in the 4.2–4.4 μs delay time range) of the approaching absorber A. As discussed in Fig. 6(a), (b), by analyzing the rate of change of the IP and Q components of the CC amplitude of the moving absorber, it turns out that A(Q)/B(IP) < 1 is the condition in this region for improving the axial resolution of the FD-PA modality even when both absorbers are located at the maximum achievable detection range of the transducer.

Similarly, Fig. 8(b) is also the collection of all 40 distances between the absorbers by combining the IP and Q components of the TC-PACT CC amplitude. Fig. 7(a) and (b) show only 7 of those distances and explain how the sensitivity of these IP and Q components changes with decreasing distance between converging absorbers as they approach the axial resolution limit of the transducer. The color scale of Fig. 8(b) has the same meaning as that of Fig. 8(a): it represents the same range of IP and Q magnitudes of the CC amplitude shown in Fig. 7(a, b).

Furthermore, in a manner similar to Fig. 8(a), (b) shows the TC-PACT CC amplitude IP and Q components using the same experimental raw data as in Fig. 8(a). The figure shows that when d decreases in the Scan ID 40 and ID 80 range (d < 500 μm), the CC amplitude Q component due to the approaching absorber A (lower red color band in the 4.2–4.4 μs delay time range) increases and becomes comparable (A(Q)/B(IP) ∼ 1) to the signal magnitude of the IP component of the fixed absorber B (upper red color band in the 4.4–4.6 μs delay time range). Compared to the FD-PA broader yellow peak of the approaching absorber A shown in Fig. 8(a), the sharper and narrower red double line band in Fig. 8(b) exhibits superior resolution of the location of both TC-PACT processed absorbers. Even though the axial resolution ranges of FD-PA and TC-PACT in Fig. 8(a) and (b) are similar, the higher TC-PACT signal magnitude of the CC amplitude Q component is the critical parameter that improves the axial resolution for quantitative localization of adjacent absorbers over FD-PA. It is important to note that the improved TC-PACT signal IP and Q components become clearly superior to the FD-PA CC amplitude Q component in the range d < 500 μm. Because of this TC-PACT axial resolution advantage, analyzing the rate of change of the CC amplitude IP and Q components suffices for the identification of the presence and location of adjacent weak absorbers when their separation lies at, or beyond, the limits of the maximum detection range of the axial resolution of the transducer. In conclusion, compared to Fig. 8(a), (b) shows how the location of the local absorbers (red lines) can be easily identified from the IP, Q components of the local absorbers when they converge close to the axial resolution limit of the transducer. As shown in Fig. 4, Fig. 5, IP and Q components are out of phase. Therefore, when the IP and Q components of each location are plotted next to each other in Fig. 8(a) and (b), the IP maximum peak appears next to the Q minimum (dip) and vice versa. These max-min sequence formations appear as vertical strips in the figures.

Even though the SNR of the IP, Q components is relatively high in TC-PACT compared to FD-PA, the magnitude of the CC amplitude is 4 orders of magnitude less than that of FD-PA. This occurs because in TC-PACT, the narrowest possible slice width (1/sampling frequency) from each pulse of the pulsed chirp waveform sequence (∼5% of the pulse width) is cross-correlated with the response signal. On the contrary, in FD-PA, the maximum possible slice widths in all pulses (100 % of the pulse width) of the excitation signal are cross- correlated with the response signal. Hence, the IP and Q components of the FD-PA CC amplitude are broader and of increased magnitude in Fig. 8(a) than the corresponding IP and Q components with narrower peaks in TC-PACT, Fig. 8(b). In view of the fact that Fig. 8 only concerns the absorption peak area for enhancing the identification and resolution of the local absorbers, only one order of magnitude of color changes is displayed in Fig. 8(a) and (b) and all values less than 1 order of magnitude are buried in the dark background color region.

Fig. 8(c) shows the TC-PACT CC amplitude () of the PA transients experimentally measured as a function of distance between the absorbers. In order to match with Scan ID labels of the IP, Q components as shown in Fig. 8(a) and (b), the corresponding identification of CC amplitudes is denoted as ½ x Scan ID on the x-axis of Fig. 8(c) and (d). Therefore, the Scan ID 80 label in Fig. 8(a), (b) and ½ Scan ID label 40 on the x-axes of Fig. 8(c) and (d) represent the same location i.e. when the absorbers are closest to each other (d < 100 μm, defined as “touching”). Decreasing Scan IDs imply that one absorber gradually moves away from the other. Scan ID 20 on the x-axis of Fig. 8(c) shows that the moving absorber A TC-PACT peak starts to deviate from the fixed absorber B peak while the absorber B position variance is buried in the CC amplitude. Fig. 8(c) further shows that, unlike the superior IP and Q component axial resolution in Fig. 8(b), the absolute CC amplitude could not resolve the local absorbers in the region d < 500 μm. An important message of Fig. 8(c) is that it shows the advantage of TC-PACT axial resolution information (double line as shown in Fig. 8(b)) could be lost, if only the absolute CC amplitude is analyzed. Taking into consideration both IP and Q components of the CC amplitude in Fig. 8(c) highlights the advantage of using TC-PACT in that raw signal processing form instead of calculating FWHM a posteriori using CC amplitude peaks.

In order to improve the axial resolution of the closely stacked local absorbers even farther, the second derivative of the CC amplitude, was computed and is shown in Fig. 8(d). The advantage of this operation lies in that the double peak/red band presence in the TC-PACT IP and Q which was highlighted in Fig. 8(b) is now seen to be the primary contrast contributor in the 2nd derivative of Fig. 8(d). In summary, Fig. 8(b–d) are generated from the same raw data after processing with TC-PACT. The two red lines in Fig. 8(b) and bright strips in Fig. 8(d) show the same information in different format with enhanced contrast useful for calculating the distance between adjacent local absorbers as they converge close to the axial resolution limit of the transducer. Nevertheless, Fig. 8(b) is the primary, direct and highly sensitive venue to be used for the resolution and identification of local absorber locations in an unknown ROI. Fig. 8(b) results were also validated at each Scan ID location by Fig. 8(d) derived from the 2nd derivative of Fig. 8(c). The contrast enhancement aspect of Fig. 8(d) clearly amplifies the axial separation between absorbers A and B especially at very close distances on the order of the transducer axial resolution limit.

It should be mentioned that the 2nd derivative of the FD-PA CC amplitude did not reveal the location of local absorbers because of inherent lower resolution than the TC-PACT modality as shown in Fig. 6(a) and (b). Using Fig. 8(d) the axial resolution can now be estimated from the distance between the valley regions to be 318 μm. On the other hand, using Fig. 8(b), the axial resolution can be directly computed from the distance between the local absorbers responsible for the double red band and was found to be 281.2 ± 29.6 μm, a number commensurate with that obtained from Fig. 8(d). The CC FWHM calculated from the 2nd derivative of the peaks in the region of Scan ID # 30–40 of Fig. 8(d) and that from Fig. 8(b) did not change. This is so because in this region where both absorbers were in near contact, both peaks shifted with distance without narrowing the measured gap between them. This might be related to reaching the axial resolution limit of the transducer. Nevertheless, even at the closest separation, the IP and Q peaks are still resolved with a well-defined valley between them. Axial resolution degradation due to background noise interference even within the near-contact region was insignificant in TC-PACT.

Pulsed PA axial resolution

Fig. 9(a) shows simulated pulsed PA signals as a function of distance d between the absorbers. A, B, C, D, E and F correspond to the points where the pulsed (ns) transients crossed the baseline (0). The challenge in determining the actual distance d from the simulated PA transients is the uncertainty in the location of the moving absorber at points A, B, C and D and the fixed absorber at points E or F of the x-axis. Simulated PA transients normalized to the 3.5 MHz transducer impulse response show deviations in the calculated distance at d < 300 μm which is less than the axial resolution limit of the transducer. This effect was already encountered with the FD-PA and TC-PACT measurements shown in Fig. 4, Fig. 5. Fig. 9(b) shows the experimental pulsed laser PA transients with inter-absorber distance d as a parameter. The parts of the transient originating with the moving absorber cross the baseline at points A, B, C and D and those due of the fixed absorber at points E or F of the x-axis. The assigned distances (d) between the absorbers at various locations d correspond to D–F (188.9 μm), C–F (287.7 μm), B–E (351.1 μm) and A–E (454.8 μm). Like the computer simulated results at d = 256.5 μm in Fig. 9(a), the experimental results in Fig. 9(b) show that due to distortion in the shape of the PA transient at point D, the calculated distance d  = 188.9 μm is not a reliable measurement because it is below the axial resolution detection range of the transducer [24]. This represents one of the critical error sources in pulsed PA detection that needs to be resolved in order to improve axial resolution accuracy of closely stacked absorbers. Background noise can be reduced with averaging over repeated PA transients. However, minimizing background noise interference to detect closely stacked multi-layers of absorbers in deep subsurface locations is a challenging task, particularly with low frequency transducers.

Fig. 9

(a) Simulated pulsed (ns) PA signal transient evolution with distance (d) between the two absorbers as a parameter. A, B, C, D, E and F correspond to the points where the transients cross the baseline (0). The assigned distances d between the absorbers at various pre-set nominal locations d, d, d, dD-F (256.5 μm), C—E (299.3 μm), B—E (391.3 μm) and A—E (496.5 μm, respectively) are subject to the 3.5 MHz transducer bandwidth limitation distortions similar to Fig. 4, Fig. 5. (b) Experimental pulsed-laser PA transients directly processed in the time domain as a function of inter-absorber separation distance d. The transients from the moving absorber at points A, B, C and D and those from the fixed absorber at points E and F are shown. Like the computed simulated results at d = 256.5 μm in (a), the experimental results show that due to the distorted shape of the PA transient at point D, the computed distance d = 188.9 μm is not a reliable measurement because it is below the axial resolution detection range of the 3.5 MHz transducer.

Fig. 10(a) shows the experimental PA transients obtained from the pulsed laser and processed in the TD as a function of distance d between the absorbers. The slope change in the valley region in-between absorbers is the critical component to determine whether the PA transient contains multi absorbers or just background noise. Determining the distance between the primary peaks of 476.9 μm and 427.6 μm can be done accurately however, that distance is not the distance d between the absorbers. The accuracy of measuring d depends on the following two factors: 1) Calculation of d between points C–E and D–E; and 2) Consideration of the shape of the PA transient in between the absorbers and the slope of the valley region A–C as opposed to the distorted trace B–D. The decision making process for determining the axial resolution with minimal measurement error might be influenced by the presence of background noise. The axial resolution of the ns pulsed laser system measured with the 3.5 MHz transducer was 287.7. In the ns pulsed laser results, we only highlighted the limits of the axial resolution calculations when the distance between the absorbers became close to the axial resolution limit of the transducer (Fig. 10a).

Fig. 10

(a) Experimental pulsed PA transients as functions of separation distance (d between the two absorbers. The axial resolution of the pulsed laser system was d 287.7 μm. (b) IP and Q results of the TC-PACT CC amplitude: Localized IP, Q components of the local absorbers with high SNR were used to calculate the axial resolution of the experimental system: 281.2 ± 29.6 μm, same as the pulsed laser results.

Fig. 10(b) shows the effects on the TC-PACT CC amplitude of the IP and Q components due to inter-absorber distance d variation when approaching its minimum value. As shown in Fig. 7(a) and (b), the Q component is more sensitive than the IP component. Fig. 10(b) results taken from the TC-PACT double red band region in Fig. 8(b) show that the closest achievable axial resolution ranges for the 3.5 MHz transducer is 281.2 ± 29.6 μm. Although the axial resolution calculated from the TC-PACT and the pulsed ns laser TD modalities are in the same range, the high SNR of the IP and Q components of the TC-PACT CC amplitude have the potential to identify local absorbers in closely stacked regions down to the detection limits of the transducer. In contrast, the TD data analysis showed that the pulsed ns laser system attained its best achievable axial resolution only if the measured PA transients were strong enough to cancel out the background wideband measurement noise. This was not the case with the TC-PACT results because the noise bandwidth was substantially less than that of the pulsed laser and so was the noise level. Furthermore, unlike the pulsed laser PA transient requirement to determine the location of the baseline crossing, TC-PACT has the potential to automate the computational process of calculating the minimum inter-absorber distances more consistently, as well as enhance the axial resolution of weaker absorber PA responses.

Conclusions

In summary, a new PA modality named truncated correlation photoacoustic coherence tomography (TC-PACT) was introduced. An improvement of the axial resolution of photoacoustic detection was introduced in the form of the novel TC-PACT modality, intermediate between frequency- and time-domain laser excitation waveform. TC-PACT was shown to be able to resolve local absorbers up to the maximum achievable detection limit of the transducer more accurately than conventional FD-PA and (ns) pulsed PA. The TC-PACT modality has an axial resolution and distance measurement accuracy advantage over (ns) pulsed PA detection in the tested case where closely stacked subsurface local absorbers probed with a low frequency transducer and located as close as possible to the axial resolution limit of the transducer. The high SNR of IP and Q components of the TC-PACT CC amplitude resulted in narrowing of the FWHM of the CC amplitude peaks so that the axial resolution used to identify the location of IP and Q components of local subsurface absorbers was enhanced compared to FD-PA with a minimum pulse width of the IP and Q components of the excitation signal equal to 1/sampling frequency of the transducer. As a result, the sensitivity of the IP, Q components of the CC amplitude using TC-PACT was shown to be significantly higher than FD-PA despite the full cross-correlation and superior SNR of the CC amplitude of the latter. In a situation of multiple subsurface absorber layers within the probed subsurface volume, conventional FD-PA can be used to target the position of the primary (uppermost) absorbing layer. Then, TC-PACT can be to generate the CC amplitude for each slice based on delay time signals and its ability to enhance the axial resolution with increased N, the number of depth-resolved slices. The maximum slice number N can be determined by the distance between the last two pulses of the PA chirp excitation signal. Instead of manually computing the FWHM of the CC amplitude of local peaks, the IP and Q components of the TC-PACT CC amplitude, Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8, show the distance effects in the actual signal generation waveforms up to the limit of adjacent local absorbers approaching each other commensurate with the axial resolution limit of the transducer.

Finally, the study showed that the axial resolution of an absorber system probed with TC-PACT and that of a pulsed ns laser PA set-up attained the highest achievable value within the response bandwidth of a 3.5 MHz transducer. It was further shown that in a specific situation where there is a need to probe weak PA transients originating from closely stacked low absorbers in a deeply seated ROI, one can directly use pulsed PA transients as long as the background noise is lower than the PA transients. However, the high SNR of the TC-PACT amplitude IP and Q components and its improved axial resolution achieves more accurate measurements of intra-absorber distance in closely stacked weak absorbers near, or at, the limits of the ultrasonic transducer frequency bandwidth. The distance between adjacent absorbers can be accurately computed using pulsed PA detection, as long as the distance between the absorbers is much larger than the axial resolution limit of the transducer.

Author contributions

A.J chose and prepared the samples, conducted the experiments, and interpreted the results along with B.L, K.S and A.M. A.J and K.S created the figures and wrote the manuscript with assistive revisions and inputs from A.M and B.L. B.L provided the simulated data. K.S and A.M developed the TC-PACT modality. A.M provided the resources, contributed to the interpretation of the data, the writing and editing of the manuscript, and conducted the overall supervision of the project.

Declaration of Competing Interest

The authors declare that there are no conflicts of interest.