Measuring isotropic subsurface light transport

Subsurface light transport can affect the visual appearance of materials significantly. Measuring and modeling this phenomenon is crucial for accurately reproducing colors in printing or for rendering translucent objects on displays. In this paper, we propose an apparatus to measure subsurface light transport employing a reference material to cancel out adverse signals that may bias the results. In contrast to other approaches, the setup enables improved focusing on rough surfaces (e.g. uncoated paper). We derive a measurement equation that may be used to deduce the point spread function (PSF) of subsurface light transport. Main contributions are the usage of spectrally-narrowband exchangeable LEDs allowing spectrally-resolved measurements and an approach based on quadratic programming for reconstructing PSFs in the case of isotropic light transport. © 2014 Optical Society of America OCIS codes: (110.4850) Optical transfer functions; (290.5820) Scattering measurements; (290.7050) Turbid media. References and links 1. J. A. C. Yule and W. J. Nielsen, “The penetration of light into paper and its effect on halftone reproduction,” Tech. Assn. Graphic Arts 4, 65–76 (1951). 2. G. L. Rogers, “Optical dot gain in a halftone print,” J. Imaging Sci. Technol. 41(6), 643–656 (1997). 3. G. Rogers, “Optical dot gain: lateral scattering probabilities,” J. Imaging Sci. Technol. 42(4), 341–345 (1998). 4. J. S. Arney, T. Wu, and C. Blehm, “Modeling the Yule-Nielsen effect on color halftones,” J. Imaging Sci. Technol. 42(4), 335–340 (1998). 5. J. A. S. Viggiano, “New models for the reflectance spectra produced by halftone-based hardcopy,” Ph.D. thesis, Center for Imaging Science, Rochester Institute of Technology, Rochester, NY, USA (2010). 6. M. Hašan, M. Fuchs, W. Matusik, H. Pfister, and S. Rusinkiewicz, “Physical reproduction of materials with specified subsurface scattering,” ACM Trans. Graphics 29, 61 (2010). 7. Y. Dong, J. Wang, F. Pellacini, X. Tong, and B. Guo, “Fabricating spatially-varying subsurface scattering,” ACM Trans. Graphics29, 62 (2010). 8. M. Goesele, H. Lensch, J. Lang, C. Fuchs, and H. Seidel, “Disco: acquisition of translucent objects,” in ACM Trans. Graphics23, 835–844 (2004). 9. X. Tong, J. Wang, S. Lin, B. Guo, and H. Shum, “Modeling and rendering of quasi-homogeneous materials,” in “ACM Trans. Graphics ”24, 1054–1061 (2005). 10. P. Peers, K. vom Berge, W. Matusik, R. Ramamoorthi, J. Lawrence, S. Rusinkiewicz, and P. Dutré, “A compact factored representation of heterogeneous subsurface scattering,” in “ACM Trans. Graphics” 25, 746–753 (2006). 11. J. Yule, D. Howe, and J. Altman, “The effect of the spread-function of paper on halftone reproduction,” TAPPI Journal7, 337 – 344 (1967). 12. P. G. Engeldrum and B. Pridham, “Application of turbid medium theory to paper spread function measurements,” TAGA Proceedings1, 339–352 (1995). #206540 $15.00 USD Received 17 Feb 2014; revised 28 Mar 2014; accepted 28 Mar 2014; published 7 Apr 2014 (C) 2014 OSA 21 April 2014 | Vol. 22, No. 8 | DOI:10.1364/OE.22.009048 | OPTICS EXPRESS 9048 13. C. Ackermann, H. Praast, and L. Göttsching, “Einfluss von Lichtstreuung und Lichtabsorption auf die Bildwiedergabe gedruckter Rasterflächen,” Tech. Rep. 12395N, AiF (2002). 14. H. Wakeshima and T. Kunishi, “Light scattering in paper and its effect on halftone reproduction,” J. Opt. Soc. Am. 58(2), 272–273 (1968). 15. S. Inoue, N. Tsumura, and Y. Miyake, “Measuring mtf of paper by sinusoidal test pattern projection,” J. Imaging Sci. Technol.41(6), 657–661 (1997). 16. G. L. Rogers, “Measurement of the modulation transfer function of paper,” Appl. Opt. 37, 7235–7240 (1998). 17. J. Arney, C. D. Arney, M. Katsube, and P. G. Engeldrum, “An mtf analysis of papers,” J. Imaging Sci. Technol. 40, 19–25 (1996). 18. K. Happel, M. Walter, P. Urban, and E. Dörsam, “Measuring anisotropic light scatter within graphic arts papers for modeling optical dot gain,” in “IS&T/SID, 18th Color and Imaging Conference,” San Antonio, Texas, 347– 352 (2010). 19. M. Ukishima, H. Kaneko, T. Nakaguchi, N. Tsumura, M. Hauta-Kasari, J. Parkkinen, and Y. Miyake, “A simple method to measure mtf of paper and its application for dot gain analysis,” IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences 92, 3328–3335 (2009). 20. M. Walter, “Methode zur Messung orthotroper Lichtstreuung in grafischen Papieren,” Master’s thesis, Technische Universität Darmstadt (2009). 21. K. Happel, “Led-based light scattering measurements of papers for printing applications,” Ph.D. thesis, Technische Universität Darmstadt, Germany (2011). 22. E. S. 1288, “Standard for Characterization of Image Sensors and Cameras,” European Machine Vision Association (2010). 23. K. Happel, P. Urban, E. Dörsam, and X. Ludewig, “Classifying Papers According to their Light Scatter Properties,” in “Midterm Meeting of the International Colour Association (AIC),” Zurich, Switzerland, 138–141 (2011). 24. B. Jähne, Digital Image Processing , 6th ed (Springer, 2005). 25. F. Berg, “Isotrope Lichtstreuung in Papier Neue Überlegungen zur Kubelka-Munk-Theorie,” Ph.D. thesis, Technische Hochschule Darmstadt (1997). 26. F. Ruckdeschel and O. Hauser, “Yule-nielsen effect in printing: a physical analysis,” Applied Optics, 17, 3376– 3383 (1978). 27. G. Fischer, J. Rodriguez-Giles, and K. R. Scheuter, “Ein physikalisches modell für die beschreibung von lichtstreuprozessen,” Die Farbe 30, 199–220 (1982). 28. S. Gustavson, “Dot gain in colour halftones,” Ph.D. thesis, Sweden (1997). 29. F. P. Callahan , “Light scattering in halftone prints,” J. Opt. Soc. Am. 42(2), 104–105 (1952). 30. A. Murray, “Monochrome reproduction in photoengraving,” Journal of the Franklin Institute, 221, 721–744 (1936). 31. J. Viggiano, “The Color of halftone tints,” in “TAGA Proceedings,” 647–661 (1985). 32. E. Marchand, “Derivation of the point spread function from the line spread function,” J. Opt. Soc. Am. 54(7), 915–919 (1964).


Introduction
When light falls onto a translucent material, a fraction of it penetrates into the material and exits at different locations of the material's surface.This light transport is mostly caused by scattering and can significantly impact the material's visual appearance.
The work described in this paper is mainly motivated by spectral modeling of halftone prints where subsurface light transport within the substrate causes a halo around the ink dots.This is known as optical dot gain because the ink dots appear larger than their physical area coverage [1][2][3][4].Neglecting optical dot gain yields poor spectral predictions of the printout.Most spectral halftone models employ various parameters including empirical parameters to account for subsurface light transport.Parameter adjustment requires many measurements of printed training samples and constrained optimization.Accurately modeling the intrinsic physical phenomena (especially mechanical and optical dot gain) could aid the creation of spectral halftone models with only few parameters requiring a small number of measurements for model adjustment [5].This is crucial for the applicability of a halftone model, particularly for printers employing many inks.Ground truth measurements of light transport within the printing substrate are an important prerequisite for modeling physical phenomena related to the interaction between the printed substrate and incident light.
Due to its impact on the visual appearance of translucent materials, modeling subsurface light transport is also necessary in computer graphics to accurately render such materials on a display.Attempts were also proposed to reproduce materials with specified subsurface light transport by 3D printers.[6,7].A prerequisite for all these applications are accurate measurements of the light transport within the material.Various approaches were proposed for this purpose.All have in common that a structured light pattern is projected onto the material resulting in illuminated and less-or non-illuminated areas: Setups proposed in computer graphics used a single laser beam [8], a laser strip [9] or a regular grid pattern emitted from a DLP projector [10].In the scope of printing, Yule and Nielsen [1] projected the filament of a halogen lamp, Yule et al. [11], Engeldrum and Pridham [12], and Ackermann et al. [13] slided a razor blade into a light path projecting an edge onto the substrate, Wakeshima and Kunishi [14] used a small pencil of light, Inoue et al. [15] used a sinusodial pattern, Rogers [16] and Arney et al. [17] employed a bar target and Ukishima et al. [19] used an octagonal pencil of light by closing the iris of a microscope.
Light transport within the material causes a blurring of any projected pattern.The emitted light distribution is used to determine the so-called point spread function (PSF) characterizing the light transport.It describes the fraction of light that exits the material at position (x, y) in direction ω o when illuminated by an infinitesimally narrow light beam incidenting from direction ω i at position (x 0 , y 0 ).In this paper, we consider only one pair of directions which both coincide with the surface normal n, i.e. ω i = ω o = n.The corresponding measurement geometry is called 0/0.
A system for accurately measuring subsurface light transport must allow for various error sources that may bias the results.These are direct reflection at the material's surface (particularly light scatter due to surface roughness), the modulation transfer function (MTF) of the measuring setup's optical system, inaccurate focusing, non-linearity of the camera sensor, noise (thermal noise, shot noise, read-out noise, etc.) and stray-light.Surface reflectance, MTF of the optical system and inaccurate focusing may cause an additional blurring of the captured light pattern resulting in an over-prediction of subsurface light transport.Noise and stray-light may affect the robustness of methods employed to compute the PSF from the captured light distribution.
For measuring subsurface light transport allowing for all error sources mentioned above, we propose in this paper an enhanced version of the apparatus presented in our conference paper [18].Similar to Ukishima et al. [19], we use a reference material to remove the MTF of the optical system for determining the PSF.We encapsulate the whole measuring setup to avoid external stray light.In contrast to other approaches, we propose a mechanical setup enabling focusing on surfaces with non-uniform microtopography (e.g.unevenness caused by paper fibers).This new apparatus improves the instrument introduced in [18] by using • cross polarizing filters to significantly reduce the impact of surface reflection, • spectrally-narrowband exchangeable LED-based illumination instead of a spectrallybroadband halogen lamp.This avoids focusing errors due to chromatic aberration and allows spectrally-resolved measurements of subsurface light transport, • a highly linear CCD camera with 12-bit signal depth instead of a CMOS RGB camera (only green pixels were used) to increase the bit and spatial resolution of captured images.
Another important contribution of this paper is a novel approach based on quadratic programming for reconstructing PSFs of isotropic light transport which is robust against noise.

Measurement setup
The measurement setup is incorporated into a standard microscope, shown on the right in Fig. 1.It allows to project a defined light pattern onto the sample and features a 0/0 measurement geometry.The latter is advantageous for adjusting the focus of the camera and of the projected pattern.The PSF characterizing the subsurface light transport is determined from the laterally resolved capture of the light emitted from the sample's surface.We divide the measurement setup in three logical units: the illumination unit, the sample unit, and the observation unit.A schematic of these units is presented in Fig. 1 on the left.It is an enhancement of the setups presented in previous work [18,20,21].The essential component of the illumination unit is the light source.A light source with a broadband spectral power distribution (SPD) may adversely affect exact focusing due to chromatic aberration.For our setup, a blue, a green, and a red light-emitting diode (LED) are used respectively.The peak emission wavelengths of the LEDs are at approximately 457 nm, 536 nm, and 632 nm respectively.Fig. 2 shows the SPDs of the three LEDs.They were measured with a spectroradiometer (Konica Minolta CS-1000A).These narrowband light sources feature a good compromise between exact focusing, spectral resolution, and luminous power as well as cost efficiency.The position of the LED can be adjusted relatively to the light path.
Another important part of the illumination unit is the setup to create the light pattern.We use the edge of a razor blade, which can be swiveled into the light path to produce a half-plane pattern.The focus of the blade is adjusted by moving a lens.
The observation unit holds a monochrome 2-megapixel CCD camera (Basler scA1600-14gm) with 12-bit signal depth and 4.4 µm pixel pitch.The signal depth of 12-bit and an adapted exposure time ensure that the signal range of the camera is best exploited.Working in the linear domain of the sensor, gray values are linearly related to captured light intensities.
For image acquisition, the manufacturer's software development kit (SDK) was connected to MATLAB, so that images can be directly accessed and processed there.
The combination of all optical parts in the observation unit will be called the imaging system in the following.It will be treated as one unit.The elements of the imaging system are the CCD camera, the beam splitter, the objective lenses, and other windows of the objective turret and the microscope that lie in the observation path.Capturing a stage micrometer revealed the imaging system's resolution of 0.54µm/Pixel.The sample holding unit consists of a special sample holder [20] that was designed for plane samples.It is attached to the microscope stage and allows for a simple adjustment of the sample focus by moving the microscope stage in vertical direction.For this task, an aluminized first surface mirror (FS mirror), its front surface leveled with the surface of the sample, is used as a reference.When the focus is optimally adjusted, the sample is moved to the position of the reference by sliding both perpendicularly to the optical axis of the microscope.This allows meaningful focusing on rough surfaces (e.g.uncoated paper).Furthermore, the sample holding unit allows the sample to be rotated around the optical axes of the microscope.A schematic of the sample holding unit is shown in Fig. 3.The 0/0 geometry of our setup has the disadvantage that signals reaching the camera are inherently biased by specular reflections from the sample's surface.The use of crosswise oriented polarizing filters reduces this effect to a minimum.The filters are located in the illumination and the observation unit as indicated in Fig. 1.
Another important improvement to similar measurement setups is an optically encapsulated system to avoid external stray-light biasing the captured images.

Deriving the measurement equation
Measuring light transport within a material by the described setup requires the detection of differences in the lateral distribution of the light emitted from the sample.
The sample is illuminated by a certain unknown distribution of irradiance E as a function of of the wavelength λ and two variables (i, j), indicating the location on the sample and the pixels of the camera respectively.
For deriving the measurement equation a series of simplifications were introduced.First, we assume that the surface reflection at the bottom interface of the sample has only negligible effect.This also means that the backing or multiple reflections inside the material do not falsify the measurement.Second, the absorption and fluorescence properties of the sample should be small for the considered range of wavelengths.The surface reflection can be reduced to a negligible amount when using polarizing filters.The exitance from the sample is detected with a CCD camera as a distribution of gray values V(i, j), where (i, j) stand for the i'th row and the j'th column.Linearity of the camera provided, the linear signal model according to EMVA [22] is: where V are the light induced gray values, V total the gray value output of the camera, and V dark the dark signal.The responsivity of the camera Kη is the product of overall system gain K and quantum efficiency η.A denotes the area of the pixel, h is Planck's constant, c the speed of light, E camera the irradiance on the sensor, and t E the exposure time.
Introducing an abbreviating variable the light induced gray value at the position (i, j) can be written as: The irradiance on the sensor is a convolution of the exitance from the sample and the PSF P o representing the optical distortions of the observation unit due to the lenses and windows in the light path.In turn, the exitance from the sample is a convolution of the irradiance on the sample E and its light scattering properties described by the PSF P.This way, the full expression for the irradiance on the camera is: The induced gray values are: Since the values of k, E, P, and P o depend on wavelength, Eq. ( 5) is restricted to cases, where the considered wavelength range is adequately small.
The irradiance on the sample E can be modified by generating a defined light pattern.The assumption is reasonable that the new patterned irradiance E M is simply the product of the original irradiance E(i, j, λ ) with a mask M(i, j), the latter generating the light pattern.As a result, the detected gray values caused by the exitance response of the sample with a masked illumination is: Generally, the mask does not need to be of a specific shape.Bar pattern, sinusoidal pattern, an edge or other shapes can be applied.
Let the irradiance on the sample E be approximately constant with respect to (i, j).Then Eqs. ( 5) and ( 6) simplify to: The ratio of relative gray values V M /V is given by: Now, we introduce an opaque reference material with negligible subsurface light transport.In our case, this is a first surface mirror (FS mirror).For such a material, the PSF P reduces to a dirac function and Eq. ( 9) simplifies to: Combining Eqs. ( 9) and ( 10), we gain a measurement equation, where all unwanted variables are removed: For determining the PSF the following four images have to be acquired: • V ref -is the captured reference using the unmasked illumination • V M,ref -is the captured reference using the masked illumination • V -is the captured sample using the unmasked illumination • V M -is the captured sample using the masked illumination Some restrictions have to be considered, when capturing the images for measurement.First, the acquisition of the images with masked and unmasked illumination should be performed using the same exposure time.This holds for both sample and reference.Second, the images of the reference (masked and unmasked) must be taken without polarizing filters, since the reference is assumed to have perfect surface reflection and polarizing filters would eliminate the captured reflectance.The sample images require the use of polarizing filters.Third, samples with high absorption in the wavelength range of the irradiance produce almost no signal and cannot be measured.Fluorescent components in the sample can alter the results.For good results, absorption and fluorescence should be negligible.
The image of the edge will always be influenced by Fresnel diffraction at the razor blade and distortions inside the microscope's windows and lenses.In order to not impair the measurements by such optical distortions we introduce the reference material enabling us to remove them from the measurement Eq. (11).

Computing the PSF for isotropic light transport
Anisotropic light transport (for the considered 0/0 measurement geometry) can easily be detected by rotating the sample around the optical axis of the camera.The sample holder of our setup allows this without removing the sample (see Fig. 3).If for each rotation angle the captured edge image V M is approximately similar, the light transport can be classified as isotropic, i.e.P(x, y) = Ṕ( x 2 + y 2 ) = Ṕ(r).Many materials possess isotropic light transport, for instance, various paper types [23].
In this case, a 1D version of the measurement equation can be obtained by averaging along the j-direction.This is only reasonable, if the mask is approximately constant in that direction, i.e. the projected edge is aligned with the j-direction.Introducing the line spread function (LSF) L (see Eq. ( 26)), we get the 1D measurement equation: where are called edge spread functions (ESFs) and n is the number of pixels in j-direction.A positive byproduct of the 1D measurement equation is that the impact of noise is reduced by averaging n pixels for each position i.The PSF P can directly be computed from an LSF L using Eq. ( 27).
One approach for calculating the LSF from Eq. ( 12) is direct deconvolution in the Fourier domain: where F is the Fourier transform and F −1 its inverse.Unfortunately, this method performs poorly because errors in focusing as well as the MTF of the optical system cause F (a ref (i)) to be zero or very small and any noise is amplified by 1/F (a ref (i)) [24].This indicates that solving Eq. ( 12) with respect to L is an ill-posed problem.Direct deconvolution in the Fourier domain was also used by Ukishima et al. [19] to solve a related 2D measurement equation.They set the fraction to zero if the denominator falls below an empirical threshold.
In order to avoid empirical parameters that might bias the results, we transform Eq. ( 12) into a variational problem where m is the number of pixels in the i-direction.In previous work [18], we used predefined functions of one parameter for substituting Ĺ , turning the variational problem into a simple parameter optimization.Some functions proposed in literature [25][26][27][28] produce small error rates for selected materials.Such an approach is simple but relies on the choice of the function that might be well suited for one material but inappropriate for another material, i.e. yielding to large residual errors.
We need to relax the stiff constraint of using predefined functions by incorporating more general properties an LSF must possess.By transforming the variational ( 16) into a quadratic programming problem we are able to incorporate suitable a priori knowledge on LSFs by equality and inequality constraints: where A ref is the convolution matrix satisfying A ref Ĺ = a ref * L .Since PSFs are positive and for isotropic materials symmetric, LSFs must also be positive and symmetric according to Eq. ( 26).This is ensured by constraint ( 18) and ( 19) where i 0 is chosen at the point of symmetry.A reasonable assumption for PSFs is that they are monotonically decreasing with increasing radius.It can be shown using Eq. ( 26) that this yields monotonically decreasing LSFs for i 0 + i, i > 0, which is forced by constraint (20).One property of LSFs is unimodality which is ensured by combining constraints (19) and (20).The Beer-Lambert law suggests an exponential decline of light transport with increasing distance from the incident position.This applies only for non-turbid materials and subsurface light scattering may cause deviations from exponential decline.For incorporating a constraint on the rate of the light transport's decay, we must formulate a less strong assumption that is applicable to most materials.For this purpose, we assume that LSFs are decreasing at a decreasing rate from the symmetry point.This implies that an LSF is convex (also called concave upward) for i 0 + i, i > 0, or its second derivative is larger or equal to zero as forced by constraint (21).We would like to emphasize that this is our weakest assumption and we do not have evidence that all LSFs must possess this property.In the case of large residual errors, this constraint can be omitted or further relaxed.
Since absorption coefficients as well as all geometry dependent effects are canceled out by the fractions shown in Eq. ( 11), the LSF describes only how incident light is distributed by subsurface light transport.Therefore, the sum of all LSF elements must be one, which is ensured by constraint (22).
The PSF P can directly be computed from the resulting LSF L by numerically evaluating the term in Eq. ( 27).For the papers, measurements were conducted for three specimen cut from different positions of each sample sheet.The opal glass was measured three times, while each time the sample was removed and replaced into the measurement setup.The measured values were than averaged for further processing.
For each sample, the ESF was captured at six angles (0°, 60°, . . ., 300°).At each angle, nine images were taken to prevent deviations due to temporal noise of the light source or the camera.
Since all captured ESFs are almost invariant to the rotation angle, the materials can be assumed to have isotropic PSFs for the 0/0 geometry.Fig. 4 (top left) shows the measured ESFs of the reference and sample materials for the green LED.The corresponding LSFs computed by the proposed quadratic programming approach are shown in Fig. 4 (top right).And the PSFs obtained by numerically evaluating Eq. ( 27) are shown in Fig. 4 (bottom left and right).
Figure 4 (bottom right) illustrates that the average distance of subsurface light transport is largest for inkjet paper followed by opal glass, coated and uncoated paper.Note that all LSFs and PSFs are normalized to integrate to one.Figure 5 shows LSFs of the opal glass for the employed LEDs.It can be seen that light transport depends on wavelendth.However, for the investigated materials such differences are rather small compared to inter-material differences.
Table 2 shows the root-mean-square errors (RMSEs) of the computed LSFs and PSFs: The low error rates indicate that the assumptions used by the quadratic programming approach as well as the computational accuracy are reasonable for the tested materials.Figure 6 shows the measured and reconstructed ESF with the largest error.

Validating the measurements
Quantifying the accuracy of the obtained PSFs is difficult, because ground truth data is not available.Therefore, we only validate whether the prediction performance of a physically-based spectral printer model employing the measured PSF is equivalent to a common empirical model -which requires more spectral measurements for parameter adjustment.
For this, we employ a slightly modified optical dot gain model proposed by Gustavson [28] based on the work of Callahan [29].This model was developed to predict the reflectance spectrum of printed halftones as measured by a 45°/0°spectrophotometer with an aperture A. It assumes that the PSF P of the substrate is known a priori.The model's concept is illustrated   following form (convolution, multiplication and integration are performed wavelength-wise): with Vol(A) = A 1 dz.
To validate our PSF measurements we compare the model predictions with real spectral measurements.For this, a halftone film was prepared containing regular screens with nominal area coverages of 20%, 30%, 40% and 50%.We used conventional amplitude modulated (AM) screens of 40 lines/cm frequency.The film was placed in full contact on top of the material and the whole stack was measured by a Techkon SpectroDens spectrodensitometer obtaining the reflectance R m .The dot gain model allows us to predict R m by knowing: Note that Eq. ( 25) is a simplified model assuming no light transport within the film and neglecting any interreflections at the material -film and -air interfaces due to different refractive indices.It is likely that this introduces additional light transport.Furthermore, the effective area coverage as well as dot-shape of the halftone screen may deviate from our assumptions and the usage of a binary halftone mask (suggesting a uniform transmittance within each dot-area) might only be a rough approximation of the real transmittance-distribution.All these quantities are difficult to estimate limiting the suitability of the model to be used for deducing the PSF from the model equation.
To account for this bias, which we assume to be almost independent of the investigated material, the effective area coverages of the halftones were adjusted to minimize the spectral RMS differences between model predictions and measurements.We used StoraEnso's Lu-miSilk 250 g/m 2 matte coated paper whose PSF was measured by our setup for this purpose.The resulting area coverages, 21%, 32%, 43% and 55%, were then used to create the halftone mask employed by the reflectance model for all investigated materials.Note that optical dot gain caused by the additional light transport mentioned above is intrinsically included in these area coverages.
Figure 8 shows the spectral RMSE of the predictions including and excluding PSF filtering.Excluding PSF filtering in Eq. ( 25) disregards any light transport in the substrate.Hence, optical dot gain is not considered.Our results show in this case that predicted reflectances are always higher than the measured ones.By considering PSF filtering (and thus optical dot gain), the accuracy of the model increases drastically.Error rates are almost similar in magnitude to those of the widely used empirical Yule-Nielsen corrected spectral Murray-Davies printer  model [1,30,31] whose parameters were adjusted to the measured test reflectances.Note that this information was not used by model (25).This is the main result of our validation experiment and indicates that the accuracy of the measured PSFs is reasonable for the purpose of spectral printer modeling.
It is worth mentioning that PSFs of the investigated materials are also well described by the modified exponential function proposed by Gustavson [28].Employing these PSFs, model (25) shows a similar prediction performance as for PSFs reconstructed by quadratic programming.

Conclusions
We presented a microscope-based setup for measuring subsurface light transport.It projects a half-plane light pattern onto the material's surface and captures the emitted light distribution allowing to deduce the light transport's Point Spread Function (PSF).In contrast to related approaches, the setup comprises a reference material to cancel out biasing effects caused by the optical components of the microscope and a special sample holding unit allowing to adjust the focus on rough surfaces (e.g.uncoated paper).Crosswise polarizing filters in front of camera and illumination are used to minimize unwanted signals caused by surface reflectance.Exchangeable spectrally-narrowband LED illumination is employed to avoid focusing errors due to chromatic aberration and to allow spectrally resolved PSF measurements.
We proposed an approach to compute the PSF for isotropic light transport using quadratic programming.This approach is robust to noise and ensures distinct properties the corresponding Line Spread Function (LSF) must likely possess, such as positivity, symmetry, unimodality and decreasing with an decreasing rate from the point of symmetry.
Finally, we show that the accuracy of the obtained PSFs enables simple physically-based spectral printer models to achieve a prediction performance almost similar in magnitude to empirical models whose parameters were fitted to the test samples.

Transformation equations between point and line spread functions
For an isotropic PSF P and the corresponding LSF L the following equations apply [32]: A ′ (v 2 + r 2 ) dv ( 27) where the origin (x, y) = (0, 0) is chosen at the point of symmetry and the prime indicates differentiation with respect to t = v 2 + r 2 .

Fig. 2 .
Fig. 2. Relative spectral power distribution of blue, green, and red LED.

Fig. 3 .
Fig. 3. Sample holding unit (schematic and CAD rendering): The focus is adjusted using the reference (first surface mirror) and is maintained for the sample material by mechanical shifting.The unit allows rotating the sample around the optical axis of the microscope to detect anisotropic subsurface light transport.

Fig. 4 .
Fig. 4. ESFs (top left), LSFs (top right), PSFs (bottom left) and cutout (bottom right) for the green LED.The circles in the cutout (bottom right) indicate the average distance of subsurface light transport.Resolution: 0.54µm/Pixel.

Fig. 8 .
Fig. 8. Spectral RMSE of model predictions including and excluding PSF filtering.The red stars indicate the spectral RMSEs of the Yule-Nielsen modified Murray-Davies model fitted to the test samples (ground truth spectral measurements).

Table 2 .
RMSE a and RMSE b .