1. Introduction

Probing optically inside opaque materials seems at first glance to be a contradiction in terms. Opacity is defined in layman’s terms as the capacity to block the passage of light [1], in which case probing inside an opaque material seems impossible. However opaque materials such as paint or human tissue do not block the passage of all of the impinging light but instead scramble the direction of light rays via multiple scattering. While this effect blocks simple direct imaging, it still allows randomized diffuse light to penetrate deep inside. This deeply penetrating light can be used to as a probe.

There currently exist a number of methods for using diffuse light to probe inside multiple scattering media. For example, diffuse imaging is a well developed technique for mapping out the structure of opaque human tissue [2, 3]. This technique typically uses fibers to create a point source of light at the the skin surface and collect the outgoing light from points on the skin of varying distance from the source. The resulting spectral and spatial information provides information about the structure of the tissue, including, for example the spatial distribution of certain biologically significant molecules such as oxy-haemoglobin (HbO) and deoxy-haemoglobin (HHb)) [4]. With enough sources and detection points, a full tomographical image can be inferred albeit with limited resolution.

In diffuse imaging, one can ask two questions about the light exiting at some distance from a focused source. First, how much light is exiting and second how long did it take for the light to exit? Answers to these two questions can be used to separate absorption from scattering effects. For example, consider two slabs, one with high absorption and one with strong scattering. Far from the source, both will show reduced amounts of light exiting the material, but for different reasons. Strong scattering implies a short mean free path and therefore, for a given detection distance, long times of flight. The reduction in intensity occurs because short paths are reflected out of the surface. In contrast, strong absorption implies shorter average times of flight, since long paths are converted to heat. Therefore both the absolute intensity and time sensitive information are useful for characterizing what is happening to the the light inside the material.

Different methods have been used for obtaining time resolved information depending on the type of material. For biological samples such as human tissue a pulsed [5] or radio frequency modulated [6] source allows for time resolved detection. In strongly scattering materials such as paint, where typical photon times of flight may be beyond the temporal resolution of electronic detectors, interferrometric methods have been used to obtain time resolved information [7]. These methods have been applied to single spatial speckle spots to give information about diffusion, absorption, as well as localization effects. Information about the time of flight is also contained in the frequency speckle for a broad band source [8]. Frequency speckle gives information about decorrelation of phases - the degree to which two photons of different frequencies decorrelate depends on the time of flight relative to the magnitude of the frequency difference. The longer the path, the smaller the frequency difference needed for decorrelation. A simple way to characterize this decorrelation is the intensity-intensity correlation function. Recently Muskens et. al. have shown that determination of this correlation function can be used to determine the optical diffusion constant over a broad band of frequencies [9]. Such measurements are direct and robust compared to interferrometric methods.

In this paper we demonstrate a simple diffuse imaging setup for strongly scattering materials. The approach gives both the average intensity and the frequency correlations for a point source of white light as a function of position across the sample. We use a low cost equipment that can be attached as a module to a standard optical microscope. As a test of this method we apply it to the low-absorbing thick white paint. We show that for such a sample, we can model the results based on an simple analytical theory for diffuse light. The technique should be potentially useful for understanding strongly scattering structured materials.

2. Experiment

2.1. Measurement

Samples each with a uniform layer of white paint layer were prepared via dip coating. A 150 micron thick microscope slide was cleaned with acetone. The paint, (AKZO Nobel Flexa 2025 wit zijdeglans) was mechanically stirred for five minutes in a small beaker. Air bubbles were removed from the surface by pricking. The microscope slide was then submerged vertically in the paint while avoiding any remaining bubbles on the surface. The slide was removed and held vertically for at least 10 seconds to allow excess paint to drip off. The sample was then left to dry horizontally in a dust free environment for twenty-four hours.

The thickness and surface roughness of the paint layer were measured. First a long (>1 cm) thin (<0.5 mm) scratch was made in the surface to expose the glass underneath. Then profilometer scans were made at several locations along the scratch. The results show that the dip coating technique produced uniform smooth layers with surface variation on the order of 1 micron. The overall thickness of the paint could be controlled via the dripping time with resulting layers from 50 to 100 microns for dripping times ranging from 10 to 60 seconds respectively.

The spectrally resolved image was created using a collimated wavelength-tunable light source and a microscope (Fig. 1). The light source is a supercontinuum laser (Fianium SC400-4), which has enough spatial coherence to allow for tight focusing on the surface of the sample. The beam passes though a Glan-Taylor linear polarizer and is then sent into a high resolution collinear beam acousto-optic tunable filter (AOTF) (Crystal Technology CBAOTF 97-02768-01) creating a collimated tunable monochromatic beam. The wavelength resolution of the AOTF is less than 0.6 nm over the entire wavelength range. This monochromatic beam is spatially filtered using a single mode polarizing preserving optical fiber. A small percentage of the resulting beam is sent to a reference photodiode via a beam splitter. The main portion of the beam is steered into a microscope where it is focused to a 3 micron spot on the surface of the paint sample. This focused light diffuses into the paint and radially away from the focus spot. The light exiting the top surface of the sample passes through a second polarizer (the analyzer) and is then imaged onto a CCD camera (Allied Vision Technologies Pike F-145). A PC controls both the CCD and the AOTF to collect images for each wavelength step.

 

Figure 1 Diagram of the experimental setup. The tunable source (inside the dotted lines) is created by combining a supercontinuum laser with a high resolution acousto-optic tunable filter (AOTF). The light from the laser passes through a polarizer (P), followed by the AOTF. The single wavelength from the AOTF is spatially filtered (SP). A small percentage is then split off to a reference photodiode (R) with the bulk of the light traveling to a polarizing microscope, where it is focussed to a spot. The diffuse light is imaged onto a CCD for each wavelength.

Download Full Size | PDF

The apparatus can scan quickly through the entire wavelength range resulting in a hyperspectral image of the diffuse light. The high resolution of the AOTF allows for frequency speckle to be resolved. The limiting speed factor of the measurements is the signal to noise ratio at the CCD and the wavelength resolution of the scan. In principle full scans of several seconds are possible. For our measurements, in which speed was not an issue, integration times and image averaging lead to full scans of several minutes. Each CCD image was taken with a 1 ms exposure time. To reduce the signal to noise, ten images were averaged together for each wavelength. The AOTF was stepped through the wavelength range in steps of 0.2 nm, intentionally below the resolution of the AOTF both to increase signal to noise and to allow us to track the effect of the AOTF resolution on the results.

Before measuring the paint characteristics two reference measurements were taken. First a background (“dark”) signal was measured by blocking the laser beam. The resulting image consists mainly of dark noise from the camera with a small amount of photo current presumably from scattered light from the setup. The camera also showed a small degree of striping. These inhomogeneities could be effectively subtracted from the paint data. Following the dark measurement a measurement was taken using a silver mirror at the sample location (“reference”). This gives the spectrum of the laser as modified by the setup and is used to normalize the paint data. The mirror measurement also allowed us to record the quality of the focus spot for all wavelengths measured. A clean, roughly circular Gaussian spot was observed for the entire wavelength range with a full width at half maximum of approximately one micron. Finally, the mirror measurement allows for an estimate of the light scattered from the setup back to the surface (as opposed to the light that diffuses through the sample). This will be discussed further in the results section.

Following the reference measurement, measurements of the paint sample were made. The measurements were taken for a) both parallel and perpendicular conditions of the analyzer b) for two different incoming intensity levels, a low intensity level to be used for small radii (ρ < 10 microns) and a high intensity level used for large radii (ρ > 10 microns). The high intensity level was roughly ten times the low level and saturated the camera for small radii. For each condition, a measurement was taken at each of three different locations of the sample to test for consistency and to estimate the uncertainty. At the end of the paint measurements a second reference measurement was made. This gave a check to ensure that no systematic drift of the laser occurred during the measurement.

2.2. Data analysis

The large quantity of data is analyzed in terms of average values and correlation of fluctuations (speckle). A frequency scan results in a matrix of intensity values I (ρ_i, ω_j) where ρ_i = (x, y)_i is the vector position of the i^th pixel measured from the center of the laser focus and ω_j is the j^th discrete optical frequency where ω_j = 2πc/λ_j. In our notation we will use bold font for vectors and normal font for the scalar magnitude. I (ρ_i, ω_j) can be obtained for both analyzer polarizations, parallel ‖ and perpendicular ⊥ to the polarization of the incoming light.

For statistical analysis, we analyze subsets of this data over given spatial and frequency ranges using Matlab routines. Each spatial range is a ring of pixels about the focal point with average radius ρ̄ and width Δρ such that $\left|{\rho }_i-\overline{\rho }\right|\le \frac{\mathrm{\Delta }\rho }{2}$ where ρ_i = |ρ_i| is radius of the i^th pixel. Each frequency range has a width Δω about a given average frequency ω̄ such that |ω_j − ω̄| ≤ Δω/2. The ranges are chosen to be large enough to provide good statistics, but narrow enough so that the optical properties of the material do not vary substantially within the range. In this way fluctuations due to speckle are assumed to arise from independent mesoscopic realizations of the same macroscopic material parameters allowing for standard statistical methods to be applied.

Two average intensity functions are useful in the analysis, the pixel intensity at a given position ρ_i averaged over a given frequency band ω̄ ±Δω/2,

The frequency autocorrelation function at a given pixel position ρ_i is calculated by

3. Theory

3.1. Analytic expressions

The equation of radiative transfer describes the temporal and spatial evolution of the radiance I_s (r, t, ŝ), the intensity per unit solid angle for a given position r at time t and direction ŝ [10]. The equation includes the propagation and scattering properties of electromagnetic waves, but does not include interference effects such as speckle, enhanced backscattering, or localization. It can however be used to deduce speckle statistics since it gives informaton about average light path lengths and times of flight.

The equation of radiative transfer can under certain conditions be simplified to a diffusion equation for the average intensity $U(\mathbf{r},t)\equiv \frac{1}{4\pi }{\int }_{4\pi }d\mathrm{\Omega }{I}_s(\mathbf{r},t,\widehat{\mathbf{s}})$ given as

Our experimental geometry can be approximated by an infinitely thick slab of thickness onto which light has been focused to a point at position ρ = 0 where ρ is the 2d position vector on the surface of the sample, and a depth of z = z_p. We therefore assume a delta function source of the form S (r, t) = δ(ρ)δ(z_p) Θ(t), where r = (ρ_x, ρ_y, z) and Θ(t) is the step function. We rescale all length units by the mean free path (ρ ⇒ ρ/l, z ⇒ z/l) and time units by $\frac{3}{4}{\tau }_{mf}\left(t\Rightarrow t/\left(\frac{3}{4}\right){\tau }_{mf}\right)$, ${\tau }_a\Rightarrow {\tau }_a/\left(\frac{3}{4}{\tau }_{mf}\right))$ and use the relation D = l²/ (3τ_mf) to simplify the expressions. The method of images yields an analytic solution for the average intensity propagator from the source to position (ρ, z):

The reflected diffuse intensity exiting the sample is given by Fick’s law $R=-D{l}^{-1}\frac{\partial }{\partial z}{U|}_{z^{\prime }=0}$ [13]. Applying this to Eq. (8) yields:

Up to this point, we have not included polarization effects. The polarization can be easily analytically incorporated into the above equations as follows. The initial polarized state is assumed to decay to unpolarized light exponentially with a time constant ofτ_p. This gives:

This substitution can be carried through the entire calculation, yielding

Note also that the entire calculation can be carried out for a finite slab of thickness L_slab by adding an additional infinite series of image charges. This is achieved with no further calculation via the substitution dz_n ⇒ dz_n,_m where dz_n,_m ≡ (1 − 2n)z_p − z −2nz₀ + 2m(L + 2z₀) and adding the infinite summation over m such that

3.2. Simulated speckle

Simulated speckle is also useful for comparing diffusion theory with the data. The simulated speckle allows for incorporating into the model two features of the real measurement that are not incorporated in the analytical derivations. The first is the finite resolution of the spectral scan. The second is the limited frequency window Δω over which the correlations can be calculated for a given value of ω̄. This second limit affects the data for small values of ρ where the frequency correlations are very wide. Simply widening Δω indefinitely is not an option, since the real sample scattering and absorption parameters are frequency dependent.

The reflectivity as a function of time can be calculated from Eq. (7) using Fick’s law. This yields

The simulated frequency speckle can be generated using Eq. (14) as follows [16, 17]. A sequence of Gaussian distributed random complex numbers G_i is generated and each random number associated with a point in a time series t_i. The simulated temporal speckle amplitude E_sim (ρ, t) is then generated for each value of ρ and t_i via $E_{\mathit{sim}}(\rho ,t)=G_i\sqrt{R{(\rho ,t_i)}_{{\tau }_a}}$. The simulated frequency speckle amplitude E_sim (ρ, ω) is then computed via numerical (fast) Fourier transform of E_sim (ρ, t). The intensity speckle is the absolute square of this I_sim (ρ, ω) = |E_sim (ρ, ω)|². This intensity speckle may be used to simulate correlation functions following the routines described in section 3.1.

4. Results

Figure 2 gives the measured intensity function I_⊥ (ρ_x, ρ_y, ω) for ω = 2.99 fs⁻¹, (λ = 630 nm). The raw data (inset) shows the spatial intensity speckle with the expected exponential intensity distribution statistics. These speckles have been digitally filtered out via a low pass digital square boxcar filter of size 2.5 × 2.5 μm² to reduce the spatial speckle. The isointensity lines are logarithmically distributed, demonstrating a large dynamic range for the measurement. They form nearly perfect circles about the focus point, indicating radially symmetric propagation of light in the plane of the sample surface and good alignment of the sample surface with the image plane of the microscope. The roughness of the isointensity lines is the result of incomplete filtering of the speckle and not, for the most part, of sample roughness or dark noise of the camera.

 

Figure 2 Reflectivity per unit μm² the sample (smoothed in contour plot and not smoothed in the inset). The isointensity lines starting from the center are 64, 32, 16, 8, 4, 2, 1, 0.5, 0.25 × 10⁻³. The main data has been filtered with a low pass digital square boxcar filter of size 2.5 × 2.5 μm² to reduce the spatial speckle. The inset shows a 2.5μm square image of the unfiltered data near the focus, which shows the speckle pattern.

Download Full Size | PDF

The isotropic nature of the data allows for azimuthal averaging. Figure 3 contains the radially averaged data Ī(ρ̄, ω) with Δρ = 1 μm and ω = 2.99 × 10³ THz for each of the two polarizations. The two polarizations show different behavior near the center and then merge completely beyond ρ ≈ 10 μm. This is expected, since the light immediately exiting the sample at ρ = 0 contains many paths that have undergone just one or two scattering events and are thus not yet depolarized. Away from the focus the all light paths contain many scattering events and the light is completely depolarized. No dramatic specular reflection effects were visible by eye for these samples. By definition, such effects would occur only at the focal spot, i.e. within a radius of 0.5 microns. At the smallest averaged radius, 1 micron, the data is slightly higher than the fit for parallel polarized light and lower for perpendicular. However this slight difference falls within the error bars and thus does not affect the fitting results. Enhanced backscattering in principle also slightly biases the total refection towards the parallel polarization configuration. However the effect is so small for mean free paths on the order of microns and angular apertures as large as ours that it should not alter the fitting results.

 

Figure 3 The diffuse intensity, averaged azimuthally, as a function of radius for parallel (circles) and perpendicular (triangles) orientation of the analyzer with the fit (line) to Eq. (12) with ω = 0. The errors bars are the standard error calculated from the azimuthal averaging. The inset shows the same plot at larger radii to highlight the quality of the fit even above 30 microns.

Download Full Size | PDF

The data was fit to Eq. (12) with ω = 0 as shown in Fig. 3. The value of z₀ was fixed at 2.42 by assuming an effective refractive index of the paint of n_eff = 1.5, while the parameters A, τ_a, τ_p, and l were allowed to float freely. The quality of the fits for ρ ≤ 30 μm were excellent as seen in the example in Fig. 3 with χ² values on the order of 1. This demonstrates the effectiveness of the diffusion model for characterizing the data. Extending the fitting range beyond ρ = 30 μm lead to increases in χ² as the fit could visibly not account for the data for the entire radius range.

The intensity plotted in Fig. 3 drops non-exponentially until roughly 15 microns, after which the behaviour is exponential (i.e. linear on the semi-logarithmic plot). The exponential behaviour at large radii indicates the point at which a large fraction of the emitted light is absorbed rather than scattered back out of the sample. Indeed the slope of the tail of Fig. 3 can be used to directly calculate the absorption depth L_a. In this article we use the fit to the full theoretical expression to calculate material parameters.

The resulting material parameters are plotted in Fig 4. The absorption time τ_a (in units of mean free time), decreases systematically while the depolarization time τ_p and the mean free path l remain constant. The bottom plot shows the goodness of fit for the entire range.

 

Figure 4 Resulting fitting parameters over the green, yellow, and red wavelengths. The absorption time τ_a decreases systematically while the depolarization time τ_p and the mean free path l remain constant. (both τ_a and τ_p are given in units of 0.75τ_mf), The bottom plot shows the goodness of fit for the entire range.

Download Full Size | PDF

Figure 5 shows the raw data vs. optical frequency at 3 spatial locations, with increasing radii from top to bottom (ρ =0 μm, 4 μm, and 10 μm for the top, middle, and bottom plots respectively). The near zero crossings indicate sufficient measurement resolution for applying single speckle statistics. It is already apparent that smaller radii show correlations over a wider range. This corresponds to our theoretical expectations. The smaller radii contain more short paths that show a high degree frequency correlation. The top plot also shows one difficulty of analyzing the data. The correlations begin to span the entire range of the scan, while we know from the fit to average intensity that the parameters (the absorption time) may change considerably over this span.

 

Figure 5 Selected speckle patterns at 3 different radii: ρ = 0, 4, and 10 microns for the top, middle, and bottom plots respectively. The frequency correlation decreases markedly as the radii increases.

Download Full Size | PDF

The correlations functions C̄_ρ̄,ω̄, (δω_k) for several values of ρ̄ and ω̄ are plotted in Fig. 6. In the top plot, these functions drop more rapidly to zero as the radius decreases signifying a reduction in correlation and therefore longer average time of flights. This is expected due to the greater distance from the focus point. By contrast, C̄_ρ̄,ω̄, (δω_k) is fairly insensitive to variations in ω̄ at constant r̄ as seen in the lower half of Fig. 6.

 

Figure 6 Top: C̄_ρ̄,ω̄, (δω) for several values of ρ̄ (0, 4, 8, and 14 microns for the squares, circles, up triangles, and down triangles respectively) at a constant average frequency of ω̄ = 3.03 fs⁻¹(λ̄ = 0.620 μm) averaged over a range of δω = 0.107 fs⁻¹(δλ = 0.022 μm). The correlation becomes weaker (narrower) as the radius increases. Bottom: C̄_ρ̄,ω̄, (δω) for several values of ω̄ (3.03, 3.14, 3.25, and 3.35 fs⁻¹) for the squares, circles, up triangles, and down triangles respectively) at a constant radius of ρ̄ = 6 μm. The frequency dependence of the correlation function at this radius is below the measurement uncertainty.

Download Full Size | PDF

In order to compare theory to data, we plot the half maximum of the correlation functions as a function of radius in Fig. 7. We have considered only the perpendicular case here. The trend towards narrower correlation functions with increasing radii is clearly demonstrated for both measurement and theory. There is also a slight dependence in the data on the average frequency. The analytical result from diffusion theory, Eq. (12), yields the dashed line in Fig. 7. The same parameters as the intensity data fit along with the value τ_mf = 16 fs were used with no free parameters. The mean free time was calculated from τ_mf = l/(c/n_eff) with the same value of n_eff as used previously (n_eff = 1.5). This calculation assumes no resonances in the material, i.e. that the energy and phase velocities are equivilant. The analytical result shows qualitatively the same trend as the data. Nevertheless, the comparison at both low and high values of the radius are off by as much as a factor of two.

 

Figure 7 The full width at half max (FWHM) of C̄_⊥ (r̄, ω̄, δω_k) as a function of radius crossed polarization case. The FWHM decreases with increasing radius due to the longer flight times in the material. The two theoretical curves are the analytical result (dashed line) from Eq. (12) and the simulated result (solid line) using the method described in section 3.2.

Download Full Size | PDF

The discrepancy can be explained as follows. At large radius, and therefore narrow correlation functions we may be approaching the frequency resolution of the setup. Once the resolution limit of the frequency scan is approached, the correlation function will flatten out faster than predicted by diffusion theory since frequencies separated by less than the resolution of the scan will be not well differentiated. The discrepancy at the low radius can be explained by the undersampling due to the finite size of the frequency window used in calculating the correlation function. This window was limited in size by the range over which the optical parameters of the material were expect to change (δλ ∼ 20–50 nm).

These conjectures about the discrepancy between the analytical expression and the data can be tested by simulating frequency speckle using the method described in Section 3.2. The results of this simulation are given by the dashed line in Fig. 7 The parameters were chosen to be reasonably close to those used for fitting the average intensity data (l = 2.8 μm, τ_a = 50, τ_p = 0.6, and, τ_mf = 20 fs) with a wavelength resolution of 0.9 nm to fit the data by eye. (Constructing a fitting routine using the simulation proved to be too time consuming). The results show that indeed the discrepancies with the analytical expression can be explained in part by resolution and processing limitations. However the data still does not overlap with the theory within the given error bars. This could point to a limitation of the theory itself. For example for short light paths, which dominate the signal at small radii, the diffusion theory itself may not be applicable and a more complete radiative transfer model may be needed. This will be an important consideration when analyzing more complex structures.

Dark current noise also influences the frequency correlation measurements. The signal from each pixel is a combination of photo-current and dark current. The dark current has zero correlation. Therefore, as the dark current overwhelms the photo current, the correlation width decreases to zero. This occurs at radii larger than 25 μm.

5. Conclusions

Frequency resolved diffuse imaging of a strong scattering material has been shown to yield results consistent with diffusion theory for white paint. The frequency resolved data can be analyzed using both the average intensity data as well as the frequency correlations. These give different information about the sample and allow the mean free path, mean free time, depolarization time, and absorption time to be probed. The current setup has been shown to be effective for efficient data collection in a modular unit and relatively inexpensive components.

Limitations to the current setup include the resolution of the frequency scan, the signal to noise of the intensity measurement (which limits the maximum radius for both intensity and frequency correlation measurements), and the limited band width for frequency correlations. The resolution of the frequency scan can be improved fairly easily by using for example a high resolution monochromator or a tunable laser. Both of these approaches could increase the frequency resolution to 0.1 nm or better, i.e. a factor of 10 improvement. The signal to noise could also be easily improved by a factor of ten by using a cooled camera and longer averaging times. The bandwidth limitation of the frequency scan is slightly trickier to overcome. We expect that averaging over multiple realizations of the sample, using for example an automated sample positioner would improve that problem.

We hope that this general approach can be applied to more complicated paint sample geometries. Colored paints with strong absorption bands should show unique radial intensity distributions and frequency correlation functions that allow for detailed characterization. Multilayered structures should be accessible to depths on the order of 20 microns, yielding depth dependent spectra. Finally, the possibility of true diffuse imaging of complex samples, such as multilayered master oil paintings [18] seems possible.

Furthermore, the technique can be applied to any multiple scattering material for which the optical mean free path is on the order of microns to tens of microns. Examples include powders, pills, ceramics, bone, and paper products to name a few.