Evaluation of the orientation distribution of fibers from reflection images of fibrous samples

We consider illumination systems and mathematical algorithms for determination of the anisotropy and topographical features of an illuminated surface from its reflection images. As a particular example we study determination of the fiber orientation of paper surface. We also consider illumination systems with multiple light sources, and introduce optimization algorithms that exploit different spectral bands of these light sources. We show that a system of three light sources, e.g., a blue, green and red LED placed in a regular triangular form, efficiently prevents distortion of the above features. It is also easy to implement in applications, e.g., of the paper industry. We furthermore show that a pentagonal illumination system would be even better. However, in this case a five-colour photographic system and the corresponding light sources would be needed.


Introduction
Wood fibers, mineral fillers and other additives form the basic structure of paper.The properties of paper depend largely on the spatial distribution of fibers.Therefore, it is important to be able to measure and control fiber orientation during, e.g., the paper-making process, and also in the manufacturing of other fiber-containing products.
This work concentrates on testing and optimizing a quantitative tool for fiber orientation analysis, based on optical reflection images of paper and statistical signal processing that involves the curvelet transform.
Fibers form a more or less random network with predominantly planar orientation of fibers.As an off-line measurement, it is possible to study the three-dimensional (3D) fiber structure of the paper by sheet splitting [1], or X-ray tomographic imaging [2,3].Both these methods are, however, slow and do not apply to online measurements.
Charge-coupled device (CCD) cameras allow fast optical reflection imaging of paper samples, although these images predominantly reveal only the orientation of fibers near the sample surface.Fortunately, it is often enough to do this for both sides of the sample.One example of such a case is the gloss variations caused by fibre oriena e-mail: jouni.j.takalo@jyu.fitation streaks in high quality paper coating [4,5].Anoher problem is the the tendency of paper sheets to curl because of the two-sidedness of fiber orientation.Let us discuss the latter example in more detail.
The global two-sidedness of fiber orientation is known to cause curvature in a sheet made of fibers (see [6]).If the orientation angles of the two sides of the sample are significantly different, there usually appears a diagonal (twisted) curl in the paper sheet.A quantitative expression for this kind of curl is given by [1,7].
where K xy is the intensity of the diagonal curl, θ TS − θ BS is the difference in the average fiber orientation angles on the top and bottom sides of the sample, and t is its caliper (thickness).H MD stands for hygroexpansion (expansion or reduction in size due to variation in the moisture content of the sample) in the so-called machine direction (MD).Hygroexpansion in the perpendicular (cross-) direction (CD) is denoted by H CD .Even if the fiber orientation were similar on both sides of the sample, difference in their anisotropies may cause a CD curl, approximated by [1,8].
In equation ( 2), K y is the intensity of the CD curl and A TS −A BS the difference in the anisotropies of the top and

10703-p1
The European Physical Journal Applied Physics bottom sides of the sample.Anisotropy is defined as the ratio of the major-axis to the minor-axis of the orientation ellipse.Expression equation ( 2) assumes that fibers are mainly oriented in the MD direction on both sides of the sheet with roughly an equal strength.
In this work we show how to use the curvelet transform for extracting A TS and A BS effectively from optical reflection images of a paper sample.The reason for using curvelets is that it has been recently shown [9] that the curvelet-based orientation analysis gives much better (more robust) results than the traditional methods, like gradient-and FFT-based orientation analysis methods.Moreover, we study what kind of illumination configuration would be optimal for such an orientation analysis in terms of stability, cost and simplicity.Our results strongly suggest that, to this end, an oblique illumination from five equiangular directions would be the best choice.
Notice that imaging technology is advancing so fast that it may soon be possible to image the whole paper web on-line in an operating paper machine with a resolution that would be high enough for, e.g., the orientation analysis.Also, the methods developed in this study will probably work in many other applications as well.One such example is the determination of the orientation of fibers or nanofibrils in fiber-reinforced composites [10].

The optical measurement system
Optical measurements on paper samples were done by illuminating them with a parallel light beam created by red light emitting diodes (LEDs), and recording the reflected light directly above the sample, see Figure 1.The arm holding the LEDs can be rotated around the sample table.Images were recorded with a Nikon D300 camera equipped with a 105 mm macro objective.

Correction of uneven illumination
Pre-processing of images for the orientation analysis consisted of estimation and removal of uneven illumination caused by the measurement system and optics.The main cause of the uneven illumination here was the discrete system of light sources, which produced a bias in the illumination at low spatial frequencies.Therefore, this kind of unevenness in the illumination can be estimated, e.g., with low-pass filtering or with a low-order polynomial fit to the grey scales of the image.
We estimated this unevennes with a set of successive one-dimensional polynomial fits.The grey scales of each column of the image was first fitted by a second-order polynomial.A new image was constructed from these polynomial fits.The grey scales of each row of this new image was fitted by second order polynomials.The final image constructed from these polynomial fits was then considered to provide an estimate for the unevennes of illumination.Removal of uneven illumination was performed as a point-wise division of the original image with the above estimate.An original image of paper surface, its estimated illumination distribution and the image resulting from removal of this illumination distribution is shown in Figure 2.

Determination of orientation
For evaluation of the fiber orientation we used the method recently introduced by Sampo et al. [9], which is based on directional wavelets, so-called curvelets [11,12].Sampo et al. have shown that curvelet-based orientation analysis is, in most cases, better than traditional methods like the gradient-based (e.g., the structure tensor (ST) method [13,14]) and direct FFT-based [15] methods.The curvelet-based orientation analysis is more robust in particular when the resolution of the image decreases [9].As a comparison of the aforementioned methods we show orientation analysis results by them for a simulated fiber network of known orientation distribution in Figure 3, where we show the simulated network and its orientation distribution as determined by the ST-, curvelet-and FFT-based method (Figs.3b, 3c and 3d, respectively).In all the result figures we also show (by dashed curve) the known orientaton distribution for comparison.Let f be a representation of a grey scale (optical) image of the sample.We can express this representation in terms of basis functions φ b and γ abθ such that Functions φ b are restricted to low frequencies and are therefore not interesting to us in this application.However, functions γ abθ are called curvelets for which we have discretizations for the scale (a), location (b) and angle (θ), which are roughly like a j = 2 −j , j = 1, 2, ...; θ j,l = 2π•l•2 − j/2 , l = 0, 1, ..., 2 j/2 −1 ; and b

10703-p3
The European Physical Journal Applied Physics Fig. 4. Four curvelets with different orientation, location and scale.This image was produced with CurveLab [16].
In Figure 4 four curvelet functions γ abθ are illustrated in both the spatial and frequency domain (zero frequency in the middle).Curvelets of the two largest scales have the same orientation, while the two curvelets of smaller scale have orientations that are different from that of the other two.In the frequency domain curvelets are smooth and compactly supported, while in the spatial domain they decay rapidly in the dominated dimension (along the long axis of its elongated shape) and oscillate in the perpendicular direction.The important difference to ordinary wavelets is that curvelets obey parabolic scaling, meaning that their width to length ratio is a 1/2 , while for wavelets this ratio is independent from scale.Another difference to usual wavelet transforms is that for curvelets, the number of orientations grows when their scale decreases.These features makes curvelets very efficient for approximating/analyzing edges [12].
At small scales a, the magnitudes of inner products f, γ abθ are big for such curvelets γ abθ that are oriented like edges in an image [9].Moreover, curvelets form a tight frame, i.e., . This property suggests that an estimate for the orientation distribution could be expressed in the form: where index set I for the scales depends on the resolution and size of the particles in the image, such that we do not capture any unwanted structures.In our work we used two smallest scales depending on the resolution.If the resolution of angle θ needed to be beyond what CurveLab [16] could provide, it was possible to generate any γ abθ by rotating γ ab0 .For each scale we had a rectangular grid of points as the translation set, with constant C 1 and C 2 [9].

Principles of analysis
Let us first study measurements with only one light source.
In Figure 5 we show four different positions of the source (marked with arrows), and the orientation distributions determined by illuminating a paper sample from these positions.Measurements were made with the setup shown in Figure 1 with an inclination angle of 14 • .Here we prefer cartesian coordinates to polar coordinates that are the traditional way of showing orientation distributions, because in a cartesian presentation differences in the distributions are displayed more clearly.It is evident from this figure that orientation distributions appear to be very different, although they are determined from the same paper sample.Illumination direction seems to affect so that it selects the features that are emphasized in the image: for each direction of the light source, the apparent orientation direction appears rather perpendicular to that of illumination.Notice that when the illumination direction is parallel to the main orientation direction, there is a minimum in that distribution and maximum is at 0 • , i.e., perpendicular to the actually known maximum of 90 • (MD direction).
If the orientation distributions determined for different illumination directions at regular intervals around the circle are combined, the apparent direction-dependent features in them cancel out such that the end result is the true orientation distribution.Even three light sources, placed at regular intervals (of 120 • ) around the circle, are enough to show in a realistic manner the main orientation direction of the sample.Furthermore, measurements for different inclination angles produce the same main orientation.However, it seems that the smaller is the inclination angle the larger is the peak to valley difference in the distribution (peak to valley difference means the same as the major-axis to minor-axis ratio in polar coordinates, i.e., the anisotropy).This is evident from Figure 6, where orientation distributions for three regularly-spaced illumination directions are shown for a varying inclination angle.Normalization of these distributions was made such that the area under each curve was unity (in radians).The inclination angle cannot, however, be too small because of distortions that arise from other features of the surface than those related directly to fibers, e.g., from its topographic features.On the other hand, too large an inclination angle leads to noisy fluctuations in the distribution.There is thus an interval of optimal inclination angles.It is sample dependent so we do not analyze it in more detail here.
In Figure 8 we show the orientation distributions for illumination systems composed of three, four and five light sources (LEDs) placed at regular angular intervals.These distributions were determined such that, after the first exposure, each illumination system was rotated about 10  counterclockwise, and a new image was taken.This procedure was continued until all independent positions of the illumination system were covered.It is evident from this figure that even the illumination system composed of three light sources gives quite reliable and stable (independent of the positioning of the system) results, and that for a system of five light sources, the results are very good.For the system of four light sources, the main orientation stays the same independent of the positioning of the system, although the overall distribution somewhat changes with repositioning.This is probably due to the fact that the angles between the light sources are right, and therefore two of the four light sources emphasize too much the cross-directional orientation.So we suggest to use an illumination system composed of three or five light sources placed at regular angular intervals for a reliable characterization of the sample.For example, an illumination system composed of six regularly placed sources will not provide any advantage over the three-source system as a hexagon can be thought to be composed of two equilateral triangles, and the results would be similar to those for a single triangle.Actually the three-source system could be arranged such that the sources are only 60 • apart (they would cover a half of the circle, the light source in the middle would be placed at the mirror-image position in comparison with its original position).This arrangement of light sources will produce similar distributions as the one based on an equilateral triangle.A schematic drawing of this arrangement is shown in Figure 7c.A similar arrangement can be made of the illumination system composed of five regularly spaced sources, as shown in Figure 7d.In this case two of the light sources will be moved to their mirror-image positions.
Fig. 8. distributions for illumination systems composed of three, four and five regularly spaced light sources at an inclination angle of 14 • .For each system this distribution was determined every time the illumination system was rotated by about 10 • from its previous position.
Notice that not all the light sources in the illumination system should be on simultaneously unless they use different spectral bands (for fast imaging).For white light (or if the light sources use the same spectral band) the orientation distribution is always determined for one light source at the time, and the final distribution is determined as the average over the individual-source distributions.Situation is very similar to determination with a telescope the surface details of the moon.The topographic features one can distinguish depend on the angle of sunlight.Although the full moon is beautiful, the topographic features cannot then be properly distinguished.One should not combine the images for varying lighting conditions (so as to get the 'full moon' image), but the topographic maps as determined in each case separately.When images are combined directly, some topographic features cancel each other, and the topography as determined from the combined image is not so good.

Circular statistics
When analyzing the statistical features of orientation, one has to keep in mind that orientation angle is a circular variable.Circular statistics must be used when we deal with angular variables.There is no unique zero and, in contrast to linear scale, there is no notion of high and low values.In fact circular statistics can be applied to any periodic phenomena, not only to data on angular variables [17][18][19].
We thus need to use circular statistics when analyzing the fiber orientation in, e.g., samples of paper.
The important quantities now include: • Standard deviation.
• Variance, which for a circular distribution does not in general equal the square of the standard deviation.
Variance has a value between zero and one.• Kurtosis, which is a measure of peakedness of the distribution.A large kurtosis indicates a spiky distribution.
has a value between zero and one.The closer is R to one, the more concentrated are the data around the mean.

Materials
The paper sample in our measurements was a sheet of a copy paper with a grammage of 83.9 g/m 2 and a caliper of about 100 μm.It was made of a furnish that contained about 77% short-fiber chemical hardwood pulp and 23% filler (calcium carbonate).

Results
All reflection images were taken with the system shown in Figure 1 (notice that in this figure only three LED sources are visible).First we analyzed the influence of the inclination angle of the illumination (the angle between the light beam and the sample surface).To this end four different inclination angles, 14, 23, 30 and 38 • , were used.The direction of illumination (by one light source) was rotated with an angular increment of about 2 • , and an image was taken in every direction.Altogether 168 images were thus taken of the sample for each inclination angle.From these images we chose subsets that defined the different illumination systems considered (shown in Fig. 7).We then analyzed the possible two-sidedness (different orientation distributions on the the top and bottom sides of the sample) of the paper sample.Here we used the five-light-source (pentagonal) illumination system in the imaging.The anisotropy of orientation was to some extent affected by the inclination angle, and we chose an inclination angle of 14 • for this study.In Figure 9a we show all the five individual (one for each position of a single light source) orientation distributions for the bottom side of the sample.The individual distributions are very different, but the average distribution (shown with a bold line in this figure) has a maximum very near (at 91 • ) the so-called machine direction (at 90 • in the figure), which is very much as expected for a paper manufactured in a paper machine with a wire section (wire-formed paper) [6].Distributions are normalized such that there is unit area under the curve.In Figure 9b we show similar curves for the top side of the sample.The main orientation angle is now 87 • , quite near that of the bottom side, but the anisotropy is much smaller on the top side.
For comparison we show the average orientation distributions for both sides of the sample in Figure 10a.In Figure 10b we show the same distributions in polar coordinates, which is more traditional in paper research.In polar coordinates the orientation distributions have rather an elliptical shape ('orientation ellipse').This shape of the distribution originates from the different speeds of the pulp suspension in the headbox and in the wire section of the paper machine [6].
Although wire forming usually produces a main orientation angle that is near the machine direction (MD), this is not always the case.That is why an illumination system is needed, which is robust enough so as to find the orientation distribution reliably also in the exceptional cases.To this end we analyzed illumination systems composed of regularly spaced (angular spacing) light sources of a varying number.The main goal was to find the number of light sources for given requirements of accuracy.As the method is based on averaging over individual distributions determined for a single light source placed in different directions, we made a dense mesh of (angular) measuring points from which subsets corresponding to different illumination systems could then be chosen.
We thus made 90 measurements (images) on the paper sample with two-degree angular intervals, corresponding to 90 different directions of a single light source.From these images we constructed 90 subsets which describe, e.g., a regularly spaced three-light system.We determined the orientation distributions for single-light illuminations for each of the 90 subsets.For each light direction we can combine the results into two-dimensional orientation diagrams in which orientation intensity (between values 0

10703-p7
The European Physical Journal Applied Physics Images were taken with a two-degree difference in the direction of the individual light source, so there were 90 subsets of the three-light illumination system (corresponding to its rotation by 0 to 178 • with two-degree increments, with respect to the sample).Normalized orientation intensities are described by gray-scale values defined by the gray-scale bar.
and 1) is described by a grey scale, vertical axis is the subset number, and horizontal axis is the angle of orientation.In Figure 11 we show the result of this kind of analysis for the regularly-spaced three-light illumination system.In three of the panels we show the single-lightillumination diagrams and in the last (bottom right) panel the average (combined) result of the three former panels.Although the three-light illumination system can find the overall orientation distribution (last panel) quite well, the diagram shows some inhomogeneity in the variation of grey scale values and, in spite of locating the maximum well, the anisotropy somewhat changes depending on the illumination subset.This is evident in Figure 12a which shows the average orientation diagram of Figure 11d in a three-dimensional representation (orientation intensity described by grey scales in Fig. 11 was made a third dimension).The maximum orientation intensity displayed periodic oscillation as a function of illumination subset.The mean value of the maximum was 0.583 with a standard deviation of 0.016.This means that orientation intensity and anisotropy (the maximum to minimum intensity ratio) somewhat depend on the direction of the illumination system with respect to that of the maximum orientation.
It is important that the intensity of the the light sources are equal in the illumination system.In order to demonstrate the effect of uneven illumination, we show in Figure 12b the result of decreasing in the three-light system the intensity of one light by 40%.The standard deviation of the maximum intensity increases to 0.053 in this case.
In Figure 13 we show the results of our analysis for the regularly-spaced five-light illumination system.In this case we show only the average result, orientation diagrams for individual lights are similar to those for the three-light illumination system above.It is evident that there is less inhomogeneous variation and oscillation in the orientation intensity than in the three-light system (the mean value of the maximum intensity was 0.583 with a standard deviation of 0.007) so orientation intensity and anisotropy did not now depend much on the relative orientation of the illumination system.
We used circular statistics (see Sect. 2.5) to characterize the two-sidedness of the present paper sample.To this end we created ten realisations of 3000 points (angles) between [0, π] of the orientation distributions for the top and bottom sides of the sample, and calculated the means of the characteristic statistical variables for both distributions.The results obtained are shown in Table 1.
The most striking difference in the characteristics shown in Table 1 is in kurtosis, which for the bottom side is about twice as big as for the top side.Variance and standard deviation are both bigger for the top side, as expected from Figure 10.

Discussion
A new method was introduced for estimation of orientation distribution from reflection images of fibrous materials.In this study paper was used to exeplify such materials.This method is based on directional wavelets, so-called curvelets.We demonstrated that, in order to measure the orientation on sample surface, a three-light or five-light illumination system should be used depending on the required accuracy.The three-light system can, e.g., be realized with red, green and blue LED lights, and using three-colour photographs in which the RGB components are treated separately.In this case all light sources could be used simultaneously so that imaging would be a fast process, suitable, e.g., for on-line applications.The five-light system would require a five-colour photographic system, and the corresponding LED sources, for having their five components separate (it is possible at present, but costly).We emphasize that a six-light illumination system in a hexagonal arrangement would not give any benefit over the three-light system because it is composed of two equilateral triangles.
It is important that the intensities of the the light sources are as equal as possible.This is evident from Figure 12b, which shows the result of decreasing the intensity of one light source by 40% in the three-light system.One must also ensure that the intensities at different spectral regions are similar.
In this work we use orientation of fibers in paper as the basic application and framework.This choice was made because data to be analyzed in this application are common and challenging.Therefore, if the methods developed work well in this case, they will probably work in many other (similar) applications such as, e.g., determination of the orientation of fibers or nanofibrils in reinforced composites.Furthermore, in papermaking industry, it would be advantageous to have a good orientationanalysis method for on-line measurements during the manufacturing process.The curvelet orientation algorithm used here can be implemented such that it only takes a few milliseconds to carry out the analysis.This would allow on-line analysis with an about 20 μm resolution in a paper machine with a web speed of about 20 m/s.

Fig. 1 .
Fig. 1.Setup of the experiment for recording reflection images.

Fig. 2 .
Fig. 2. (a) An original image of a paper sample.This image contains light reflected from the sample and effects of uneven illumination, (b) the estimated distribution of illumination and (c) the image resulting from removal of the estimated unevennes in the illumination.

Fig. 3 .
Fig. 3. (a) A simulated fiber network, (b) the result of a structure tensor -based orientation analysis of the network of figure (a), (c) the result of a curvelet-based orientation analysis of this network, and (d) the result of an FFT-based orientation analysis of this network.The known orientation distribution is shown by the dashed line for comparison.

Fig. 5 .
Fig. 5. Orientation distributions for four different directions of illumination by a single light source (alternating places of the light source are shown by arrows) of the same paper sample.

Fig. 6 .
Fig. 6.Orientation distributions with three light sources placed at regular angular intervals for different inclination angles. •

Fig. 7 .
Fig. 7.A schematic layout of the illumination systems composed of three (a, c) and five (b, d) regularly spaced light sources.

Fig. 9 .Fig. 10 .
Fig. 9. Fiber-orientation distributions in a paper sample (top and bottom surface) for different orientations of the light source.

Fig. 11 .
Fig. 11.The orientation distribution diagrams for individual lights in a three-light system (panels a-c) and the average diagram (panel d).Images were taken with a two-degree difference in the direction of the individual light source, so there were 90 subsets of the three-light illumination system (corresponding to its rotation by 0 to 178 • with two-degree increments, with respect to the sample).Normalized orientation intensities are described by gray-scale values defined by the gray-scale bar.

Fig. 12 .
Fig. 12.(a) The average orientation distribution of the three-light illumination system shown in Figure 11d in a threedimensional representation in which the orientation intensity appears as the third dimension (b) the three-dimensional orientation diagram (average) for a three-light system in which the intensity of one light source was decreased by 40%.The grey scale bar shows how the normalized orientation intensities are related to the grey scales.

Fig. 13 .
Fig. 13.The three-dimensional average orientation diagram for the five-light illumination system.The grey scale bar (see Fig. 12) shows how the normalized orientation intensities are related to the grey scales.

Table 1 .
Statistical characteristics for the top and bottom sides of the paper sample.