Fluence compensation in raster-scan optoacoustic angiography.

Modern optical imaging techniques demonstrate significant potential for high resolution in vivo angiography. Optoacoustic angiography benefits from higher imaging depth as compared to pure optical modalities. However, strong attenuation of optoacoustic signal with depth provides serious challenges for adequate 3D vessel net mapping, and proper compensation for fluence distribution within biotissues is required. We report on the novel approach allowing to estimate effective in-depth fluence profiles for optoacoustic systems. Calculations are based on Monte Carlo simulation of light transport and account for complex illumination geometry and acoustic detection parameters. The developed fluence compensation algorithm was tested in in vivo angiography of human palm and allowed to overcome significant in-depth attenuation of probing radiation and enhance the contrast of lower dermis plexus while preserving high resolution of upper plexus imaging.

Introduction

Biomedical optoacoustic (OA) or photoacoustic imaging is a hybrid modality combining the benefits of optical contrast and ultrasonic resolution [1]. Raster scan optoacoustic mesoscopy (RSOM) [2], [3] is a rapidly emerging technology [4] which employs confocal optical illumination and acoustic detection [5]. Due to high absorption of blood as compared to surrounding biotissues in visible range OA imaging is extremely efficient for in vivo angiography [6] with spatial resolution below 100 μm at millimeter depths. However, strong scattering and absorption in biotissues result in rapid in-depth attenuation of probing laser radiation and limit the diagnostic depth of RSOM imaging [7].

In order to overcome this disadvantage of RSOM bioimaging, the compensation of fluence distribution is required [8], [9]. Typical approach to fluence correction consists in calculation or measurement of the probing radiation distribution within medium and calculating the ratio of raw optoacoustic signal from the image voxel and probing radiation fluence within this voxel [10] thus providing the actual local absorption coefficient [11].

The optical fluence can be reconstructed from diffuse optical tomography measurements [10]. Obviously, such approach allow to accurately account for optical macroinhomogeneitites within object, however, it requires additional object inspection by an independent imaging modality. Alternativelly, diffuse reflectance spectroscopy [12], and acousto-optical techiques [13] can be employed, providing, however, smaller accuracy in optical fluence distribution.

Model-based correction for optical fluence is efficiently implemented in OA tomography [14], when illumination system is stationary with respect to the probing tissue [10], [15], [16] and a single calculation can provide distribution of the probing intensity in the entire measurement volume. Model-based calculations can be performed either with Monte Carlo technique [13], [17] or using diffuse equation [18]. Obviously, phantom mesurements ususally demonstrate excellent performance of fluence compensation as optical properties of the inspected objects are controlled and well-known, on the contrary to measuring real biological tissues. In this respect, concept of multiple irradiation sensing [19] is a promising approach since it does not require a priory knowledge of tissue optical properties.

Regarding fluence compensation in RSOM imaging techniques, previously quantitative multispectral optoacoustic measurements were performed employing an analytical model based on spectral optoacoustic amplitudes at shallow imaging dephs [20] or in the optical diffuse regime [21].

In the diffusive regime modified Monte Carlo convolution approach [22] has been employed for calculating fluence in axially symmetrical illumination configuration for planar-parallel geometry [15]. Since a scanning head of RSOM system moves along the investigated object, fluence distribution should be calculated or measured for each head position. Alternatively, lateral inhomogeneities rarely distributed within object could be neglected, and a single calculation for planar geometry can be employed for all the positions of scanning optoacoustic head.

In this paper we report on the novel approach to fluence compensation in RSOM phantom measurements and in vivo angiography. Each A-scan of the resulting in vivo 3D OA image is normalized by 1D array of effective fluence distribution derived from full 3D Monte Carlo simulation with optical parameters close to those of inspected biotissue accounting for the geometry of optical illumination and ulstrasonic antenna detection diagram.

Materials and methods

RSOM system and measurements

Fig. 1 presents the schematic of RSOM system for in vivo imaging. The geometrical configuration of scanning optoacoustic head provides confocal laser illumination and ultrasonic detection (Fig. 1). The probing fiber bundle consists of 70 optical fibers with 100 μm core diameter and 0.12 numerical aperture used for light delivery to the object. Spherically focused polyvinylidene fluoride detector (35 MHz central frequency, 30 MHz bandwidth, F = 6.7 mm focal distance, and 0.6 numerical aperture) provided 50 μm lateral resolution in phantom experiments with the same system [23].

Fig. 1

Schematic of RSOM system and geometry of the scanning head.

Lateral movements of OA head across the investigated object (Fig. 1) were controlled by computer-driven XY stages M-664 (PI Micos, Germany). Probing laser pulses were generated by HB532 (BrightSolutions, Italy) laser system at λ = 532 nm with maximum repetition rate of 2 kHz. Radiant exposure at the palm’s surface was kept below 7.5 mJ/cm2 level in accordance to 20 mJ/cm2 maximum value allowed by ANSI Z136.1 standard for laser safety. A-scans corresponding to discrete XY positions of the scanning head were recorded by analog-to-digital converter (ADC) CSE1622 (Gage, USA). In order to minimize errors in fluence compensation originating from the biotissue surface roughness, the boundary of the medium was automatically recognized for each B-scan corresponding to the scanning along the X-axis with fixed Y position. Then each A-scan starting from the determined boundary position was compensated for calculated distribution of probing radiation fluence within tissue. Finally, 3D data sets were processed using the reconstruction algorithm [24].

All phantom experiments were performed in a phantom consisted of 14 copper wires with diameter of 30 μm placed in 2% water solution of Intralipid. RSOM B-scans of wires were obtained by mechanical scanning within 7 mm area in X direction with lateral scanning step of 5 μm, used for all the phantom experiments. According to this each B-scan consists of 1400 A-scans. The images were acquired for two different focus position FD = 2.21 mm and FD = 1.1 mm.

All in vivo experiments were conducted at the researchers’s human palms with caucasian skin and different ages: a 20 year-old male, a 33 year-old male and a 24 year-old female. Arms of healthy human volunteers were fixed in palm-up position and pressed by immersion chamber filled with water. Rectangular hole in the bottom of immersion chamber allowed to provide RSOM scanning within 7.5 mm × 7.5 mm area in XY plane. Lateral scanning step of 25 μm was used for all the in vivo experiments. The images were acquired for three different focus positions (on the top of the human palm, 0.5, and 1.5 mm in depth).

Monte Carlo simulation of fluence distribution

To simulate fluence distribution provided by illumination of RSOM system previously developed Monte Carlo code [25], [26] was adapted. We assumed uniform distribution of incident photons within the infinitely thin ring of radius 4.6 mm corresponding to radial distribution of optical fibers in fiber bundle (Fig. 1). The finite 125 μm thickness of the ring was neglected as it was small compared to millimeter distances of the considered geometry (Fig. 1). Further, it was assumed that illumination system produces a uniformly illuminated circle with radius of 1.283 mm in the focal plane which corresponds to the profilometry of the experimental probing beam.

The Monte Carlo simulations were performed in a medium with a total thickness of 6 mm, the step in Z direction was 7.5 μm, that corresponds to the time step in OA-scan. The number of the launched photons in each Monte Carlo simulations equaled 107.

Since the acquired A-scans correspond to the independent scanning head positions, it is convenient to normalize them by a one dimensional profile of the effective fluence. On the other hand, one should account for detection diagram of the employed focused acoustic antenna, for which collection area strongly depends on the distance from the antenna focal plane. In this respect, it is not correct to average the fluence over the equal areas for different probing depths. To derive the 1D fluence distribution we averaged the fluence at a given depth over the area equal to the cross-section of the antenna detection diagram at this depth. The shape of the averaging area in lateral direction was assumed to be Gaussian:where z is the in-depth distance from the focal plane, r is the radius of the corresponding cross-section, r0 is the waist of the antenna diagram, and λ is the effective detection wavelength. In simulations the parameters were chosen in accordance with characteristics of the employed antenna: r0 = 25 μm, and λ = 62 μm.

Optical properties of phantom and human skin

An important part of fluence compensation is the choice of proper optical properties of the inspected medium. They can be easily obtained in model experiments with calibrated media, however, in in vivo measurements distribution of the optical properties within medium is quite complex and even knowing the exact tissue geometry may produce discrepance due to inaccuracy in choosing particular optical properties value. In this situation averaged optical properties may produce the desirable result. Optical properties of skin in the spectral range of 500–600 nm from two review papers are summarized in Table 1.

Table 1

Optical properties of human skin in the spectral range of 500–600 nm taken from [27], [28].

Tissue type	λ, nm	μa, cm−1	μs, cm−1	μs´, cm−1	g
Dermis	500	0.837	9.224	–	0.715
	550	0.786	7.722	–	0.715
	600	0.694	6.309	–	0.715
Dermis	500	0.336 (0.43)	–	4.62 (4.6)	–
	514	0.31 (0.4)	–	4.32 (4.1)	–
	520	0.30 (0.4)	–	4.20 (3.9)	–
	570	0.22 (0.3)	–	3.50 (3.1)	–
	600	0.172 (0.24)	–	3.22 (2.9)	–
Caucasian bloodless dermis	500	0.345	11.99	–	0.120
	550	0.228	10.81	–	0.288
	600	0.181	9.738	–	0.410
Dermis	514	0.3	25	58	0.77
	585	0.3	19.6	41	0.79
Dermis	517	0.22	21	–	0.787
	585	0.22	20.5	–	0.790
	590	0.22	20	–	0.800
	595	0.22	20	–	0.800
	600	0.22	20	–	0.800

The input parameters for Monte Carlo simulations of fluence distribution for human skin were chosen as averaged values of those summarized in Table 1 optical properties. Basing on these reported values the following optical properties of dermis were selected for calculation: μs = 10 mm−1, μa = 0.3 mm−1, g = 0.8, n = 1.37.

For Monte Carlo simulations in phantom the optical properties of 2% aqueous solution of Intralipid [29] were reconstructed from spectrophotometry measurements with Analytik Jena Specord 250Plus system (Germany). The obtained values at the wavelength of 532 nm are μs = 5.37 mm−1, μa = 0.057 mm−1, g = 0.631.

Results and discussion

Monte Carlo simulation of axial fluence distribution

The simulations of effective fluence distribution were performed for various focus positions (FP) used during in vivo experiments. In order to demonstrate the importance of the introduced features of the developed approach, we compared the resulting 1D profiles of fluence with simplified results corresponding to point beam illumination, flat beam illumination, and situation when fluence averaging is performed over the equal areas at all considered depths. Considered combinations of optical illumination geometry and averaging determined by the acoustic antenna diagram are shown in Fig. 2.

Fig. 2

Considered configurations of the simulation geometry. Arrows show probing optical beam geometry (point beam (a-c), planar wave (d-f), and RSOM illumination (g-i)), dashed line shows acoustic antenna receiving diagram determining averaging area at given depth (circular area with radius of 0.05 mm (antenna diagram waist radius) (a,d,g), with radius of 1.283 mm (beam waist of the RSOM illumination system) (b,e,h) and with variable radius in accordance with formula (1) (Gaussian shape diagram) (c,f,i)).

Fig. 3 demonstrates axial distributions of fluence within model medium for configurations g,h, and i from Fig. 2 for two FPs: focus is on the surface of the medium and focus is 2.5 mm under the surface of the medium. In the case when the optical focus of the scanning system is located on the sample surface the in-depth fluence decay is very close to exponential, which is usually employed for compensation [30]. However, when the focus is within the object, in the superficial layer the fluence distribution drastically differs from the exponential. It is also evident, that avegaring area significantly affects the result when the focus is 2.5 mm under the surface, hence, account for anetenna diagram shape is important in this situation.

Fig. 3

Axial fluence distribution for different focus positions when averaging is performed over the area with radii 0.05 mm (conf. g), 1.283 mm (conf. h), and over variable area in accordance with antenna diagram (conf. i).

The fluence distribution calculated for configuration g (further addressed as RSOM distribution) was also compared to fluence calculated for point-size monodirectional beam (conf. a-c) and for flat distribution of the probing intensity equivalent to probing with a planar wave (conf. d-f). The results aiming to demonstrate the effect of illumination geometry are shown in Fig. 4. Point beam illumination, as expected, produces extremely high values close to the surface, especially when averaging is performed over the small area concentric with the beam axis (conf. a). In general, results for point-beam illumination geometry strongly depend on averaging and significantly differ from RSOM distribution. Fluence distribution for flat illumunation (conf. d-f) does not depend on averaging and is much closer to RSOM distribution, however, it still differs indicating that RSOM distribution employment is critical for proper quantitative correction.

Fig. 4

Axial fluence distributions for point-beam (conf. a–c) and planar wave (conf. d–f) illumination compared to RSOM distribution (conf. i).

Fig. 5 demonstrates RSOM fluence distributions (conf. i) for different in-depth positions of the illumination system focus. As it was already demonstrated earlier, when the focus is close to the surface the dependence is extremely close to exponential decay, however, starting from FP = 1.5 mm the distribution significantly differs from the exponential for smaller depths under the object surface. Additionally, distribution slope at larger depth also differs for different FPs.

Fig. 5

Axial RSOM fluence distribution (conf. i) for different FPs.

Fluence compensation in OA imaging of phantoms

Firstly, the developed approach was applied to measurements of the phantom consisting of Intralipid solution with embedded absorbing wires. The results of application of fluence compensation to phantom images are demonstrated in Fig. 6: optoacoustic raw B-scan of phantom for focus position FD = 2.21 mm is shown in Fig. 6a, while the results of fluence compensation are shown in Fig. 6b–f. From raw B-scans one can clearly see that the signals from superficial wires are high while deep laying wires cannot be distinguished due to strong in-depth signal attenuation (Fig. 6a). However, application of fluence compensation significantly improves images resulting in much higher contrast to deeper laying wires (Fig. 6b).

Fig. 6

B-scans of the phantom obtained by RSOM system with focus positions at the depth of 2.21 mm under the surface: uncompensated image (a) and images compensated for Monte Carlo simulated RSOM fluence distribution for different optical parameters (b–f) correspond to Table 2. All bars are 1 mm.

Table 2

Optical properties used for Monte Carlo modeling for RSOM phantom measurements.

	Fig. 6b	Fig. 6c	Fig. 6d	Fig. 6e	Fig. 6f
μa, mm−1	0.057	0.0285	0.114	0.057	0.057
μs, mm−1	5.37	5.37	5.37	2.685	10.74
g	0.631	0.631	0.631	0.631	0.631

The effect of compensation depends on the chosen optical properties. For the most efficient fluence compensation the native values of tissue optical properties should be used for calculating fluence distribution. In order to demonstrate how the image quality is changing with variations of optical properties we performed fluence compensation with Monte Carlo-derived fluence in-depth distributions for 50% variations of μa and μs above and below the reconstructed values for phantom. The chosen optical parameters are presented in Table 2: the native optical parameters are shown as bold.

As we can see from Fig. 6, the Monte Carlo-derived fluence distributions allows to improve contrast of deep laying wires, even if the selected optical parameters differ from the ones considered as native (Fig. 6c–f), however, variations of μa and μs distorts the image of wires located on the surface. Meanwhile, Monte Carlo derived fluence distribution for the native values of optical properties allows to obtain the contrast simultaneously for superficial and deep laying wires (Fig. 6b). Note that the understated optical properties with respect to the native values produce smaller effect to the reconstructed image in comparison to the exceeding values.

In order to demonstrate the role of focus position in fluence compensation we acquired RSOM images of the phantom (Fig. 7a, d) for two different in-depth positions of the illumination system focus FD = 1.1 mm (Fig. 7a) and FD = 2.21 mm (Fig. 7d). The central column of Fig. 7(b, e) demonstrates images compensated for exponential fluence distribution that corresponds to the modeling for FD = 0 mm when the dependence is close to exponential decay. The Figs. 7b, e demonstrates that compensation for exponential dependence significantly improves images giving much higher contrast to deeper laying wires at depths of 3.5 mm and 4.5 mm (Fig. 7b, e) which are not distinguished in uncompensated images. However, exponential compensation significantly decreases the signal from superficial layers resulting in almost complete disappearance of the superficial wires (Fig. 7e) as compared to (Fig. 7d). The third column of Fig. 7(c, f) demonstrates the compensated images by RSOM fluence distribution calculated for corresponding focus positions FD = 1.1 mm (Fig. 7c) and FD = 2.21 mm (Fig. 7f). As it was previously shown in Fig. 4, fluence distribution for focus positions (FP = 0 mm, FP = 1 mm) are extremely close to exponential, which leads to the similar results of compensation for exponential function (Fig. 7b) and compensation with fluence calculated for corresponding focus depth (Fig. 7c). However, for deeper focus position (Fig. 7f) compensation to RSOM fluence distribution provides much better results compared to compensation for exponential: there are no distortions of the wires located at smaller depths (Fig. 7f).

Fig. 7

B-scans of the phantom obtained by RSOM system with focus positions at the depth of 1.1 mm (a–c) and at the depth of 2.205 mm (d–f) under the surface: uncompensated (left column), compensated for exponential fluence distribution (central column), and compensated for Monte Carlo simulated RSOM fluence distribution (right column). All bars are 1 mm.

Fluence compensation in OA mesoscopy imaging of human skin

Next, the fluence compensation approach was applied to OA images of human palm in vivo. The first column of Fig. 8(a, d, g, j) demonstrates uncorrected depth-color-encoded OA images of human palm vasculature with FP on the tissue surface (a), at the depth of 0.5 mm (d) and at the depth of 1.5 mm (g, j) obtained by the RSOM system. Note, that Fig. 8 a and d are obtained from the same site with different FPs. One can see that in these images the signal from superficial skin layers manifested by blue palette vessel net prevails due to strong in-depth signal attenuation. Note, that for smaller focus depths the entire image depth range is smaller.

Fig. 8

Depth-color-encoded RSOM images of human palm vasculature with focus position on the tissue surface (a–c), at the depth of 0.5 mm (d–f) and at the depth of 1.5 mm (g–l): uncompensated (left column), compensated for exponential fluence distribution (central column), and compensated for Monte Carlo simulated RSOM fluence distribution (right column). All bars are 1 mm.

The second column of Fig. 8(b, e, h, k) demonstrates corresponding images compensated for exponential distribution (Fig. 5, FP = 0 mm). One can see that this compensation significantly improves images giving much higher contrast to deeper laying vessels. However, for deeper focus position of 1.5 mm (Fig. 8h, k) exponential compensation produces noise manifested by white points originating from amplification of noise signal from greater depths after fluence compensation. Additionally, exponential compensation significantly decreases the signal from superficial layers resulting in almost complete disappearance of this upper vasculature as compared to uncompensated images. It should be noted that due to surface roughness small blood vessels displayed in yellow in uncompensated images Fig. 8(a, d, g, j) are actually closer to the surface than the depth reads in the OA image. However fluence compensation in combination with automated boundary detection, similar to approach reported in [30], allowed to correctly display the depth of blood vessels.

The third column of Fig. 8(c, f, i, l) demonstrates the images compensated by RSOM fluence distribution calculated for corresponding illumination system focus position. As it can be seen from Fig. 4, fluence distribution for superficial focus positions (FP = 0 mm, FP = 1 mm) are very close to exponential, which results in images extremely similar to ones obtained for exponential compensation (Fig. 8c, f). However, for deeper focus position (Fig. 8i, l) compensation to RSOM fluence distribution provides much better results compared to compensation for exponential: noise component is significantly suppressed and upper vasculature becomes more visible.

To conclude, we demonstrated that for considered RSOM configuration fluence compensation allows to significantly improve image quality and enhance contrast of deep laying vessels. For focus position near surface fluence distribution has exponential character, however, when focus is embedded deeper into the tissue, the fluence distribution significantly differs from exponential and accurate Monte Carlo calculation accounting for illumination geometry and acoustic antenna configuration is critical. More than that, considered configuration does not employ sharply focused optical beam (probing beam waist is 1.283 mm in radius); in the case of a sharper focused beam accurate account for beam and antenna geometry becomes even more critical, especially for quantitative optoacoustic imaging.

In order to more evidently demonstrate the results of the proposed compensation algorithm we separately built the images of human vessel net at particular depths before and after fluence compensation for the case when the focus position is 1.5 mm. The depth-color-encoded RSOM images of human palm vasculature are shown in the entire depth range (Fig. 9a, b), in the depth range of 0.5–1.5 mm (Fig. 9c, d) and in the depth range of 1.0–1.5 mm (Fig. 9e, f). One can see that when shifting to deeper range the contrast of deeper laying vessels enhances significantly thus providing additional information about the entire vessel net.

Fig. 9

Depth-color-encoded RSOM images of human palm vasculature with focus position at the depth of 1.5 mm at depths from 0 to 2.1 mm (a, b), at depths from 0.5 to 1.5 mm (c, d) and at depths from 1 to 1.5 mm (e, f): uncompensated (left column), compensated for Monte Carlo simulated RSOM fluence distribution (right column). All bars are 1 mm.

Finally, it is worth mentioning an advantage of the presented RSOM angiography approach against optical angiography. One of the rapidly emerging technologies for high-resolution blood vessel mapping is optical coherence tomography angiography (OCT-A) alose referred as optical coherence angiography [31], [32], [33] based on optical coherence tomography technique providing angiographic images with resolution down to units of microns at depths up to 2 mm. However, a drawback of this approach consists in so-called shadowing effect when a shadow appears in the image below blood vessel originating from stronger absorption in blood vessels with respect to neighbouring tissues [34], as the technique is based on detection of backscattered probing radition. In RSOM angiography the illumination is perfomed by a much wider beam as compared to that in OCT-A and more uniform volumetric distribution of probing radiation in combination with acoustic detection allows to avoid this effect.

Conclusions

In this paper we report on development of novel algorithm for compensation for probing fluence distribution within inspected object in RSOM angiography. Compensation for each A-scan is performed for 1D fluence distribution array derived from Monte Carlo simulation. Monte Carlo simulation is performed in full 3D accounting for actual complex geometry of the OA illumination system and geometry of the detecting high numerical aperture acoustic antenna. Simulations were performed for medium mimicking human dermis. Application of fluence compensation allowed to significantly increase the contrast of deep vessels thus providing better visualization of the vessel plexus in the entire measurement volume. We demonstrate that even simple medium geometry provides significant enhancement of the image conrtast at higher depths, however, accounting for local inhomogeneitites within tissue has high potential for further contrast improvement.

Funding

Russian Science Foundation project (14-15-00709-П).

Conflicts of interest

The authors declare that there are no conflicts of interest.