1. Introduction

Light-matter interactions are mainly correlated with absorption, emission, and scattering phenomena. The attenuation, i.e., the total loss of intensity during the interaction, is related to the absorption and represented by the absorption coefficient µ_a. The scattering, in which the incident light phase and direction are changed due to the light-matter interactions along with the light propagation in the media, is represented by the reduced scattering coefficient µ’_s. Since the absorption and the scattering phenomena are frequency-dependent, their coefficients are wavelength-dependent as well [1]. Most of the optical imaging techniques can acquire high-resolution images of the surface. However, imaging inside turbid medium is challenging even when it does not absorb the light. The phase, together with the direction of light, are lost and imaging becomes more difficult as the light optical length within the media is longer. Despite the challenge, different methods for imaging under scattering conditions have been developed. The field of underwater imaging has made a significant progress, with growing necessity, especially for the growing field of autonomous underwater vehicles [2] mainly used for navy purposes. During the last years, various imaging technologies have been trying to overcome the scattering in underwater images; some are based on Sonar acoustic imaging, a technology that is constantly developing [3]. In parallel, other methods were developed based on light illumination [4], requiring compensation algorithms [5]. Another type of method uses the movement of the object [6]. Recent progression in the last two methods includes machine learning tools [7–9]. Yet, the above techniques are limited too; Sonar is limited to big or stiff objects and far distances [10], and optical methods are limited to short distances [11]. Indeed, the imaging techniques can acquire important data which is also intuitive for understanding. However, in addition to their limitations, imaging techniques do not profess to achieve all the information contained in the lost phase. Thus, some information seems to be lost together with the phase [12] and so there are constant attempts to retrieve the phase in various fields such as biology [13], material science [14], astronomy [15] etc. Although the phase cannot be measured directly, some of the information can be obtained from proper analysis of the medium's response to light, which is defined by the medium’s optical properties. Therefore, a sensing technique, the iterative multi-plane optical property extraction (IMOPE) [16–21] technique, was suggested, aiming to detect changes in the reduced scattering coefficient. The IMOPE is a non-invasive technique for the detection of media scattering. It reconstructs and analyzes the reemitted light phase from the irradiated medium, based on the relation between the medium's scattering properties and the reemitted light phase. For reconstructing the phase, the IMOPE uses a multi-plane version of the iterative Gerchberg-Saxton (GS) algorithm [22]. The GS algorithm is an error reduction algorithm [23] for phase retrieval [24] and image reconstruction [25]. The IMOPE technique does not use the retrieved phase itself but its root mean square (RMS). After the RMS of the reemitted light phase is obtained, it is compared to the theoretical model and allows the extraction of the medium’s µ’_s. As mentioned above, the IMOPE technique was originally developed in the red regime of the electromagnetic (EM) spectrum (633nm $\textrm{in wavelength})$ [18,19,26], for medical applications using small distances. Later on, we have extended the IMOPE to the other edge of the visible spectrum, to the blue regime, applying it for 473nm in wavelength [27]. The red wavelength is mainly targeted to biological applications. The blue wavelength however has higher potential for underwater research since the blue wavelength has the lowest absorption undersea [28,29]. Since underwater research will require significant magnification of the IMOPE technique, we show here the potential of the technique to be modified for different, yet small, magnifications. However, for more significant magnifications future work is required.

In this work we have extended the detection range of µ’_s using the IMOPE technique at the blue regime. This extension allows tailoring the linear detection range of the technique to a desired µ’_s range. First, we show an adaptation of the theoretical model to the blue regime. Then, we change the magnification of the image captured by the detector, showing how different magnifications affect the theoretical model. We also observe the change in the linear technique’s detection range for each magnification. Last, we present the magnification effect on the IMOPE using phantoms with known µ’_s, and prove how the experiments support the theoretical work.

2. Theoretical model

Inside a turbid medium, the radiation’s propagation is commonly described by the Radiative Transport Equation (RTE) [30]. It is an integro-differential equation which concerns the different displacements and directions, for a photon, then calculates the losses and benefits in its energy due to scattering and absorption within the medium. Different solutions for the RTE were suggested over the years [31]. The RTE is a complex equation, which has been suggested with different solutions, and the most common one is the diffusion approximation (DA) [32]. The DA describes the diffusion reflection (DR), i.e., the intensity of the scattered light reflected from the surface, as well as the diffusive transmission. The DA uses the method of images with fixed boundary conditions [33] where the boundary (red dotted line in Fig. 1) can be defined either on the surface [32] or higher [34]. The method of images is used such that the light source is described as a real isotropic light source beneath the surface and its image light source symmetrically above the chosen boundary (black circles in Fig. 1). The real source is located at a depth of 1/µ’_s beneath the surface, and its image is reflected from the chosen boundary. The resulted intensity is represented as a function of ρ, the cylindric coordination of the distance from the center of the source along the surface (green arrow in Fig. 1). This method works well for relatively higher ρ, but is limited at smaller distances where it cannot describe the reflectance accurately [33]. To address this, Piao et al. [34] have published a more accurate model for the DR at short paths of ρ and for a low scattering semi-infinite homogeneous medium. According to Piao’s model, the boundary is located at a distance of 2AD where A=(1+R_eff)(1-R_eff) (with R_eff=0.477 [32]) is the mismatch factor and D=(1/3)/(µ_a+ µ’_s) is the diffusion coefficient. In Piao’s model, the real and image sources are referred as master sources (black circles in Fig. 1), where another slave source is beneath the surface, and its image reflects on the other side of the boundary (gray circles in Fig. 1). This model, for shorter light paths and lower scattering coefficients, is referred as the dual source configuration (when the simpler DA model, that includes the master sources only, is referred as the single source configuration). The IMOPE’s theoretical model combines the intensity that is described by the single and double source configurations, I_Rdual(ρ), with a phase model [19] φ(ρ). This, in order to describe the electromagnetic field for calculating its phase RMS as a function of µ’_s:

The distance ρ represents the distance between the center of the pencil beam illumination and a specific pixel. However, the actual pathlength of the light beneath the surface is longer than ρ and therefore the phase is accumulated along this larger distance. The ratio between the actual optical pathlength and ρ is the differential pathlength factor (DPF):

 

Fig. 1. A schematic illustration of the image method as used for the DA model. The medium (light blue area) is illuminated by the laser beam (blue arrow). The light is scattered in the medium (fading blue half a circle). To model the light propagation, the DA model uses the images method, and refers to the field of each point ρ (green circle) as if it was created by the two master sources (black circles) that are located symmetrically above and beneath a chosen boundary (red dotted line). The DA model does not describe well the area close to the center and therefor, Piao’s model describes the field using two additional sources- real and image slave sources (gray circles) located symmetrically above and beneath the chosen boundary (red dotted line).

Download Full Size | PDF

The DPF calculation refers to the intensity of the real and master source (S), the steady state photon fluence rates (Ψ) at the detector.

Thus, the phase itself is:

Where the multiple scattering area is obtained at ρ>MFP’ (referred as ‘ring’), and the single scattering area is the area closer to the light source (blue arrow in Fig. 1), where ρ<MFP’ (referred as ‘center’). The phase RMS of each ROI is calculated separately.

The basic IMOPE technique uses the dual source configuration (Piao’s model) to describe both intensities at the single and multiple scattering areas. This is due to the improvement that this model suggests for the areas that are close to the illumination point and since for optical magnification of M=1, both ROIs are relatively close to the center. Magnifying the field of view, we found that calculation of the intensities for the multiple scattering area (i.e., the ring) with simply the single source configuration fits better with the experimental results; this is the model that better describes the larger distances. However, for the single scattering area (i.e., the center), the dual source configuration remained.

3. Materials and methods

3.1 IMOPE technique

The IMOPE is a technique developed, originally, to extract the reduced scattering coefficient. It evaluates µ’_s from the phase image reconstructed by the GS algorithm. The basic GS algorithm requires two intensity images taken from two planes of the electromagnetic field. At the entrance plane the intensity is notated as I_en, and at the exit plane as I_ex. The intensity of each plane is captured by a camera. The mathematical expression for the field propagation is the Fresnel transform (FRT) [35,36]. From the intensities and the field propagation equations, the algorithm aims to reconstruct the field's phase that was lost when the intensities were captured. Since the GS algorithm converges to the intensity’s RMS, but not to the global minima [22], a multi-plane improved version was suggested [37,38]. The multi-plane GS technique uses N planes, rather than two planes only, along the propagation axis. Thus, the IMOPE uses N light intensity images taken at N planes along the z-axis (Fig. 2(a), along with the light-blue arrows). The method starts with applying the multi-plane GS algorithm. The result is a phase image at the desired N^th plane which is the reconstructed phase of the reemitted light. The average value of the phase image ${\varphi }$ is then subtracted from the received phase image. As mentioned, the received phase images can be separated to their ROIs, when the border between them is the MFP’ (Eq. (4), and the orange circle in Fig. 2(b)). However, in order to explore materials with unknown µ’_s, we have automized the process which extracts the border between the ROIs, and experimentally confirmed it using phantoms with known µ’_s. By knowing the border location, the RMS of each ROI in the phase image can be calculated and compared with the theoretical model (Fig. 2(c) the blue line is the center ROI and the red line is the outer ring ROI), thus the reduced scattering coefficient, µ’_s, has been extracted.

 

Fig. 2. The IMOPE technique’s process. (a) The experimental setup consists of a laser source (1) is connected to an attenuator and a linear polarizer (2). In front of them sits the sample (3). The detector (4), the second linear polarizer (5) and the lens (6) are sitting on a motorized stage. The stage moves while taking the N intensity images required for the GS algorithm. In order to change the magnification, we change the distances between the detector and the lens and between the lens and the object. These distances are represented by the light blue arrows. (b) An example of a phase image extracted by the GS algorithm from a phantom with µ’_s≅0.75 mm^-1. The phase image is divided by the orange circle to the two ROIs: the center (single scattering area, ρ<1/ µ’_s) and the ring (multiple scattering area, ρ>1/ µ’_s). (c) The RMS of each ROI of the phase image is calculated and compared to the theory graph of phase RMS as a function of µ’_s. The suitability of the two calculated values to the ring (red line) and center (blue line) graphs yields the µ’_s of the measured object.

Download Full Size | PDF

3.2 Optical setup

The experimental setup for light intensity images (Fig. 2(a)) is composed of a laser with a wavelength of λ=473nm, with a beam diameter of 1.2mm and power of 100mW, an attenuator, that transfer 1% of the laser's intensity, polarizers- for optical clearing purposes, and a lens (focal length of f=75mm) in order to focus the illumination beam. A CMOS camera was used for the intensity images acquisition. The lens, polarizer and camera are set on a moving stage with a small angle θ from the laser source, hence the distances between the images are corrected appropriately. In order to study the effect of the magnification M on the IMOPE technique, we have changed the distances between the lens, the camera and the sample (blue arrows in (Fig. 2(a)). The samples are set on a 3-axis micrometer stage to enable their fine-tuning during the experiments [37].

3.3 Agar-based phantoms

Agar-based phantoms played a significant role in this research. The suitability of the phantoms for this purpose is since we are able to control their µ’_s accurately, using a known lipid concentration, while having a relatively low absorption coefficient. According to Mie theory, the size of the unit cell and material type determines the scattering coefficient of a substance; the diameter of the intralipid’s unit cell varies between 25-625nm [39]. The phantoms were first used for the validation of the IMOPE technique following a 473nm [27] laser illumination and the different magnifications. The phantoms are composed of Intralipid (IL) (Intralipid 20% Emulsion, Sigma-Aldrich, Israel) with varying concentrations and 1% Agarose powder (Agarose- low gelling temperature, Sigma-Aldrich, Israel) and water [21]. The µ’_s varies as function of the intralipid concentration [1,21]. The IL concentrations of the phantoms used in the following experiments are 0.23%, 0.495%, 0.71%, 1.06%, 1.27%, 1.48%, 1.85%, 2.13%, 2.39%, 2.62% respective to µ’_s values of 0.25, 0.5, 0.7, 1.03, 1.23, 1.4, 1.77, 2.03, 2.28, 2.49 mm^-1 [40].

4. Results and discussion

The theoretical model was first applied to the blue regime, using different magnifications (M: 1,$\frac{2}{3}\; ,\frac{1}{2}$, 0.4 and $\frac{1}{3}$ represented in Fig. 3(a), by red, yellow, green, blue and purple, respectively). For each magnification, the solid lines describe the center ROI, and the dashed lines describe the ring ROI. Our results suggested that different magnifications changed the theoretical graphs, as for a lower M the graphs presented steeper slope and reached saturation earlier. Also, the ring graphs start for lower µ’_s values. For each magnification M, the curves of the ring and the center have a crossing point around the phase RMS value of 1; We noticed that the crossing point of the ring’s graph and the center’s graph shifts with M. Extracting the µ’_s value in which this crossing occurs, µ’_s,cp, for each M (asterisks in Fig. 3(b)) yields the linear fit (line in Fig. 3(b)) µ’_s,cp=M.

 

Fig. 3. The results of the theoretical model’s adaptations. (a) The theoretical model for the different magnifications $M$: $1,\frac{2}{3}\; ,\frac{1}{2}\; ,\; 0.4\; $ and $\frac{1}{3}$ are represented by red, yellow, green, blue and purple, respectively. For each magnification, the solid lines describe the center (i.e, the single scattering ROI) and the dashed lines describe the ring (multiple scattering ROI) . (b) With the change in the the magnification, we have noticed that the crossing point of the ROI’s RMS graphs shifts along the µ’_s axis. The black asterisks show the extracted µ’_s,cp, the µ’_s value in which the crossing point occurs for every magnification M and presented them VS M. The error bars are defined according to the simulation resolution.

Download Full Size | PDF

Next, we applied the measurements on the phantoms using the basic M=1 setup. i.e., for each phantom N intensity images were taken using the optical setup described above and their phase RMS was calculated. Since the µ’_s values of the phantoms were known, we could calculate a theoretical radius that discriminates between the ROIs of each phase image obtained for every phantom; This border is the phantoms’ MFP’ (Eq. (4)). We aimed to find the range of MFP’ which can be detected by our method. Therefore, we used an algorithm that scans the phase image and finds the border radius of each phantom. This algorithm is useful while measuring samples with unknown µ’_s, here however, we used it to examine the capabilities of the IMOPE.

For M=1, the extracted border radius values (the asterisks in Fig. 4(a)) do not match the theory (line in Fig. 4(a)) for the lower µ’_s values. As expected, the IMOPE cannot detect the radiuses of the lower µ’_s since they are larger than the size of the image, and thus cannot be captured by the detector. Adjusting the optical setup to M=$\frac{1}{2}$, we saw a significant improvement: we were able to reveal the radiuses with the algorithm (asterisks in Fig. 4(b)), achieving a good fit between them and the theoretical MFP’ curve (solid line in Fig. 4(b)) even for low µ’_s values.

 

Fig. 4. The influence of magnification on the radius that discriminates between ROIs (theoretical and experimental results are the lines and asterisks, respectively). For (a) M=1 lower µ’_s deviant from theory and for (b) M=$\frac{1}{2}$ lower µ’_s fit theory. Inset: The optical setup of each configuration is presented above each graph, where the magnification M is changed by adjusting the distances before and after the lens.

Download Full Size | PDF

Next, by changing the distances between the camera, the lens and the phantoms, we created several magnifications (M=1, $\frac{1}{2},\;\frac{1}{3},\;$2 and 3, presented in Fig. 5(a)-(e), respectively). We present the theoretical phase RMS graphs of the ring and the center of the phase RMS (red and blue solid lines, respectively), together with the experimental results (red and blue asterisk, respectively). We see that the basic IMOPE technique, for M=1 (Fig. 5(a)), cannot yield ring’s values for low µ’_s. In addition, the high µ’_s values suffer from saturation. The limitations at both edges of the µ’_s axis narrow the linear ange of informative µ’_s in which the analysis is possible (between 0.7mm^-1 to 1.77mm^-1). These limitations were expected, due to the image size limitation; When $\mu _s^{\prime}$ is low, the image captured by the detector is smaller than the MFP’ (Fig. 4(a)). When the $\mu _s^{\prime}$ values are high, they are limited due to low number of pixels which fulfill the condition ρ<MFP’, causing the differences between the phase RMS values become negligible. For M<1 (Fig. 5(b) and (c)), we see that the experimental results, i.e., the phase RMS of the ring and the center (red and blue asterisks, respectively) are in high agreement with theoretical model. The µ’_s values in which the agreement occurs, are lower for M<1 than for M=1. However, for higher µ’_s we achieve saturation at lower values. For M=$\frac{1}{2}\;$ (Fig. 5(b)) the informative µ’_s range is between 0.5mm^-1 to 1.4mm^-1, and for M=$\frac{1}{3}\;$ (Fig. 5(c)) between 0.25 mm^-1 to 1.23mm^-1.

 

Fig. 5. The influence of magnification on the phase RMS. The above figures present the phase RMS versus µ’_s graphs for image magnification M of (a) $1$, (b) $\frac{1}{2}$, (c) $\frac{1}{3}$, (d) $2$ and (e) $3$. The solid lines present the theories graphs, blue for the center ROI and red for the ring ROI. The experimental results of the center and ring ROIs are represented by blue and red asterisks.

Download Full Size | PDF

Zooming-in i.e., configurating the setup to M>1 (Fig. 5(d) and (e)) solved the sensitivity problem expressed by the saturation at the higher µ’_s values as it shifts the linear range to the left. For M=2, the ring’s phase RMS is in good correlation with the theory, in which the µ’_s ranges between 1.03mm^-1 to 2.03mm^-1 and for M=3 it seems that the informative µ’_s range is above 2.03mm^-1. The phase RMS graphs of the centers in this range show a linear behavior with some offset from the theoretical values, as expected since the diffusion-based theory is inaccurate when getting closer to the light source [33]. In other words, this deviation does not indicate a mismatch of the method, but an inaccuracy of the DA based theory describing the area close to the center. Nonetheless, in the low µ’_s values the ring ROIs were not captured by the detector (red points are missing) and the center values are saturated. Table 1 summarizes the µ’_s ranges that were found to be most informative due to their linear behavior. As mentioned previously, the GS algorithm requires N images of the same phantoms, taken from N planes with a distance Δz between them. All the above experiments were done for different magnifications, based on the same constant distance Δz=0.635mm. We wanted to verify that the deterioration to the fit between the experimental results and the theoretical graphs for M>1 (Fig. 5(d)) and (e)) does not point that the Δz should have been magnified according to M. Hence, for $M = 1,\;\frac{1}{2},\;\frac{1}{3}$ magnifications (Fig. 6(a), (b) and (c) respectively) we compared the theoretical results with the experimental results of Δz (Fig. 6, blue and red asterisks for center and ring, respectively) and 2Δz (Fig. 6, blue and red circles for center and ring, respectively). Note that the above analysis presents a case in which the reconstruction depth (i.e. (N-1)Δz) of each magnification remains the same when changing the number of images. The results show that there is no significant difference between the phase RMS accepted for the different analysis methods (Δz vs. 2Δz).

 

Fig. 6. Influence of distance between image planes on phase RMS. Note that we scan the same total distance, hence the double distance analysis (blue and red circles for the double distance analysis of the center and the ring respectively) includes half the pictures then the regular analysis (blue and red asterisks for the single distance analysis of the center and the ring respectively). The figures present different magnifications: (a) $M = 1$ (b) $M = \frac{1}{2}$ and (c) $M = \frac{1}{3}\; .$ For all cases, the blue solid lines present the center ROI theory and the red solid lines present the ring ROI theory.

Download Full Size | PDF

Table 1. A summary of the $\mu _s^{\prime}$ ranges that were found to be most informative due to their linear behavior.

View Table

In this paper, we examined the influence of the optical magnification of the IMOPE’s optical system on its linear range. The change in magnification effectively changes the pixel size, and hence influences the RMS statistics of the phase. A theoretical DA based models for the intensity profiles, combined with Eq. (4)) for the phase theoretical model, predicted a shift in the linear range of the phase RMS values in both, the center (ρ<MFP’) and ring (ρ>MFP’) ROIs of the phase image. Experimental results from phantoms with varying scattering coefficients confirmed this behavior. However, some offset from the expected theoretical values were observed. The change of the optical magnification therefore extends the range of detection or enhances the detection accuracy for different µ’_s ranges and may also enable scanning deeper into the medium. To verify this, we had to first check that the µ’_s we have used meet the conditions required for the phase reconstruction. The change in the phase is extracted from the intensity image captured by a detector with a pixel size Δx. Each pixel contains data for the phase change along Δz. Thus, Δz must be large enough to enable a phase change along it, but small enough so there will not be too many phase changes averaged by the pixel. The condition on Δz is given by:

In the beginning, we reconstructed the phase from the phantoms using the same analysis for each magnification; The analysis used the same distance Δz=0.635mm between the intensity planes and the same number of planes (N=7) recorded by the detector scanning a total depth Δz⋅(N-1)=D=3.81mm for each phantom. Later, when we changed the number of images N and the distance Δz scanning the same depth (D=3.81mm), the phase RMS of each magnification remained. In other words, we have shown that the results accepted for different Δz and N but with constant reconstruction depth D yield the same phase RMS experimental results (Fig. 6). In addition, we would like to examine how by changing the reconstruction depth, taking different Δz with the same number of images, result with different reconstruction depths for the same magnification, meaning, different scanned depth D. Therefore, for M=$\frac{1}{3}$ we compared two reconstruction depths based on the same number of images with different distances between them: single and triple distance Δz between the image planes (asterisks and circles respectively in Fig. 7). The distances between the images are therefore 0.635 and 1.905mm for single, and triple distances. Hence the total depths D are: 3.81 and 11.43mm. The reconstruction along the different depths had yield the same phase RMS values for the same phantoms, confirming the Δz we chose fulfills the condition of Eq. (5)). Note, that the results shown in Fig. 7 were achieved from a M=$\frac{1}{3}$ configuration. The center is therefore smaller, and the average of the pixels is too coarse. This explains the deterioration of the fit between the theory and the experimental results. In contrary, the ring theory and results highly agree along the µ’_s axis.

 

Fig. 7. Influence of reconstruction depth on the phase RMS. The reconstruction depth was changed by taking the same number of images and changing the distance between them; single and triple distances Δz between the images (corresponding to 0.635 and 1.905 mm represented by asterisks and circles, respectively). The theoretical graphs (lines) and experimental results (points) of image magnification M=$\frac{1}{3}$ are presented. The red and the blue lines and points represent the ring and the center ROIs, respectively.

Download Full Size | PDF

5. Conclusion

The IMOPE technique aims to extract the reduced scattering coefficient µ’_s from a phase RMS reconstructed from a turbid medium. It is based on a multi-planar version of the Gerchberg Saxton algorithm for phase retrieval. The optical setup takes several images that are inserted to the algorithm that results in a reconstructed phase image. The phase image is analyzed and from a comparison with the theoretical model, µ’_s can be extracted. In this research, we have searched for the IMOPE’s optimal range of detection, for different magnifications. We show that in different ranges of µ’_s the technique’s linear range is limited: at the lower values of µ’_s it is limited by the size of the image captured by the detector, where at the higher µ’_s values it is limited due to the saturation the phase RMS reaches at low number of pixels. The work presented here aims to improve the capabilities of the IMOPE technique for the detection of lower and higher reduced scattering coefficients by changing the magnification of the image. In fact, we suggest adaptation of the setup in order to tailor the technique for desired µ’_s range (summarized in Table 1). We did so by varying the magnification of the image and finding the linear µ’_s range for each magnification. The phantoms experiments have supported this theoretical work as well. For a magnification of M<1, the sensitivity for lower µ’_s is improved since they are included at the linear range. For M>1 magnification, the linear range includes higher µ’_s values. Since these values do not suffer from saturation anymore, the sensitivity of the technique for these ranges has improved. In addition, we notated the cross-point between the center and ring phase RMS curves that shift with the magnification M. Overall, this paper shows that the IMOPE, designed to detect changes in biological tissues, can be extended to other fields beyond its original purpose; It can be tailored to detect wider ranges of µ’_s than necessary for biological research as well as larger distances along the turbid media. In the blue wavelength, for example, it could be utilized for underwater research. The IMOPE itself is a promising sensing method, that breakthrough the imaging limitations and achieves information about invisible area. Here we show the great potential of the IMOPE to be extended to various applications.