## Abstract

A new probe of multiple scattering material is demonstrated experimentally. Light from a tunable wavelength source is focused to a point on the surface of an opaque slab. A fraction of this light penetrates into the slab, is multiply scattered, and reemerges at the surface creating a surface speckle pattern. The full spatial and frequency speckle can be easily and quickly recorded using a CCD and an acoustooptical tunable filter. Both the average intensity and frequency correlations of intensity are analyzed as a function of the distance to the source. This method is demonstrated experimentally for white paint. The resulting model yields information about both the static and dynamic transport properties of the sample. The technique has prospects for both static and time resolved diffuse imaging in strongly scattering materials. The setup can be easily used as an add-on to a standard bright field microscope.

© 2014 Optical Society of America

## 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.

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}*,

*ω*) where

_{j}*= (*

**ρ**_{i}*x*,

*y*)

*is the vector position of the*

_{i}*i*pixel measured from the center of the laser focus and

^{th}*ω*is the

_{j}*j*discrete optical frequency where

^{th}*ω*= 2

_{j}*πc/λ*. In our notation we will use bold font for vectors and normal font for the scalar magnitude.

_{j}*I*(

*,*

**ρ**_{i}*ω*) can be obtained for both analyzer polarizations, parallel ‖ and perpendicular ⊥ to the polarization of the incoming light.

_{j}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}* = |

*| is radius of the*

**ρ**_{i}*i*pixel. Each frequency range has a width Δ

^{th}*ω*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,

*r̄*±Δ

*r*/2 and frequency band

*ω̄*±Δ

*ω*/2,

*Ī*(

*ρ̄*,

*ω̄*) is its standard deviation divided by the square root of the number of speckle spots in the ring defined by

*ρ̄*± Δ

*ρ*/2. The area of the speckle spot in pixel units is determined from the full width half maximum of the spatial autocorrelation each image.

The frequency autocorrelation function at a given pixel position * ρ_{i}* is calculated by

*δω*is given by the discrete frequency difference

_{k}*ω*

_{j+k}−

*ω*. This function is then averaged over all pixels in the range Δ

_{j}*ρ*to give the average autocorrelation function,

*δω*is frequency independent. By definition,

_{k}*C̄*

_{ρ̄,ω̄}, (0) = 1. When presenting the experimental data, some of the above expressions will be subscripted with ‖ and ⊥ for the parallel and perpendicular analyzer orientations respectively. Nonlinear least squares fitting routines from Matlab were used to fit the processed data to theoretical predictions.

## 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

*D*is the diffusion constant,

*τ*is the absorption time, and

_{a}*S*(

**r**,

*t*) is the source function [11]. This simplification is valid when

*U*(

**r**,

*t*) is much larger than the diffuse flux of energy in any direction and the absorption time is small compared to the mean free time between scattering events. The diffusion equation has been applied effectively to a wide variety of multiple scattering problems providing that the surface behavior is treated with care. For short paths near the surface and for extreme absorption this approximation breaks down, but for the low absorbing samples and long paths, as will be important for this article, the diffusion approximation is valid.

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*. We therefore assume a delta function source of the form

_{p}*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*

^{2}/ (3

*τ*) to simplify the expressions. The method of images yields an analytic solution for the average intensity propagator from the source to position (

_{mf}**,**

*ρ**z*):

*dz*≡ (1 − 2

_{n}*n*)

*z*−

_{p}*z*− 2

*nz*

_{0}, and

*z*

_{0}is the extrapolation length in units of

*l*[12]. Equation (7) is a scalar quantity. The frequency dependent propagator can be found by Fourier transforming Eq. (7) yielding:

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:

*Ī*(

*ρ̄*,

*ω̄*), the average (static) intensity from Eq. (2) can be directly compared with

*R*(

*ρ*,

*ω*= 0). The correlation function

*C̄*(

*ρ̄*,

*ω̄*,

*δω*) from Eq. (5) can be directly compared to |

_{k}*R*(

*ρ*,

*ω*= 0)|

^{2}.

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

*C*

_{‖,⊥}(

*ρ*,

*ω*) = |

*R*

_{‖,⊥}(

*ρ*,

*ω*)

_{τa,τp}|

^{2}. This simple analytic dependence on

*τ*allows analytic calculation of any polarization dependent quantity and is therefore somewhat simpler than previous descriptions of the polarization dependence of diffuse light [14, 15].

_{ap}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}*where*

_{m}*dz*,

_{n}*≡ (1 − 2*

_{m}*n*)

*z*−

_{p}*z*−2

*nz*

_{0}+ 2

*m*(

*L*+ 2

*z*

_{0}) 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*. The simulated temporal speckle amplitude

_{i}*E*(

_{sim}*ρ*,

*t*) is then generated for each value of

*ρ*and

*t*via ${E}_{\mathit{sim}}(\rho ,t)={G}_{i}\sqrt{R{(\rho ,{t}_{i})}_{{\tau}_{a}}}$. The simulated frequency speckle amplitude

_{i}*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}*ρ*,

*ω*)|

^{2}. 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

^{−1}, (

*λ*= 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

^{2}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.

The isotropic nature of the data allows for azimuthal averaging. Figure 3 contains the radially averaged data *Ī*(*ρ̄*, *ω*) with Δ*ρ* = 1 *μ*m and *ω* = 2.99 × 10^{3} 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.

The data was fit to Eq. (12) with *ω* = 0 as shown in Fig. 3. The value of *z*_{0} was fixed at 2.42 by assuming an effective refractive index of the paint of *n _{eff}* = 1.5, while the parameters

*A*,

*τ*,

_{a}*τ*, and

_{p}*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

*χ*

^{2}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

*χ*

^{2}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

*τ*and the mean free path

_{p}*l*remain constant. The bottom plot shows the goodness of fit for the entire range.

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.

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̄*

_{ρ̄,ω̄}, (

*δω*) is fairly insensitive to variations in

_{k}*ω̄*at constant

*r̄*as seen in the lower half of Fig. 6.

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*) with the same value of

_{eff}*n*as used previously (

_{eff}*n*= 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.

_{eff}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,

*τ*= 0.6, and,

_{p}*τ*= 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.

_{mf}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.

## Acknowledgments

This work was supported by the Dutch Technology Foundation STW and by FOM. We also received support from Akzo Nobel, Fianium Ltd., and the Netherlands Institute for Cultural Heritage.

## References

**1. **Merriam Webster online, http://www.merriam-webster.com/.

**2. **T. Durduran, R. Choe, W. B Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Reports on Progress in Physics **73**(7), 076701 (2010). [CrossRef]

**3. **A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. **50**, R1–R43 (2005). [CrossRef] [PubMed]

**4. **F. F. Jobsis, ”Noninvasive, infrared monitoring of cerebral and myocardial oxygen, sufficiency and circulatory parameters,” Science **198**(4323), 1264–1267 (1977). [CrossRef] [PubMed]

**5. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt **28**(12), 2331–2336 (1989). [CrossRef] [PubMed]

**6. **M. S. Patterson, J. D. Moulton, Br. C. Wilson, K. W. Berndt, and J. R. Lakowicz, “Frequency-domain reflectance for the determination of the scattering and absorption properties of tissue,” Appl. Opt. **30**(31), 4474–4476 (1991). [CrossRef] [PubMed]

**7. **P. M. Johnson, A. Imhof, B. P. J. Bret, J. Gómez Rivas, and A. Lagendijk, “Time-resolved pulse propagation in a strongly scattering material,” Phys. Rev. E **68**(1), 016604 (2003). [CrossRef]

**8. **A. Z. Genack and J. M. Drake, “Relationship between optical intensity, fluctuations and pulse propagation in random media,” EPL (Europhysics Letters) **11**(4), 331–336 (1990). [CrossRef]

**9. **O. L. Muskens and A. Lagendijk, “Method for broadband spectroscopy of light transport through opaque scattering media,” Opt. Lett. **34**, 395–397 (2009). [CrossRef] [PubMed]

**10. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic Press, 1978).

**11. **M. C. W. van Rossum and Th. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion,” Rev. Mod. Phys. **71**(1), 313–371 (1999). [CrossRef]

**12. **J. F. de Boer, Optical fluctuations on the tranmsission and reflection of mesoscopic systems (PhD thesis, University of Amsterdam, 1995).

**13. **D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. i. theory,” Appl. Opt. **36**(19), 4587–4599 (1997). [CrossRef] [PubMed]

**14. **L. Fernando Rojas-Ochoa, D. Lacoste, R. Lenke, P. Schurtenberger, and F. Scheffold, “Depolarization of backscattered linearly polarized light,” J. Opt. Soc. Am. A **21**(9), 1799–1804 (2004). [CrossRef]

**15. **E. Akkermans, P. E. Wolf, R. Maynard, and G. Maret, “Theoretical study of the coherent backscattering of light by disordered media,” J. Phys. France **49**(1), 77–98 (1988). [CrossRef]

**16. **J. Aulbach, Spatiotemporal Control of Light in Turbid Media (PhD thesis, University of Twente, 2013).

**17. **D.J. Young and N.C. Beaulieu, “The generation of correlated rayleigh random variates by inverse discrete fourier transform,” Communications, IEEE Transactions on **48**(7), 1114–1127 (2000). [CrossRef]

**18. **M. Alfeld and J. A.C. Broekaert, “Mobile depth profiling and sub-surface imaging techniques for historical paintings - review,” Spectrochimica Acta Part B: Atomic Spectroscopy **88**(0), 211–230 (2013). [CrossRef]