An experimental set-up was devised to record the transmission of red and green HeNe lasers through different types of paper. The measured data was compared with data obtained using the Henyey-Greenstein function (often employed in paper optics models to represent the bulk scattering of material samples) and data obtained using an alternative exponentiated cosine function. The comparisons are used to qualitatively assess the degree of fidelity of the bulk scattering approximations provided by both functions.
©2008 Optical Society of America
In the field of grading paper, many appearance related attributes are evaluated using optical techniques. Properties like brightness and opacity can be determined through the measurement of light reflection and transmission responses of material samples . Additionally, structural characteristics like roughness and thickness can be evaluated using optical styli  and lasers . Recently, computer simulations are also being used to predict paper appearance attributes before manufacturing to mitigate production costs. For example, the effect of dyes and pigments on color and brightness can be predicted before physical samples are produced.
It is important to note that different approaches can be used to simulate paper optical properties, and no single approach is superior in all cases. While stochastic approaches (e.g., Monte Carlo based models [4, 5]) are known for their flexibility, deterministic approaches (e.g., Kubelka-Munk and discrete ordinate based models ) usually require less computational time. Both approaches, however, employ functions to describe the bulk scattering of material samples.
In practice, the choice of a function to describe the bulk scattering of a given material is associated with data constraints. For example, after Bruls and van der Leun  suggested that their measurements of the scattering profile of skin tissues could be approximated by curves obtained using the Henyey-Greenstein function (henceforth referred to as HG) , it started to be employed in studies involving light transport in tissues for which measured scattering data is scarce, either by plugging it into the radiative transfer equation , or using it to derive warping functions within Monte Carlo modeling frameworks .
One technique often employed in Monte Carlo integration is importance sampling . If the integrand is a product of two functions and one of them is known, this function, called importance function or probability density function (PDF), can be used to guide the sampling strategy. In Monte Carlo simulations involving ray optics techniques, the direction of a light ray interacting with a material’s internal structures is given by warping functions derived from the integration of the PDF. In short, these functions provide the scattering angles of the propagated ray.
Clearly the selection of an inappropriate PDF to approximate the scattering profile of a material may introduce error in the modeling of its appearance attributes. For example, the HG function is neither based on the mechanistic theory of scattering , nor was it originally proposed to represent the bulk scattering of multilayered materials . Furthermore, as stated by Jacques et al. , “the use of the HG function to specify radiant intensity for thicker samples is only descriptive, and should be distinguished from the customary use of the HG function to describe single particle phase function”. Despite these issues, due to its convenient mathematical tractability, it is often employed to approximate the bulk scattering of complex media, including stochastic (e.g., using it as PDF in a Monte Carlo integration ) and deterministic (e.g., expanding it using a Legendre polynomial technique in a discrete ordinate formulation ) frameworks used to model paper optics properties. As stated by Neuman , there is a need for more experimental data that can allow the selection of functions that adequately described the light scattering in paper.
In this work, the fidelity of HG based approximations is investigated. More specifically, data generated using HG based warping functions (HGWF) are compared with scattering data measured for four different paper samples illuminated by two different HeNe lasers, as well as data computed using a warping function derived from an exponentiated cosine function (henceforth referred to as ECWF and EC respectively). Initially, the physical measurement set-up used to collect the actual scattering data as well as the warping function formulations are presented. The closest approximations between the physical data sets and the curves obtained using the warping functions are then reported. Finally, the main conclusions derived from this investigation are summarized and directions for future research are outlined.
2. Experimental methods and materials
2.1. Physical measurement set-up
The set-up used in this investigation was selected because of its low cost/effectiveness ratio. It allows us to obtain reliable data to be used in the qualitative comparisons depicted in this investigation. It is also relatively simple to assemble, which facilitates the reproducibility of the results. In fact, similar set-ups based on the use of a CCD camera and lasers have been used in other studies tackling related optical questions. For example, Oksman  used a CCD camera to capture images of paper to determine reflectance under different pressure. Carlsson  used a CCD device to determined the time delay of the laser transmission at different wavelengths and inferred the spectral reflection and transmission of different paper samples. Béland and Bennet  used a CCD camera in conjunction with scatterometry techniques to assess the effect of microroughness on glossy papers. The reader interested in more complex set-ups used to obtain highly accurate scattering measurements is referred to a comprehensive book on this topic .
The raw measurement data were collected as images of laser transmitted through pieces of different paper samples. For each type of paper, multiple images were obtained, each from a different random position from the paper sample. These photos were all taken in a dark room to eliminate contributions introduced by unwanted light sources. The lasers employed in this investigation were a red 1.07 mW low-power HeNe laser at 633 nm and a green 0.5 mW low-power HeNe laser at 543 nm. To prevent saturation in the images, the intensities of the lasers were limited by a neutral density filter placed at the exit window. A Cannon PowerShot S5 IS digital camera was then used to capture images of the transmitted light. The settings used when taking the photographs were 2.48 second shutter speed, f/3:5 lense aperture and 400 ISO. The set-up of the experiment is illustrated in Fig. 1.
The paper samples used in this investigation were a sheet of lined paper from Hilroy white wide ruled letter pads, a sheet of Staples multiuse 20 lb acid free paper, a sheet of Heinz Jordan Acryloil Series 401 acid free illustrator’s paper and a sheet of 140 lb acid free watercolor paper. These paper samples will henceforth be referred to as lined paper, printer paper, sketch paper and water paper respectively. These papers were selected for their differing characteristics with thickness being the most qualitatively apparent. With this characteristic in mind, from the thinnest to thickest, the samples are ordered as follows: lined paper, printer paper, sketch paper and water paper.
Since the primary interest is in the spatial distribution of propagated light, the transmitted laser intensities were radially collected and averaged. The resulting values were then normalized and plotted as a function of the distance (in pixels) from the center of the laser beam. This in effect calibrated the data against the dark frame of the camera, and reduced the effect of any noise that might be present in the images. Figures 2 and 3 show these plots for the red and green lasers respectively. Due to the finite width of the beam, the central portion of the intensity profile, represented as the partition left of the vertical line in the plots, is characterized by relatively small fluctuations. Hence the focus of this investigation is directed toward the decay in intensity near the edges of the laser beam, marked as the region right of the vertical dashed line, which is associated with the most significant changes in the spatial distribution of propagated light. For each laser, the regions are kept with the same boundaries to ensure a fair comparison within each laser’s trials. The plots were then normalized to the maximum and minimum value in these regions. We remark that the results obtained using the two different lasers should be analyzed independently of one another since their corresponding physical measurement set-ups may contain minor inconsistencies due to the handling of the different emitters.
2.2. Warping functions
Probability density functions can be sampled by selecting N random variables ξi(i=1 to N) uniformly distributed over the interval [0;1], and transforming them using an importance sampling procedure . In the case of stochastic ray optics simulations, this procedure consists in the integration of PDFs to obtain the correspondent cumulative distributions, or warping functions, given in terms of the polar and azimuthal (ϕ∈[0,2π)) angles of the propagated (scattered) rays . For this investigation, we selected two PDFs, namely HG and EC.
The HG can be written as :
where g is the asymmetry factor. After performing the transformation (warping) mentioned above, the corresponding polar angle of the propagated rays is given by :
The EC can be written as :
where n is the directional (specular) exponent. Similarly, after performing the transformation (warping) mentioned above, the corresponding polar angle of the propagated rays is given by :
Sample rays (107 rays per simulation instance) were generated using these warping functions, namely HGWF (Eqs. (2) and (5)) and ECWF (Eqs. (4) and (5)), and collected using a virtual goniophotometer  to determine the corresponding intensities.
2.3. Comparison procedure
The simulated data was collected in terms of intensity as a function of angles. A conversion was required from angles to distance in order to compare against the measured data, which was collected in terms of intensity as a function of distance. Additionally, a horizontal scaling factor was introduced to the simulated data in order to ensure the data were of the same magnitude. The maximum value in the simulated data was normalized to 1 as well.
The selection of the parameters n and g involved two stages. Initially, a coarse exploration of their space of values was performed to determine the ranges that yield the closer approximations with the data obtained in the physical measurement set-up. The parameter searching process was then refined within these ranges to obtain the best approximation. At the end, the parameter space explored for n was 1 to 500 with refined steps of 1, and the parameter space explored for g was 0.0001 to 0.9999 in refined steps of 0.0001. The best approximation was the one that yielded the smallest root mean square (RMS) error  defined as:
where a i and b i are the ith value of the two data sets, each with N elements.
Figures 4 and 5 present the graphs depicting the closest agreements obtained between the measured and simulated data sets. These graphs show that none of the approximations completely agrees their corresponding measured data, with the larger deviations occurring in the region near the inflexion point of these curves. The discrepancies between the data sets are further illustrated in the orthogonal projections of the profile curves shown in Figs. 6 and 7. In general, the HGWF curves do not conform to the relatively sharp bend of the measured curves at the normalized distance between 0.2 and 0.5. This is seen as an overall brighter image and more diffuse central disk in the orthogonal plots. On the other hand, the ECWF curves fall off quickly, especially for normalized distances greater than 0.4. Their deviations from the measured curves, however, occur on a scale considerably smaller than the one observed in the HGWF curves.
Comparing the results provided by the two warping functions, one can verify the ECWF data sets present the closer qualitative and quantitative agreement with the actual sets. The later observation is confirmed by the RMS errors presented in Tables 1 and 2.
Also note that the parameters resulting in the best matches do not correlate with the thickness of the paper samples. As the paper samples become thicker, the absolute intensity of the measured data decreases, and the parameters would be expected to decrease as well. However, neither g nor n decreases monotonically across the measurements for the same laser.
It is worth mentioning that the ECWF is approximately 40% more expensive than the HGWF. This figure is based on the elapsed CPU time required for their computation on the SGI Onyx 3200 used in this investigation.
Usually the choice of a function to approximate the bulk scattering of a paper’s internal layers, either in stochastic  or deterministic  simulation frameworks, represents a compromise between accuracy and mathematical tractability. Such a compromise has motivated the use of the HG function in these frameworks. Although the experiments described in this paper were performed on whole samples, they suggest that, in the absence of measured data, the EC function may be a viable alternative for representing the bulk scattering of paper layers, notably for applications requiring a low level of fidelity (e.g., visual inspection of computer generated images of translucent paper). However, for applications requiring a high level of fidelity (e.g., prediction of paper properties using high sensitivity measurement data in conjunction with simulations), the degree of accuracy provided by both approximations may be below acceptable limits.
We remark that the parameters of both functions, n and g, have no direct connection with the structural characteristics of paper. Additionally, the results show that the trend of the parameters could not be predicted based on the qualitative properties of the papers, namely thickness and transmission intensity. This limitation may prevent their selection to be quantitatively guided by the characterization data describing different samples, which, in turn, is likely to affect the predictability of models employing these functions.
One of the strategies to increase the predictability of the current modeling frameworks may involve the development of data driven approaches in which physically measured data can be directly incorporated into the simulations. Such approaches, however, are bound by data scarcity, which is one of the main reasons for the incorporation of bulk scattering approximations into models in the first place. Therefore, there is a need for more comprehensive measurements of paper optical properties that take into account not only the degree of non-uniformity in different types of paper, but also a broad range of measurement and environmental conditions. Such an undertake will likely require a close collaborative effort between industry and academia. Future research initiatives in this area should also include thorough studies on how the errors introduced by these approximations propagate through complex modeling frameworks.
References and links
1. J. Casey, Pulp and Paper, 2nd ed. (Interscience Publishers, Inc., New York, 1961)Vol. 3,.
2. A. Oksman, R. Silvennoinen, K. Peiponen, M. Avikainen, and H. Komulainen, “Reflectance study of compressed paper,” Appl. Spectrosc. 58, 481–485 (2004). [CrossRef]
4. K. Green, L. Lamberg, and K. Lumme, “Stochastic modeling of paper structure and Monte Carlo simulation of light scattering,” Appl. Opt. 39, 4669–4683 (2000). [CrossRef]
6. P. Edström, “Comparison of the DORT2002 radiative transfer solution method and the Kubelka-Munk model,” Nordic Pulp and Paper Research Journal 19, 397–403 (2004). [CrossRef]
8. L. Henyey and J. Greenstein, “Diffuse radiation in galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]
9. S. Jacques, C. Alter, and S. Prahl, “Angular dependence of HeNe laser light scattering by human dermis,” Lasers Life Sci. 1, 309–333 (1987).
10. S. Prahl, M. Keijzer, S. Jacques, and A. Welch, “A Monte Carlo model of light propagation in tissue,” SPIE Institute Series IS 5, 102–111 (1989).
11. J. Hammerley and D. Handscomb, Monte Carlo Methods (Wiley, New York, 1964). [CrossRef]
12. M. Neuman, “Anisotropic reflectance from paper - measurements, simulations and analysis,” Master’s thesis, Department of Physics, UMEA University (2005).
13. M. Béland and J. M. Bennett, “Effect of local microroughness on the gloss uniformity of printed paper surfaces,” Appl. Opt. 39, 2719–2726 (2000). [CrossRef]
14. J. C. Stover, Optical Scattering: Measurement and Analysis (McGraw-Hill, Inc., 1990).
15. G. Baranoski and J. Rokne, Light Interaction with Plants: A Computer Graphics Perspective (Horwood Publishing, Chichester, UK, 2004).
16. G. Baranoski, A. Krishnaswamy, and B. Kimmel, “Increasing the predictability of tissue subsurface scattering simulations,” The Visual Computer 21, 265–278 (2005). [CrossRef]
17. E. Lafortune and Y. D. Willems, “Using the modified Phong reflectance model for physically based rendering,” Tech. rep., Department of Computer Science, K.U. Leuven (1994).
18. P. Edström, “Fast and stable solution method for angle-resolved light scattering simulation,” Tech. Rep. R-0235, Mid Sweden University (2002).
19. A. Krishnaswamy, G. Baranoski, and J. Rokne, “Improving the reliability/cost ratio of goniophotometric comparisons,” Journal of Graphics Tools 9, 1–20 (2004).
20. M. Spiegel, Schaum’s Outline of Theory and Problems of Statistics, 2nd ed. (McGraw-Hill, Inc., 1991).