Fiber orientation is an important structural property in paper and other fibrous materials. In this study we explore the relation between light scattering and in-plane fiber orientation in paper sheets. Light diffusion from a focused light source is simulated using a Monte Carlo technique where parameters describing the paper micro-structure were determined from 3D x-ray computed tomography images. Measurements and simulations on both spatially resolved reflectance and transmittance light scattering patterns show an elliptical shape where the main axis is aligned towards the fiber orientation. Good qualitative agreement was found at low intensities and the results indicate that fiber orientation in thin fiber-based materials can be determined using spatially resolved reflectance or transmittance.
© 2014 Optical Society of America
Deeper knowledge of light propagation in fibrous structures is fundamental in order to understand the optical appearance as well as for industrial applications such as fiber alignment measurements. The anisotropic behavior of light scattering in paper has been known for a long time even though the underlying physics is not well investigated, due to the difficulties of solving Maxwell’s equations for such complex materials. Work on diffraction from paper sheets have shown that transmitted light scatters in an elliptical shape where the major axis is aligned perpendicular to the fiber orientation main axis . Thicker sheets, where multiple scattering dominates, can instead scatter light in an elliptical shape where the major axis is aligned towards the main fiber orientation direction . It is believed that the ellipticity for optically thicker materials originates from either multiple scattering by fibers or from light guided inside the fibers.
Paper can be regarded as a three-dimensional but mostly planar stochastic network of interconnected fibers. For light propagations simulations in such materials, it is possible to describe the cellulosic fibers as long cylinders. By using an analytical solution to Maxwell’s equations for scattering by an infinitely long cylinder  it is possible to take both the wave nature of light and structural properties of the fiber network into consideration. It was originally suggested by Kienle et al. that phase functions from infinitely long cylinders can be used to describe multiple scattering in fiber-based materials [4,5]. They based their simulations on the Monte Carlo (MC) method which is often used to solve the radiative transfer equation numerically. Similar models have since then successfully been used to model light scattering in anisotropic structures such as biological tissue [6–8], softwood  and textile . We have previously used a similar model to compare scattering by spherical and cylindrical particles  as well modeling anisotropic scattering in paper [12, 13].
The focus of this work is to investigate the anisotropic light scattering effects in paper for both spatially resolved reflectance and transmittance. Light scattering patterns from a set of paper samples are obtained experimentally from papers with different degrees of in-plane fiber alignment. Experimental results are compared to simulations where multiple scattering and directional dependent scattering is considered. The scatterers, or fibers, are described as infinitely long cylinders and their orientation has been estimated from x-ray computed tomography images.
A set of paper sheets was made by MoRe Research, Örnsköldsvik, Sweden, where fiber material, area density and thickness was kept constant and only the fiber alignment was varied. The degree of fiber alignment was varied from randomly oriented fibers, i.e., isotropic structural properties, to high fiber orientation along one axis. The sheets were made from unbleached chemical pulp from Swedish pine and in order to obtain a purely fibrous material, the pulp was filtered to remove any fine material. Furthermore, the sheets were not calendared or compressed in any way which leaves the fibers at their natural cross-section and therefore we assume that the fibers will retain a circular shape. The average sheet thickness is 0.19 mm and the average area density 102.5 g/m2, which correspond to a sheet density of ρs = 539.5 kg/m3. This in turn corresponds to a porosity of 64% if the density for cellulose is assumed to be ρc = 1500.0 kg/m3. Measurements on the tensile strength of the paper sheets, which closely relate to the fiber orientation, where also provided by the manufacturer. Table 1 shows the differences in tensile strength in the machine-direction (MD) and cross-direction (CD) along with the ratio between them. The paper sheets span from nearly isotropic fiber orientation (sample 1) to highly aligned fiber orientation (sample 5).
2.2. Fiber orientation estimation using x-ray computed tomography
To be able to correlate the simulations to the experimental results we need to determine the papers three-dimensional micro-strucure. In the presented work we obtained three-dimensional images of the paper by using x-ray computed tomography (CT) where x-ray absorption is directly related to the local density. The CT-images were captured using a Bruker Skyscan 1172 at the division of applied mechanics, Uppsala University. The voxel size was 2 × 2 × 2 μm3 which enabled us to resolve the spatial extent of single fibers as well as their lumen.
Samples from each paper sheet were cut to roughly 3 × 10 mm2 and then folded into cubes to fit the field of view of the CT-scanner. After imaging and reconstruction, the orientation in the samples was estimated. The procedure is based on the image gradients, which in turn are averaged using the structure tensor approach. Such methodology has previously shown to be suitable for cellulose fibers . For each voxel, the orientation is “at last” found as the direction of the eigenvector corresponding to the smaller eigenvalue of the structure tensor. The orientation is finally summarized for regions of about 1.3 × 1.3 × 0.1 mm3.
The distribution functions of the fiber orientation can be seen in Fig. 1(a) and a sample image from the reconstruction can be seen in Fig. 1(b) where pseudo-coloring is used to visualize the estimated fiber orientations.
2.3. Light scattering measurements
To record light scattering patterns we used a 10 mW 633 nm HeNe laser source focused on the surface of the paper using a lens system, see Fig. 2. For both the transmittance and reflectance measurements a numerical aperture of NA = 0.005 was used for the incident laser source. The focal length was set to 200 mm, which provides enough room for a beam splitter in the reflectance measurement along with a theoretically calculated beam waist of 81 μm and a depth of focus of 16 mm. The spatially resolved transmittance and reflectance is detected using a 16 bit CMOS camera, a PCO Edge 4.2 with a Micronikkor 55 mm lens. The exposure time was adjusted for each measurement point to utilize the full dynamic range of the camera and each recorded image was normalized with the maximum pixel value. The intensity patterns were measured on each of the five paper sheets by capturing 12 points along a 12 mm line which was placed at six arbitrary positions, totaling 72 measurement points for each paper sheet.
2.4. Monte Carlo simulations
The simulation tool used in this work is based on the Monte Carlo method. It is a stochastic numerical method based on the assumption that the electromagnetic field can be regarded as a large number particles that move between scattering events until they all have left the medium or been absorbed. Phase functions from an analytical solution of scattering by infinite cylinders are utilized to describe the scattering behavior of the fiber structure. A detailed description of the implementation of the MC simulation technique with infinite cylinders as scatterers can be found in Yun et al. .
A cylinder illuminated at perpendicular incidence will scatter light in the shape of a disc around the longitudinal axis of the cylinder. Oblique incident light scatters in the shape of a cone around the longitudinal axis of the cylinder and the phase function describes the intensity distribution in these discs or cones. The cones half-angle toward the longitudinal axis of the cylinder equals that of the incident light . This behavior is the origin of a directional dependency for multiple scattering observed in structures where the fibers are aligned .
The sequential random walk between scattering events follow the Beer-Lambert law and directly relate to the concentration of scatterers in the medium. As the paper sheets where custom made, with a very low amount of small scattering particles, we assume that the only scatterer is the cellulose fibers. Hence, the cylinder density is calculated as16]. This would correspond to a cylinder density of Ca = 1.41 · 109 m2 which is consistent with the values for Ca achieved in our previous work  where Ca was estimated iteratively using a Gauss-Newton method to match total transmission and reflectance values.
The refractive index of the fibers was assumed to be n = 1.55  and the background medium air is n = 1.0. The parameters used to calculate the phase functions and scattering efficiency Qs of the cylinders were the diameter d = 35.0 μm, the refractive index missmatch n = 1.55 and the wavelength λ = 633 nm. The scattering coefficient is computed asFig. 1(a). The main fiber axis, ϕ = 0°, was set to follow the x-direction. The distribution in the depth, z-direction, was modeled following a Gaussian distribution with standard deviation σz = 10° to relate to the interconnectivity between fiber layers. Further it is assumed that each fiber scatters light independently meaning that the mean free path between scatterers is assumed to be much longer that the wavelength.
The obtained phase functions and cylinder density Ca together with the fiber alignment distribution functions obtained from the CT-images are then used to simulate the spatial intensity for each of the measured samples. The light source was modeled as a circular Gaussian profile with width 80 μm and the thickness was set to 0.19 mm. Only scattered particles with a propagation direction within 1.2° of the z-direction are detected. This value is close to the angle detected by the camera which is estimated to be between 1.4–2° depending on the aperture.
2.5. Ellipticity of iso-contour patterns
The measured and simulated iso-intensity curves are evaluated by fitting intensity levels into ellipses using the least squares criterion. The equation of an ellipse isFig 3. The ratio between major axis, a and the minor axis, b will be compared to the MD/CD-ratio of the tensile strength.
Another important property in for example paper is the amount of fiber flocs, or small-scale area density variations . This can for example be measured by studying the amount of transmitted light through the paper sheet . As the exposure time of the camera was adjusted for each measurement image, we chose to consider the area of the ellipses, Aellipse = πab, at the different relative intensity levels. The area of the ellipses thus provides an indication of the optical thickness of the material.
Figure 4 shows the iso-intensity curves for both reflected and transmitted light for the paper sheet with the highest degree of fiber alignment, sample 5. Note that the main fiber direction, or machine-direction, is referred to as the x-direction and the cross-direction is referred to as the y-direction. An overall similar characteristic behavior can be observed in the scattering patterns for both the measurements and MC simulations. At high relative intensity levels the measured iso-intensity patterns exhibit a circular shape, while simulations have an elliptical shape elongated in the y-direction, perpendicular to the main fiber direction. The simulated elongation in the y-direction originates from incident light being scattered by cylinders whose orientation is distributed toward the x-direction. When the incident light is scattered for the first time the angle between the cylinder and light propagating direction has a high probability of being close to ζ = 90°. So when the cylinders have an alignment toward the x-direction, scattering will be dominant in the y–z-plane due to the scattering characteristics of a cylinder. At lower intensity levels an increasing elliptical elongation toward the x-direction is observed in both the experimental and simulated intensity patterns. Simulated multiple scattered light propagates foremost in the direction of the main fiber direction, x-direction. This is due to light scattering in the shape of a cone around the cylinder. This aligns the light propagation along the cylinder axis, which in turn produces the elliptical iso-intensity patterns elongated in the x-direction at low relative intensity levels. Even though it is not specifically shown here it was found that lower values of σz < 10° will enhance the elliptical shape in the y-direction at high intensity levels. This in turn slightly shifts the elliptical shape in the x-direction to occur at even lower intensity levels. For values σz > 10° the elliptical shape in both the y-direction and x-direction will decrease. This is obvious as lateral diffusion will decrease if the alignment of the cylinders are not in the plane, i.e., the scattering cones will not sustain propagation in the x–y-plane. While this parameter does affect the appearance of the scattering pattern it still does not do so to a very high degree and for values close to σz = 10°.
In order to compare experimentally determined and simulated scattering patterns more quantitatively the spatially resolved transmittance and reflectance along the x- and y-directions from the center of the light source is show in Fig. 5. The experimental curves are averaged over all 72 measurement points and both the experimental and simulated reflectance and transmittance curves intensity maxima were normalized to enable comparisons with the lower intensity levels. The intensity of the measurements and simulations do not agree well at high relative intensity levels close to the light source. At small distances from the center of the source, between 0 and 0.2 mm, the reflectance determined from simulations Rs is smaller than that of the corresponding experiment Rm. The results on the spatially resolved reflectance, Fig. 5(a), are very similar to work done by Kienle et. al. . Note that their work considered solid wood (softwood) which has a much more prominent fiber alignment where the spatially resolved reflectance spanned several millimeters whereas the work presented here is limited to distances below one millimeter. Their simulations also predicted good agreement to measurements at lower intensity levels even though there were slight differences close to the laser source. They argued that this is due to reflections from the surface which was not considered in the simulations as the topology of that surface was unknown. This can possibly apply to this work as well, however the distance over which the theoretical and experimental data does not fit indicate that the diameter of the incident light source was larger than our estimations. Figure 5(b)) show the relationship between measured and simulated results for the transmitted intensity. The simulated intensity Ts close to the light source, between 0 and 0.2 mm, is a lot larger than the measured intensity Tm. A reasonable explanation for this large difference is that the wavefront of the incident light source is not considered in the simulations. The small numerical aperture of the focusing lens induces small oblique angles for the incident photons, this would in turn greatly reduce the number of detected ballistic photons. Note also that in the measured transmittance in Fig. 5(b) it is possible to see indications of a elliptical shape in the y-direction at high intensity levels, between 0.1 and 0.2 mm from the source. It is however not nearly as apparent as in the simulated transmittance. The scattered spatial intensity at distance greater 0.2 mm from the center of the incident light source was found to agree very well between simulations and measurements. This indicates that the model correctly describes multiple scattering in the fiber structure.
Figure 6(a) shows the calculated area from the ellipse fitting in relation to the tensile strength MD/CD-ratios for the five different paper sheets. The center value is the average of 72 measurement points and the error bars are two standard deviations wide for each of the paper sheets. The area of the chosen intensity levels fluctuates a lot more for transmitted light than the reflected area. This is most likely due to transmitted light being more affected by thickness variations in the paper sheets. This was confirmed by measurements on two paper sheets stacked on top of each other (not specifically shown) where the size of the transmittance light intensity patterns increased in size while no apparent difference was found in the size of the reflectance light intensity patterns.
Figure 6(b) shows how the ellipse elongation angle θ vary in relation to the x-direction. This provides information about variations of the fiber alignment in the structure. A decrease in fluctuations of θ can be observed for samples with higher degree of fiber alignment. This is reasonable since a smaller elliptical elongation of the scattered intensity pattern is more sensitive when the ratio a/b is close to unity. This also partly explains why the standard deviation is lower in the transmittance measurements. Another explanation is that the reflectance measurements are more localized, fewer fibers scatter the light in comparison to the transmitted measurements.
The measured elliptical shape of the iso-intensity curves for each of the different paper sheets was found to show a strong correlation to the fiber alignment. Figure 7 show the measured and simulated a/b ratios in relation to the tensile strength MD/CD ratios. The analyzed relative intensity levels were 0.001 for the iso-contours in reflectance Fig. 7(a) and 0.01 for transmittance Fig. 7(b). The corresponding intensity levels obtained from simulations are different due to the differences in intensity near the source between simulations and experiments. This is most apparent in the transmittance measurements, however, the experimentally determined elliptical scattering patterns at relative intensity level of 0.01 and the simulated at relative intensity 0.003 have nearly the same area. The overall behavior of the experiments and simulations is very similar, a close to proportional relationship between tensile strength and elliptical shape is observed. Both measurements and simulations based on CT-images found little difference between the two samples with the highest degree of fiber alignment. This indicates that the tensile strength measurement of sample 4 possibly can be slightly off. Note however that the CT-images were restricted to small measurement volumes and can therefore posses a fairly large error in relation to the fiber orientation of the whole paper.
4. Discussion and conclusion
Measurements on spatially resolved transmittance and reflectance patterns on paper sheets were found to have a strong anisotropic shape elongated toward the main fiber orientation at low intensity levels. The experimental results were compared to Monte Carlo simulations that considered the micro-structure of the fiber material. In the theoretical model we consider fibers to be described by long cylinders distributed in the structure following probability distributions estimated using CT-images. By using the average diameter of a pine fiber and the refractive index of cellulose we calculated the scattering behavior of the cylinders with an analytical solution of Maxwell’s equations. The only free parameters, that were not estimated in any way, are the fiber distribution in the thickness direction (z-direction) which was modeled as a Gaussian distribution following a standard deviation with σz = 10° and the absorption coefficient which was set to μa = 0.
Good agreement between experimental and simulated results was found for diffused light sufficiently far from the laser source. Large differences can however be observed close to the light source which at least partly are due to the profile and width of the incident laser beam. Measurements on the ellipticity of relative intensity patterns at level 0.01 in transmittance and 0.001 in reflectance show a strong elongation toward the main fiber axis. The iso-intensity contours followed the same overall characteristics in anisotropy both in measurements and simulations far from the laser source. These results indicate that anisotropic light diffusion in paper originates from multiple scattering between fibers and not from light guiding within the hollow lumen of the cellulose fibers. This has also been suggested for other fiber-based materials such as biomedical tissue  and wood . The results can be considered good as paper is a relatively complex material and the model is relatively simple where for example effects from contactpoints between fibers and size variations of the fibers are ignored. Further it must also be mentioned that fibers are not ideal homogeneous cylinders, they can be curved, have irregular shapes or surface roughness likely to cause some of the scattering to occur outside the theoretical scattering cone. This is a reasonable explanation to why the elliptical intensity contour patterns in the y-direction are harder to observe in the measurements. Note that the paper sheets in this work had a porosity of about 64% and contained very low amounts of small scattering particles. The addition of small particles or other filler materials is expected to reduce the anisotropic diffusion observed in this work.
The small variation in area of the elliptical scattering pattern for the reflected light indicates that most scattering occur in the layers of fibers closest to the light source. This due to the diffused reflected light propagating foremost in the fiber layers close to the light source while the transmitted light propagates through all layers. Simultaneous measurements of both spatially resolved transmittance and reflectance from the same laser source therefore has a potential to detect fiber orientation variations between the top and bottom layers of the sheets. This would apply for both the ellipse axis ratio a/b and the angle between the ellipse major axis and the machine-direction, θ. Another extension to such a technique could be to use the degree of depolarization in reflectance measurements to retrieve information about variations in the different layers.
In conclusion, we have shown a strong relation between tensile strength and anisotropic shape in spatially resolved reflectance and transmittance patterns at small relative intensity levels. Both measurements and simulations agree well for multiple scattered light even though differences close to the light source were found. These findings indicate that parameters such as fiber orientation and variations in fiber orientation can be measured on-line by monitoring scattered light in various fiber-based materials using a camera with a high dynamic range. Also, as the reflected spatial intensity only starts to show the elliptical shape at very low relative intensity level, none or very little anisotropic influence should be expected for optical dot gain in paper. Further, the method of estimating the fiber orientation distribution functions from CT-images proved excellent to get valid input parameters for the simulations. A suggestion for a possible extension of the present work is to use CT-images to estimate fiber dimensions, cylinder density and fiber alignment in the thickness direction to validate light scattering models even further. It is also of interest to investigate the anisotropic light diffusion in fiber-based materials mixed with smaller particles.
This work was supported by grants from The Kempe Foundations as well as from Process IT Innovations at LTU. Thanks also to Mr. Thomas Joffre, Division of Applied Mechanics, Uppsala University for lending us the CT-scanner.
References and links
1. C. F. Yang, C. M. Crosby, A. R. K. Eusufzai, and R. E. Mark, “Determination of paper sheet fiber orientation distributions by a laser optical diffraction method,” J. Appl. Polym. Sci. 34, 1145–1157 (1987). [CrossRef]
2. W. Blecha and H. Kent, “On-line fiber orientation distribution measurement,” (1990). US Patent 4,955,720.
3. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley and Sons, New York, 1983).
7. A. Kienle, C. Wetzel, A. Bassi, D. Comelli, P. Taroni, and A. Pifferi, “Determination of the optical properties of anisotropic biological media using an isotropic diffusion model,” J. Biomed. Opt.12(2007). [CrossRef] [PubMed]
11. T. Linder and T. Löfqvist, “Monte Carlo simulation of photon transport in a randomly oriented sphere-cylinder scattering medium,” Appl. Phys. B. 105, 659–664 (2011). [CrossRef]
12. T. Linder and T. Löfqvist, “Anisotropic light propagation in paper,” Nord. Pulp Paper Res. J. 27, 500–506 (2012). [CrossRef]
14. M. Axelsson, “Estimating 3d fibre orientation in volume images,” in “Pattern Recognition, 2008. ICPR 2008. 19th International Conference on,” (2008), pp. 1–4.
16. C. Fellers and B. Norman, Pappersteknik, 3rd ed. (Department of Pulp and Paper Chemistry and Technology, Royal Institute of Technology, Stockholm, 1996).
17. K. Saarinen and K. Muinonen, “Light scattering by wood fibers,” Appl. Opt. 40, 5064–5077 (2001). [CrossRef]
18. M. Deng and C. T. J. Dodson, Paper: An Engineered Stochastic Structure (Tappi Press, Atlanta, GA, 1994).