2.1.

Combination of Fluence Compensation and SAFT Accounting for Spatial Response of Acoustic Antenna

The proposed complementary approach includes previously developed fluence compensation technique accounting for complex illumination geometry¹⁴ and SAFT technique additionally accounting for the spatial response of the acoustic antenna. The ASR was derived from the method proposed in Ref. 19 consisting in calculation of point spread function of acoustic antenna using the k-Wave MATLAB Toolbox software.²⁰

The first step is the acoustic reconstruction by SAFT. In conventional SAFT, the focal point of the ultrasonic transducer is treated as a VPD. The conventional VPD-based SAFT consists in calculation of appropriate signal delays¹ relative to the VPD for the adjacent scan lines and then summation of the delayed signals. We employed SAFT to produce reconstructed image in 2-D realization [Fig. 1(b)] as

Fig. 1

(a) Schematic of AR-PAM system, geometry of the scanning head, and 3-D SAFT and (b) schematic of 2-D SAFT.

However, the signal from the source located apart from the $A$-scan axis ($Z$-axis) will be detected only if the source is located within the ASR map, and this condition should be implemented in the SAFT algorithm.

Similar to Ref. 10, we multiplied each $A$-scan line of a raw $B$-scan to a weighting mask of ASR corresponding to the detector current position. As 2-D SAFT is applied to each $B$-scan, we performed calculations of the ASR map¹⁹ using the k-Wave toolbox in 2-D computational grid to reduce calculation time. The modeling was performed for the spatial–temporal grid (7.75 by 10.7 mm) × ($8.8\mu \mathrm{s}$) with the steps $\mathrm{\Delta }x=\mathrm{\Delta }z=50\mu \mathrm{m}$ and $\mathrm{\Delta }t=10\mathrm{ns}$. The maximum frequency of the signal that can be detected by the antenna was equal to 15 MHz. We modeled OA signals from a point source recorded by the arc with the focal distance $R=6.7\mathrm{mm}$ and the aperture $a=5\mathrm{mm}$, which corresponds to our spherically focused detector [Fig. 1(b)].

In 2-D computational grid for the current position $x_f$ of the detector, we placed a point source at the positions $(x_s,z_s)$ and registered a signal from this source at the detector represented by the arch [Fig. 1(b)]. Figure 2(a) shows impulse responses of acoustic detector from point sources located at different depths below and above the focal positions $(x_f,z_s)$.

Fig. 2

(a) Impulse responses of acoustic detector (different colors represent different source depths) and ASP (black line) calculated as maxima of impulse responses, (b) matrix of ASR for current position of the detector, and (c) schematic of detection cone mask for the spherically focused acoustic detector.

The ASR $(x_f,x_s,z_s)$ matrix was formed from maxima of OA impulse responses obtained from the corresponding points $(x_f,x_s,z_s)$. Therefore, the ASR matrix is formed only as the amplitude matrix. The resultant ASR for the current detector position is shown in Fig. 2(b). We applied modified 2-D SAFT technique with ASR weighting to $B$-scans as

As accounting for ASR in SAFT has been implemented in a matrix form [Eq. (3)], its drawback is the large computational time. To overcome this limitation, accounting for ASR was substituted to combination of two procedures: applying a detection cone binary mask (DCM) corresponding to the geometrical shape of ASR [Fig. 2(c)] and image normalization in accordance with the axial sensitivity profile (ASP) of antenna [Fig. 2(a) (black line)]. Both these characteristics were taken from the simulated ASR [Fig. 2(b)]: the DCM is determined by the angular aperture [$2\alpha$ in Fig. 1(b)] of the acoustic antenna and the waist of ASR in the focal plane, whereas the ASP is the $z$-profile of ASR at the detector axis. Thus, the binary mask (DCM) is an additional part that allows into account only for the most principal points. The decrease in the number of summation points leads to significant reduction in the computation time.

The second step of the reconstruction is the fluence compensation technique consisting in normalization of the reconstructed OA signal in a given point for the value of local fluence in this point. To produce this compensation, each $A$-scan of the resulting 3-D OA image is normalized by one-dimensional (1-D) array of effective fluence distribution derived from full 3-D Monte Carlo simulation with optical parameters close to those of imaged biotissue accounting for the geometry of complex optical illumination [Fig. 1(a)].

To compare conventional 2-D SAFT reconstruction with the proposed methods, we introduced image contrast as

2.2.

Phantom and In-Vivo Studies

In Fig. 1(a), the block diagram of the AR-PAM system for in-vivo imaging²¹ is presented. The geometrical configuration of scanning OA head provides confocal laser illumination and ultrasonic detection [Fig. 1(a)]. The illumination fiber bundle consists of 77 optical fibers with $100\text{-}\mu \mathrm{m}$ core diameter and 0.12 NA: 70 of them are oriented at 37 deg to the $Z$-axis and 7 are used for beam sampling. A custom-made spherically focused PVDF detector with 35-MHz central frequency, 30-MHz bandwidth, $F=6.7\text{-}\mathrm{mm}$ focal distance, and 0.6 NA provides $50\text{-}\mu \mathrm{m}$ lateral resolution in the focal area in phantom experiments.

In all phantom and in-vivo experiments, we used Wedge HB 532 (BrightSolutions, Italy) laser providing 1-ns laser pulses at 532-nm wavelength with the repetition rate of 2 kHz. Each laser pulse (illuminates the medium at a fixed position of OA head along the $X$- and $Y$-axes) produces simultaneously US and OA $A$-scans, respectively.²² As a result of continuous scanning, $A$-scans corresponding to discrete $XY$ positions of the scanning head were acquired by 200-MHz 16-bit analog-to-digital converter Razor16 (GaGe) and formed 3-D OA image.

All model experiments were performed with a phantom consisted of 14 copper wires with diameter of $30\mu \mathrm{m}$ embedded in a 2% water solution of Lipofundin. The OA B-scans of the phantom were obtained by mechanical scanning within a 7-mm area in the $X$-direction with lateral scanning step of $5\mu \mathrm{m}$, resulting in a $B$-scan consisting of 1400 $A$-scans. The phantom images were acquired for focus position $z_f=1.7\mathrm{mm}$ under the surface.

Lateral scanning step of $25\mu \mathrm{m}$ was used for all in-vivo experiments. All the images were obtained for focus position $z_f=1.5\mathrm{mm}$ in the $Z$-direction. After acquisition, the raw data were interpolated to achieve a transversal step of $5\mu \mathrm{m}$ with further high-pass filtration. All in-vivo experiments were conducted with Caucasian skin of human palm of a healthy female volunteer aged 24 from the research team.

2.3.

Monte Carlo Simulation of Optical Fluence

To simulate fluence distribution provided by the ring-shaped illumination of the AR-PAM system, previously developed full 3-D Monte Carlo code that is allowed to account for the exact geometry of the probing beam¹⁴ is used. The simulation produces in-depth fluence distribution at the detector axis. The optical parameters of the phantom are obtained from spectrophotometry measurements with a Analytik Jena Specord 250 Plus (Germany) system and properties of skin in vivo are taken as typical for human skin.^23,24 The input values in Monte Carlo code are shown in Table 1.

Table 1

Optical properties of phantom and skin in vivo at the wavelength of 532 nm.

The simulations of fluence distribution were performed for focus position at 1.5 mm under the object surface for in-vivo experiments and at 1.7 mm under the surface for phantom experiments. The results of Monte Carlo simulated fluence profiles for phantom (red line) and skin in vivo (blue line) studies are shown in Fig. 3(a).

Fig. 3

Monte Carlo simulated in-depth distribution of optical fluence normalized by the number of launched photons: (a) for skin phantom ($z_f=1.7\mathrm{mm}$) (red line) and skin in vivo ($z_f=1.5\mathrm{mm}$) (blue line), (b) for different focus positions in phantom: $z_f=0\mathrm{mm}$ (red line); $z_f=1.7\mathrm{mm}$ (blue line), and (c) for varied phantom optical properties shown in Table 2.

Table 2

Optical properties used for Monte Carlo modeling for phantom measurements.

To demonstrate the comparison of the fluence correction for traditional exponential dependence and the proposed Monte Carlo derived fluence distribution with accounting for the complex illumination geometry, we performed Monte Carlo simulations for two different focus positions: $z_f=0\mathrm{mm}$ [Fig. 3(b) (red line), which provides the exponential fluence decay¹⁴], and $z_f=1.7\mathrm{mm}$ [Fig. 3(b) (blue line)].

The result of the fluence compensation depends on the chosen optical properties of the medium. To evaluate the effect of optical properties variations in the compensated image quality, we performed Monte Carlo simulations for 50% variations of ${\mu }_a$ and ${\mu }_s$ above and below the basic phantom values employed for reconstruction. The considered sets of input optical parameters are shown in Table 2: the basic optical parameters reconstructed from spectrophotometry measurements are shown as bold. The results of Monte Carlo simulations for phantom measurements for these optical properties sets are shown in Fig. 3(c).

3.1.

Phantom Experiments

The developed approach was applied to OA data from the phantom consisting of Lipofundin solution with embedded absorbing wires. As it was mentioned in Sec. 2.1, the developed acoustic reconstruction approach consists of two parts: acoustic and optical reconstructions. The raw OA B-scan of phantom for focus position $z_f=1.7\text{-}\mathrm{mm}$ under surface is shown in Fig. 4(a); the same image after fluence compensation is shown in Fig. 4(b). Note that for comparison all the images of the phantom are normalized to the signal intensity of the second wire. Fluence compensation allows to increase the signal level from deeper layers. Corresponding images after employment of 2-D conventional SAFT technique when the reconstructed image is calculated from the raw one in accordance with Eq. (2) are shown in Figs. 4(c) and 4(d) without and with fluence compensation, respectively. One can see that the application of the traditional SAFT approach allows to diminish the effect of distortions originating from acoustic detection while fluence compensation significantly increases signal level from large depths; however, the main drawback of traditional SAFT approach is the low signal from the focal plane of the system. Accounting for ASR has high potential in overcoming this drawback. The results of application of 2-D SAFT accounting for ASR are shown in Figs. 4(e) and 4(f), respectively. One can see a significant enhancement of the signal level from the focal plane; however, accurate account for ASR is very time-consuming (corresponding processing times are shown in Table 3).

Fig. 4

Two-dimensional reconstruction performed for the phantom with embedded absorbing wires: raw B-scans obtained by OA microscope (a, b); 2-D SAFT (c, d); 2-D SAFT + ASR (e, f); 2-D SAFT + DCM (g, h); 2-D SAFT + DCM + ASP (i, j) with (b, d, f, h, j) and without (a, c, e, g, i) fluence compensation. All bars are 1 mm. Solid and dashed green ellipses represent integration areas $\sigma$ and $\mathrm{\Sigma }$ in Eq. (4), respectively.

Table 3

Processing times, axial, and lateral spatial resolutions for different SAFT modifications applied to compensated OA B-scan of the phantom.

To reduce computation time, accounting for ASR was substituted to combination of DCM determining ASR geometry shape and normalization to ASP. Application of the binary mask [shown in Fig. 2(c)] alone resulting in accounting the contribution of the sources within angular aperture is shown in Figs. 4(g) and 4(h), respectively. This approach allows to significantly increase the signal from the focal point, however, at the expense of unbalanced signal in the out-of-focus region. To reduce this effect, accounting for ASP is used. Combination of these approaches for uncompensated and fluence-compensated images is shown in Figs. 4(i) and 4(j), respectively. The resultant image [Fig. 4(j)] demonstrates balanced signal level for superficial, in-focus, and deeper objects in combination with high spatial resolution simultaneously providing 35-fold speedup in processing time with respect to account for ASR (Table 3). Processing time can be further improved by performing SAFT reconstruction in a Fourier domain.²⁵

To compare the efficiency of the applied methods, we calculated the image contrast by Eq. (4) [the areas $\sigma$ and $\mathrm{\Sigma }$ are shown in Fig. 4(c) and equal to 0.008 and $0.44{\mathrm{mm}}^2$, respectively] for conventional 2-D SAFT reconstruction, for compensation of ASR, and for the proposed acceleration approach accounting for DCM and ASP. For each of the six wires located at different depths varying from 0.5 to 3.5 mm in the $Z$-direction, we calculated image contrast Fig. 5(a), axial [Fig. 5(b)], and lateral [Fig. 5(c)] resolutions. Accounting for ASR provides a twofold contrast increase in focus versus conventional SAFT [Fig. 5(a)]. Application of SAFT with DCM alone significantly increases contrast in the focal plane; however, leads to an unbalanced image contrast. Combination of DCM and ASP leads to smaller differences in image contrast at different depths. SAFT with DCM (black lines) and ASP (magenta lines) improves twice the axial resolution for reconstructed OA images in the focal area in comparison with conventional SAFT (blue lines) [Figs. 5(b)]. The considered approaches do not significantly affect lateral resolution calculated as the full width at half maximum of the image profile over the $X$-axis because the compensated acoustic distortions originating from the spherical detector shape are manifested in other directions. The obtained values of image contrast, axial and spatial resolutions of compensated OA $B$-scan for wires at focal point and at the depth of 3.5 mm under surface are shown in Table 3.

Fig. 5

Comparison of image (a) contrast, (b) axial, and (c) lateral resolutions for OA images reconstructed by 2-D SAFT (blue lines), 2-D SAFT + ASR (red lines), 2-D SAFT + DCM (black lines), and 2-D SAFT + DCM + ASP (magenta lines) with fluence compensation.

Application of conventional SAFT technique [Figs. 4(c) and 4(d)] provides high-contrast reconstruction of the signal from wires above and below the focus; however, signal amplitude in the focal area is reduced. This effect is clearly seen in the $B$-scans projection to the $Z$-axis [Fig. 6(a)]. The proposed approach (SAFT + DCM + ASP) allows to enhance the signals in focus while fluence compensation increases the signals from the depth [Fig. 6(b) (red line)].

Fig. 6

Comparison of maximum values of OA signals from absorbing wires with and without fluence compensation in projection to the $Z$-axis for two reconstruction methods: conventional 2-D SAFT (a) and 2-D SAFT + DCM + ASP (b).

As one can see from Fig. 4(j), the fluence compensation can enhance the contrast of deep wires compared with uncompensated image [Fig. 4(i)]; however, the effect of compensation depends on the chosen optical properties. Therefore, we performed SAFT + DCM + ASP reconstruction with Monte Carlo-simulated fluence distribution for 50% variations of ${\mu }_a$ and ${\mu }_s$ above and below the basic values for phantom (Table 2).

From Fig. 7, one can see that the Monte Carlo-derived fluence distribution allows to improve contrast of deep wires, even if the selected optical parameters differ from basic values; however, a strong increase of ${\mu }_{\mathrm{a}}$ [Fig. 7(d)] or ${\mu }_{\mathrm{s}}$ [Fig. 7(f)] slightly decreases the contrast of wires located close to the surface.

Fig. 7

B-scans of the phantom with embedded absorbing wires before (a) and after 2-D SAFT + DCM + ASP reconstruction with compensation for Monte Carlo-simulated fluence distribution for different input optical parameters (b–f) shown in Table 2. All bars are 1 mm.

In contrast to Ref. 13 where the fluence distribution is considered as an exponential, we calculated it accounting for the complex illumination beam shape. To compare the results of compensation techniques, we performed reconstruction and further compensation both to exponential decay [Fig. 3(b) (red line)] and complex Monte Carlo-simulated fluence [Fig. 3(b) (blue line)]. Figure 8(b) shows the image compensated for exponential fluence distribution [Fig. 3(b) (red line)] that corresponds to the Monte Carlo modeling for $z_f=0\mathrm{mm}$. This figure demonstrates that compensation for exponential dependence significantly improves the contrast of deep wires at depths from 2 to 3.5 mm, which are not distinguished in the uncompensated image [Fig. 8(a)]. However, such compensation significantly decreases the signal from superficial layers resulting in almost complete disappearance of the superficial wires [Fig. 8(b)] as compared with [Fig. 8(a)] in the image. Figure 8(c) shows the image compensated for Monte Carlo-derived fluence distribution calculated for corresponding focus position $z_f=1.7\mathrm{mm}$ [Fig. 8(c)] proving that such compensation provides much better results as compared with compensation for exponential due to better image balance.

Fig. 8

Reconstructed B-scans by 2-D SAFT + DCM + ASP algorithm without (a) and with fluence compensation for exponential fluence distribution (b), and fluence distribution accounting for complex illumination geometry (c). All bars are 1 mm.

3.2.

In-Vivo Experiments

The developed approach of SAFT with account for DCM and ASP and fluence compensation was applied to OA images of human palm in vivo. First, to each B-scan from 3-D matrix (which corresponds to a fixed position in the $Y$-direction), we employed 2-D SAFT reconstruction accounting for DCM and normalizing to ASP. Second, to each $A$-scan from OA 3-D data, we applied Monte Carlo-derived fluence compensation for focus position $z_f=1.5\mathrm{mm}$ [Fig. 3(a) (blue line)]. 3-D reconstructed image was formed from the set of reconstructed $B$-scans as full 3-D reconstruction approach leads to artifacts originating from the object motion during scanning.

In Fig. 9, one can see two raw [Figs. 9(a) and 9(e)] and fluence-compensated [Figs. 9(b) and 9(f)] $B$-scans with corresponding SAFT + DCM + ASP reconstructed [Figs. 9(c), 9(d), 9(g), and 9(h)] $B$-scans for two different $Y$-positions (dashed and solid frames, respectively) without [Figs. 9(c) and 9(g)] and with [Figs. 9(d) and 9(h)] fluence compensation.

Fig. 9

Two-dimensional reconstruction performed for OA B-scans of human palm vasculature: raw B-scans of human palm obtained by OA microscope (a, b, e, and f); results of 2-D SAFT + DCM + ASP (reconstruction time 41 s) (c, d, g, and h) with (right column) and without (left column) fluence compensation. All bars are 1 mm.

Figure 10 shows raw depth-color-encoded en face OA images of human palm vasculature, built in Fiji,²⁶ with $z_f=1.5\mathrm{mm}$ obtained by the AR-PAM system without [Fig. 10(a)] and with [Fig. 10(b)] fluence compensation together with reconstructed (SAFT + DCM + ASP-processed) images without [Fig. 10(c)] and with [Fig. 10(d)] fluence compensation. The images prove enhancement of OA signals from deeper layers.

Fig. 10

Depth-color-encoded AR-PAM images of human palm vasculature with focus position at the depth of 1.5 mm: raw OA 3-D image (a, b) and reconstructed OA 3-D image by 2-D SAFT + DCM + ASP (c, d) (Video 1) uncompensated (a, c), and compensated for Monte Carlo-simulated AR-PAM fluence distribution (b, d). All bars are 1 mm (Video 1, MPEG, 7.1 MB [URL: https://doi.org/10.1117/1.JBO.23.9.091414.1]).