Image Anal Stereol 2014;33:147-155 doi:10.5566/ias.v33.p147-155 Original Research Paper INCLUSION RATIO BASED ESTIMATOR FOR THE MEAN LENGTH OF THE BOOLEAN LINE SEGMENT MODEL WITH AN APPLICATION TO NANOCRYSTALLINE CELLULOSE MIKKO NIILO-RÄMÄC,1, SALME KÄRKKÄINEN1, DARIO GASBARRA2 AND TIMO LAPPALAINEN3 1Department of Mathematics and Statistics, P.O. Box 35 (MaD), FI-40014 University of Jyväskylä, Finland; 2Department of Mathematics and Statistics, P.O. Box 68, FI-00014 University of Helsinki; 3VTT Technical Research Center of Finland, Koivurannantie 1, P.O. Box 1603, FI-40101 Jyväskylä, Finland e-mail: mikko.niilo-rama@jyu.., salme.m.karkkainen@jyu.., gasbarra@mappi.helsinki.., timo.lappalainen@vtt.. (Received October 25, 2013; revised May 10, 2014; accepted May 26, 2014) ABSTRACT A novel estimator for estimating the mean length of .bres is proposed for censored data observed in square shaped windows. Instead of observing the .bre lengths, we observe the ratio between the intensity estimates of minus-sampling and plus-sampling. It is well-known that both intensity estimators are biased. In the current work, we derive the ratio of these biases as a function of the mean length assuming a Boolean line segment model with exponentially distributed lengths and uniformly distributed directions. Having the observed ratio of the intensity estimators, the inverse of the derived function is suggested as a new estimator for the mean length. For this estimator, an approximation of its variance is derived. The accuracies of the approximations are evaluated by means of simulation experiments. The novel method is compared to other methods and applied to real-world industrial data from nanocellulose crystalline. Keywords: Boolean model, exponential length distribution, line segments, mean length, minus-sampling, nanocellulose crystalline, plus-sampling, ratio of estimates, variance. INTRODUCTION Fibrous structures are common in natural objects such as muscle .bres and wood .bres. Currently, increasing research efforts have been directed to isolation, production and characterization of novel nanocelluloses, being .brous structures in nanoscale. Nanocellulose can be used in food, pharmaceutical, and medical industries, which explains the importance of those new materials. Nanocelluloses may be classi.ed in three main subcategories: micro.brillated cellulose (MFC), bacterial nanocellulose (BNC) and nanocrystalline cellulose (NCC) (Klemm et al., 2011), of which the latest one is of our interest (Fig. 1). Rod-shaped NCC, also known as whiskers, is prepared from natural cellulose by acid hydrolysis. The morphology and dimensions of the whiskers depend on the native cellulose source, hydrolysis time and temperature. The analysis of particle size distribution of nanocellulose is needed mainly for two purposes: to compare and learn about different isolation/production mechanisms and certain applications may require more speci.c information about size distribution of nanocellulose. Using high-resolution microscopy techniques such as atomic force microscopy (AFM) (Pöhler et al., 2010) together with image processing techniques (Kärkkäinen et al., 2012), the information on the structure of nanocellulose can be obtained (Fig. 1). In the current work, our objective is to introduce an estimation method for the mean length of whiskers observed in nanoscale. Fig. 1. An image of rod-shaped NCC particles with identi.ed and coloured whiskers. The size of the image is 5 µm × 5 µm. In practice, the spatial system of .bre-like objects is observed through a bounded observation window (Fig. 1). When estimating the individual model parameters such as the intensity, the mean number of objects per unit area or the .bre length distribution, their estimators are often degraded by edge effects: censoring effect and spatial sampling bias (Baddeley, 1999). Censoring happens if the lengths of .bres can only be observed inside the observation window. Spatial sampling bias results from an “unfair” sampling, i.e., the longer .bres are sampled more probably than the shorter ones. Examples of unfair sampling rules are plus-sampling and minus-sampling. In plus-sampling, the .bres hitting the observation window are sampled, whereas in minus-sampling only the .bres lying completely in the observation window are sampled (Miles, 1974). The biases of plus-sampling and minus-sampling can be tackled by weighting the observations, which results in Horvitz-Thompson type estimators (Miles, 1974; Baddeley, 1999). The use of (weighted) minus-sampling requires that the observation window is large enough when compared to individual .bre length, i.e., every segment must .t in the observation window (Baddeley, 1999, p. 50). Alternatively, unbiased sampling rules such as the associated point rule can be used (Miles, 1978). Then, the .bres having their associated point (for example the northern end) in the window are sampled. Note that both the weighted plus-sampling and the associated point rule assume that the parts of .bres lying outside of the observation window can be recorded. The idea of the present paper is to introduce an alternative method for estimating the mean length of the .bres by using the ratio of two biased intensity estimators, based on plus-and minus-sampling. There is no need for the measurement of lengths, only the ratio of these two intensity estimators is needed. Therefore this method can be useful especially in such cases, where we cannot see outside the sampling window, and the data contains censored lengths. When using minus-sampling, we, however, require that the observation window is large enough. The novel estimator of the mean length is a function of the ratio. In our work, the function is determined for the Boolean model (Matheron, 1972; 1975) of line segments. We also need to make assumptions about the line segments: the direction distribution is assumed to be uniform and the length is modelled by the exponential distribution, which is a common choice in industrial processes. When compared to our previous paper (Niilo-Rämä and Kärkkäinen, 2011), we further determine the approximate analytical variance of the estimator. The accuracy of the method is evaluated analytically and with simulation experiments, where the intensity of the Boolean model and the mean length of .bres are varied. Further, the method is compared to other methods and applied for estimating the mean length of nanocrystalline cellulose (Fig. 1). DATA DESCRIPTION AND IMAGE PROCESSING The nanocrystalline cellulose was prepared from ground Whatman 541 ashless .lter paper (Kontturi et al., 2007). Sample preparation for atomic force microscopy (AFM) was conducted by spin-coating (Kontturi et al., 2007; Ahola et al., 2008) to achieve a uniformly scattered .bril layer. In spin coating, solid .lms are prepared from a dissolved or dispersed substance by removing the solvent with high-speed spinning. The substrate used in our study was a commercial silicon oxide wafer (pieces with size ~ 1 cm2), which was cleaned by rinsing with MilliQ-water and placed into UV-ozonator. Firstly, cellulose nanocrystal suspension was diluted to the concentration of approximately 0.005 m-%. Secondly, suspension was spin-coated with a speed of 4000 rpm for 30 s. Finally, the substrate was carefully rinsed to remove unattached .brils and oven-dried (80.C, 10 min). The samples were stored in desiccator before imaging. Images were taken with Nanoscope IIIa multimode scanning probe microscope from Digital Instruments Inc at Aalto University. The images were scanned in tapping mode in air. The size of the images was 5 µm × 5 µm, see an example in Fig. 1. For analyzing the structure of individual whiskers in the obtained image, a set of image processing techniques needs to be performed in order to form a pixel sequence for each whisker. In the image processing the image was .rst .ltered using bandpass .lter (Gonzalez and Woods, 2002) and median .lter (Nisslert et al., 2007) in order to reduce noise. Second, the image was binarized using isodata thresholding (Ridler and Calvard, 1978). Then, the image was dilated to remove small particles and .nally skeletonized (Gonzalez and Woods, 2002). The following steps of the image processing are described in more detail in Kärkkäinen et al. (2012). Starting from the skeletonized image, each pixel point of the image was .rst categorized into one of four classes: background point, end point, branch-intersection point and normal skeleton point with a given rule. Second, some detached but close intersection areas were merged in order to form real physical intersection areas. Short connecting parts of whiskers and also intersection areas without any connected part of Image Anal Stereol 2014;33:147-155 whiskers were removed. Then, in a certain intersection area a weight for each pair of whiskers was calculated, having a low value in the case of similar curvatures and different directions of the whiskers. The pair with the lowest weight was connected to form a single whisker in favour of straighter and longer whiskers. The whiskers with random colours are illustrated in Fig. 1. THE BOOLEAN MODEL Let us consider a marked point process . = {xi,Ki}, where the locations {xi} follow a point process in R2 and the marks {Ki} are random compact sets in R2 . A germ-grain model with germs {xi} and grains {Ki} (Hanisch, 1981) is the union . =(xi + Ki) . i The germ-grain model is assumed to be stationary, in which case we can de.ne a “typical” grain K0, which is a random closed set. Its distribution Q is the mark distribution of . on the space Kof compact sets in R2 (Stoyan et al., 1995, p. 216; Chiu et al., 2013). In this work, we assume that {xi} form a stationary Poisson point process in R2 with intensity . and Ki is a line segment with random length Li and angle Ai from common distributions fL(l) and fA(.), respectively. The line segments are independent of each other and, further, independent of the points. This type of model is called a Boolean model (Matheron, 1972). In addition, we assume that the line segments are separable. THE RATIO ESTIMATOR Next we introduce a novel method for estimating the mean length of line segments observed in a square shaped observation window W, which is a convex and compact set in R2 . The idea is to use the ratio of two biased intensity estimators in the estimation of the mean length of the segments (Niilo-Rämä and Kärkkäinen, 2011). MINUS-AND PLUS-SAMPLING Let us .rst recall the basic results of minus-and plus-sampling when estimating the intensity of a stationary germ-grain model with separable segments {Xi = xi + Ki}. Without loss of generality, we can assume that W is a unit square. When using plus-sampling, the estimator of the intensity . is the number of .bres hitting W, i.e., #{i : Xi . W /The expectation of the estimator = 0}. is obtained by the Campbell-Mecke theorem (cf. Baddeley, 1999), E[#{i : Xi .W / = 0}] = E.1(Xi .W = 0/)i = .E0 1((K0 + x) .W / = 0)dx R2 2 = .E0|W . K¡0| , (1) where E0 is the expectation with respect to Q and W › W . K¡0 = {x . R2: (K0 + x) .W = 0/} is the dilation of W , and |·| denotes the surface area. Consequently, using minus-sampling, the expected number of .bres included in W is given by E[#{i : Xi . W }]= .E0 [|W 8 K0|], (2) where W › W 8 K0 = {x . R2: (K0 + x) . W } is the erosion of W . Instead of having . on the right-hand sides of Eqs. 1 and , we are dealing with sampling biases E0 |W . K¡0| 2 22 and E0 [|W 8 K0|]. Our target is to calculate both of them and use their ratio in the estimation of the mean length of the segments. CONSTRUCTION OF THE ESTIMATOR Let us assume a Boolean model with line segments having a random length L ~ Exp(1/.) with E[L]= . and a random direction A ~ U[0,2.) with the horizontal axis. In addition, the length and the direction are assumed to be independent of each other. Using our assumptions, we derive P(Xi . W |Xi .W = 0/), the conditional probability that a line segment Xi is included in W given it hits W , i.e., a line segment sampled using plus-sampling would also be sampled in minus-sampling. Finally, we relate this conditional probability to the mean length of the segments with a given length density. Mathematically, the conditional inclusion probability of a line segment conditional on hitting the window equals the ratio of the expected sample sizes of minus-and plus-samplings: E[#{i : Xi . W }] P(Xi . W |Xi .W /. (3) = 0)= E[#{i : Xi .W / = 0}] Using Eqs. 1 and 2, we obtain from Eq. 3 E0[|W 8 K0|] P(Xi . W |Xi .W /. (4) = 0)= E0[|W . K¡0|] The left-hand side of Eq. 4 can be estimated from the data with the ratio of the intensity estimates of minus-and plus-sampling. In order to use that, we need to calculate the right side of Eq. 4 and its relation to the length distribution. Recall that W is a unit square. With a .xed length l and a direction ., the area for the erosion can be written in the form |W 8 K0| = 1(l|sin.|. 1)1(l|cos.|. 1) × (1 - l|cos .|)(1 - l|sin.|) and for the dilation |W . K¡0| = 1 + l|cos.| + l|sin.| . Next, let us assume that the direction and the length are random, with density functions fA(.) and fL(l), respectively. In that case, W 8 K0 and W . K¡0 are random compact sets. Then, the expected area for the erosion can be given by E0[|W 8 K0|]= |W 8 K0|dQ(K0) K . 2. = fL(l) fA(.)(1 - l|cos.|)+ 00 × (1 - l|sin.|)+ d. dl . (5) The solution of the integral in Eq. 5 is not available in closed form but can be approximated by . 2. E0[|W 8 K0|] . fL(l) fA(.)(1 - l|cos.|) 00 × (1 - l|sin.|)d. dl . (6) This approximation is quite accurate when l is small enough (Fig. 2 and the simulation examples later). Applying some algebra and properties of trigonometric functions, the right-hand side of Eq. 6 yields Recall that L ~ Exp(1/.) with E[L]= ., in which case E[L2]= 2.2. Combining Eqs. 6–8, we can write Eq. 4 as a function of .: 2.2 1 - 4. + .. P(Xi . W |Xi .W = 0/) . =: p(.) . (9) 4. 1 + . Note that we are now dealing with a unit square, i.e., |W | = 1. The estimation method is, however, not restricted to unit squares, as we can adjust ., which is actually the ratio of the mean length of the .bres and the length of the window side. The ratio of Eq. 9 can be used for all sizes of squared observation windows if we can assume the restriction . . 1. Then, the function p(.) is continuous and strictly monotonic. In practice, when having a realization of a line segment process, we get an estimate p^ for the inclusion probability p(.) by counting the ratio of number of segments lying completely inside W and number of segments intersecting W (assuming that the denominator is greater than zero). Then the estimator of . is obtained by using the inverse function p-1, i.e., .^ = p-1(p^) . (10) .. . fL(l)dl - 4 l fL(l)dl + 1 l2 fL(l)dl 0 . 0 . 0 = 1 - 4 E[L]+ 1 E[L2] . (7) .. For the dilation we are able to calculate the accurate expectation, which is E0[|W . K¡0|]= |W . K¡0|dQ(K0) K . 2. = fL(l) fA(.)(1 + l|cos.| + l|sin.|)d. dl 00 4 = 1 + E[L] . (8) . Fig. 2. The graph of the function p(.) (solid line) together with a scatter plot of simulated inclusion ratios using different values between 0.1 and 1.0 for . and the intensity . = 30 for line segments. APPROXIMATE VARIANCE OF THE ESTIMATOR Next we are going to derive the theoretical variance of the inclusion ratio estimator. The assumptions are the same as before, i.e., a Boolean model with intensity . , exponentially distributed .bre lengths and uniformly distributed directions. Image Anal Stereol 2014;33:147-155 When estimating the mean length of the .bres variance of the mean length estimator .^ = p-1(p^) we .rst have to estimate the theoretical inclusion probability p(.) from our sample. To be more precise, Var .^ 2 2 1 . Var[p^] the estimator is pI(.) 2 p^ = N- N+ 16.2 + 8.. + .2 1(N+ > 0)+ p1(N+ = 0) , = p(1 - p) 8.2 + 4.. - 8. 1(N+ > 0) where N- is the number of .bres lying completely × E . (13) inside our window (sample size using minus- N+ sampling) and N+ is the number of .bres hitting The expectation factor in Eqs. 12 and 13 is given by the window (sample size plus-sampling) and p is the probability shown in Eq. 3. Then . 1(N+ > 0) -.+ .+ k E = e. N- N+ k=1 k!k 1(N+ > 0) Var[p^]= Var 2 N+ = .+E (1 + N+)-2+ p2Var[1(N+= 0)] = 2F2(1,1;2,2;.+)exp(-.+).+ , N- + 2pCov 1(N+ > 0),1(N+= 0) . where N+ (11) kF (a1,...,ak;b1,...,bF; .+) Using the law of total variance, the .rst term of the ..(a1 + n)....(ak + n) . sum on the right-hand side of Eq. 11 can be written as := n=0 .(a1)....(ak) .(b1)....(bF) .n + N- × 1(N+ > 0) E Var N+ .(b1 + n)....(bF + n) n! N+ N+ . N- is the generalized hypergeometric function. In the numerical computations we have used the expansion 1(N+ > 0) + Var E N+ 2 E N-1 + 2 1(N+ > 0) m (k - 1)! Since N- N+ ~ Bin(N+, p), this equals + O(.-(m+1) + ) . . k + k=1 N+ p(1 - p) pN+ E 1(N+ > 0)+ Var 1(N+ > 0) N2 N+ + as .+ › ., given in Jones and Zhigljavsky (2004). 1(N+ > 0) = p(1 - p)E N+ SIMULATION EXPERIMENTS 2 - E[1(N+ > 0)] 2 2 E 1(N+ > 0)2 + p . The estimator .^ = p-1(p^) and the theoretical formula for its variance, Eq. 13, are based on 2 Now N+ ~ Poisson(.+), where .+= .E0 |W . K¡0|the expected sample size of plus-sampling (Eqs. 1, 8). , approximations. Therefore, we examined the accuracy of the estimator and the analytical approximation of its variance by simulation experiments using R- Since E[1(N+ > 0)] = P(N+ > 0)= 1 -e-.+ , the .rst term of Eq. 11 reduces to software (R Core Team, 2013). With the same model assumptions as before, we simulated 10 000 realizations of a Boolean line segment model in a 1(N+ > 0) 2 p(1 - p)E + pe-.+(1 - e-.+) . N+ unit square using four different intensities (. = 20, . = 30, . = 50 and . = 100) and three different mean lengths (. = 0.1, . = 0.2 and . = 0.5). From these Applying similar computations, the second term of 2 e-.+(1 - e-.+) and the third term of Eq. 11 equals p Eq. 11 equals -2p 2realizations we computed the empirical means and e-.+(1 - e-.+). Hence standard errors and also the theoretical approximations for the standard errors of the estimator using the square 1(N+ > 0) Var[p^]= p(1 - p)E . (12) N+ root of Eq. 13 (Tables 1–4). According to Tables 1–4, the inclusion ratio based As .+ › . and p = p(.), by using the delta estimator seems to work quite accurately. As was method (Davison, 2003), we obtain the asymptotic expected, the bias is small, as is the variance, when the intensity is large and the line segments are short. This means that the observation window should be chosen to be large enough in order to obtain an unbiased estimator for . with a tolerable standard error. Table 1. Simulation results using intensity . = 20. . = 20 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) theoretical S.E. (.^) 0.1018 0.0487 0.0471 0.2039 0.0726 0.0692 0.5233 0.1868 0.1299 Table 2. Simulation results using intensity . = 30. . = 30 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) theoretical S.E. (.^) 0.1012 0.0389 0.0381 0.2012 0.0571 0.0560 0.5009 0.1366 0.1055 Table 3. Simulation results using intensity . = 50. . = 50 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) theoretical S.E. (.^) 0.1002 0.0294 0.0294 0.1994 0.0560 0.0432 0.4880 0.1055 0.0807 Table 4. Simulation results using intensity . = 100. . = 100 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) theoretical S.E. (.^) 0.0990 0.0207 0.0206 0.1984 0.0308 0.0302 0.4826 0.0563 0.0571 REAL DATA The novel method was applied to the processed nanocrystalline cellulose image (Fig. 1). From the image 794 .bres were detected, 757 of them lying completely inside the window, i.e., not touching the edge. As seen in Fig. 3, the assumption about the exponential length distribution seems to hold quite well. Only the portion of very short .bres seems to be too small, this might be due to the fact that some of the shortest existing .bres (< 1 pixel in the image) were not detected. Fig. 3. Histogram of the measured .bre lengths in the nanocrystalline cellulose image together with the graph of the exponential distribution (solid line) with parameter . = 9.57. The estimated inclusion probability was p^ = 757/794 = 0.9534. Further, using the inclusion ratio based estimator, the obtained estimate for the mean length was .^ = 0.0188. The approximated standard error was 0.0032. Scaled to the image size, the estimated mean length was about 9.57 pixels (0.09 µm) and the standard error 1.67 pixels (0.02 µm). For comparison purposes, the mean length of .bres identi.ed from the picture was 0.0210, that is about 10.70 pixels (0.11 µm). The standard error of the sample mean estimator was 0.0006, which is about 0.3 pixels (< 0.01 µm). Next we made simulations using those parameter values estimated from the real data, i.e., . = 773 (estimated using the associated point rule) and . = 0.0210 (see an example in Fig. 4). Using the novel method for 10 000 realizations, the following results were obtained: sample mean(.^)= 0.0207 and empirical S.E.(.^)= 0.0032, which agrees with the analytically computed standard error. Image Anal Stereol 2014;33:147-155 Fig. 4. Simulated realization of a line segment process with intensity . = 773 and mean length . = 0.0210. COMPARISON WITH OTHER METHODS We compared our novel method with two approaches. First, we used a stereological approach, where the .bre system is intersected with sampling lines. In the isotropic case, we have a stereological formula (Mecke and Stoyan, 1980) . LA = PL , 2 where LA is the expected total .bre length per unit area (in our case LA = ..), and PL is the expected number of intersection points with .bres per unit length of a sampling line in any direction. Now the estimator for PL is, for example, the number of intersection points between the .bres and the boundary of the sampling window W , divided by the length of the boundary .W . Then the unbiased estimator for the total length in W is . ^L = 2 ^PL|W | . Further, an estimator for the mean .bre length suggested by the anonymous referee is .^ = 2^L Q , (14) where Q is the number of all endpoints of .bres in W . When using the associated point rule (Miles, 1978), the expected number of northern points of segments lying inside the observation window is . |W |. Consequently, the expected value for the number of all end points Q is 2.|W |. It can be shown that the estimator of Eq. 14 is ratio-unbiased. Eq. 14 is also applied when the total .bre length is measured from known .bre lengths in the image of real data and in simulation experiments. For the real data, with the stereological approach, the obtained estimate for the mean length was 10.01 pixels (0.10 µm), whilst with the total length based estimator, the estimate was 11.55 pixels (0.11 µm). SIMULATION EXPERIMENTS For comparison purposes, we made some simulation experiments and estimated the mean length using the two alternative estimators. With chosen parameter values, the simulation results show that the mean length estimator based on the total length measure (Tables 6 and 8) has smaller standard error than the ratio estimator (Tables 2 and 4). In Tables 5 and 7, where stereological approach is used, we obtained approximately the same accuracy as with the ratio estimator (Tables 2 and 4). These methods are both based on indirect measurement of the total length. Table 5. Simulation results using intensity . = 30 and stereological estimator . = 30 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) 0.0995 0.0385 0.2000 0.0548 0.5006 0.0956 Table 6. Simulation results using intensity . = 30 and total length based estimator . = 30 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) 0.0999 0.2000 0.4984 empirical S.E. (.^) 0.0179 0.0353 0.0816 Table 7. Simulation results using intensity . = 100 and stereological estimator . = 100 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) 0.0992 0.0204 0.1989 0.0425 0.4990 0.0739 . = 100 . = 0.1 . = 0.2 . = 0.5 sample mean (.^) empirical S.E. (.^) 0.0998 0.0097 0.1996 0.0265 0.4980 0.0632 Furthermore, we made 10000 simulations using the estimated parameter values of the nanocellulose data, i.e., intensity . = 773 and mean length . = 0.0210. For the total length based estimator we obtained sample mean 0.0209 with sample S.E. 0.0007. For the stereological estimator, the sample mean was 0.0207 and the sample S.E. 0.0033, which are almost exactly the same as with the ratio estimator. CONCLUSION A novel estimation method for the mean length of line segments was proposed together with accuracy studies. Moreover, the method was applied to nanocrystalline cellulose being currently material of great interest in industry. Classically, the mean length of line segments observed in an observation window has been estimated using such methods as weighted minus-or plus-sampling. Assuming a Boolean model, we introduced an alternative method based on the ratio of the random sample sizes of plus-and minus-samplings. The novel estimator is a function of the ratio. We determined the function for the Boolean model of line segments with an exponential length distribution and a uniform direction distribution and further the approximate variance of the estimator. The method is approximate as well as the obtained theoretical variance of the estimator. Therefore, the accuracy of the inclusion ratio based estimator was evaluated both theoretically and by simulations, which gave promising results. The method was also compared to the estimators based on the stereological approach and the total length measurement. The novel method was comparable with the stereological estimator which is based on indirect measuments as our method. An advantage of our method is that we have an approximate formula for the variance. The methods based on the exact length measurements had smaller variance in simulation experiments with chosen parameter values and in real data, if the formula is available. Besides the variance, other advantages of the novel method may be the following: it is simple and requires possibly less work since there is no need to measure the lengths of individual segments. This is especially an advantage, when the exact identi.cation of segments is challenging, but the ratio may be easily available. As a disadvantage, it should be noted that the novel method is based on the minus-sampling and therefore it is required that the observation window is large enough when compared to the individual .bre length. In that case, the censoring effect is handled automatically. The expansion of the method for other length and direction distribution models may be one of the challenges in the future. In addition, the method could possibly be generalized into objects or sampling windows with different shapes as considered in this work. ACKNOWLEDGEMENTS The authors would like to thank Prof. Antti Penttinen for fruitful discussion, the referees for professional remarks, Eero Kontturi, Dr., for image data, and Arttu Miettinen, M.Sc., for image processing. The research was partly funded by STATCORE (S. Kärkkäinen) and COMAS (M. Niilo-Rämä). REFERENCES Ahola S, Salmi J, Johansson, LS, Laine J, Österberg M (2008). Model .lms from native cellulose nano.brils. Preparation, swelling and surface interactions. Biomacromolecules 9(4):1273–82. Baddeley AJ (1999). Spatial sampling and censoring. In: Barndorf-Nielsen OE, Kendall WS, van Lieshout MNM, Eds. Stochastic geometry: likelihood and computation. London: Chapman and Hall. pp 1–78. Chiu SN, Stoyan D, Kendall WS, Mecke J (2013). Stochastic geometry and its applications, 3rd Ed. Chichester: Wiley. Davison AC (2003). Statistical models. Cambridge: Cambridge University Press. Gonzalez RC, Woods RE (2002). Digital image processing, 2nd Ed. Upper Saddle River: Prentice Hall. pp. 245–6, 523–5, 534–45. Hanisch KH (1981). On classes of random sets and point processes. Serdica 7:160–6. Jones CM, Zhigljavsky AA (2004). Approximating the negative moments of the Poisson distribution. Stat Probabil Lett 66:171–181. Kärkkäinen S, Nyblom J, Miettinen A, Turpeinen T, Pötschke P, Timonen J (2012). A stochastic shape and orientation model for .bres with an application to carbon nanotubes. Image Anal Stereol 31(1):17–26. Klemm D, Kramer F, Moritz S, Lindström T, Ankerfors M, Gray D, Dorris A (2011). Nanocelluloses: A new family of nature-based materials. Angew Chem Int Ed 50(24):5438–66. Kontturi E, Johansson LS, Kontturi KS, Ahonen P, Thüne P, Laine J (2007). Cellulose nanocrystal submonolayers by spin coating. Langmuir 23(19):9674–80. Matheron G (1972). Ensembles fermés aléatoires, ensembles semimarkoviens et polyedres poissoniens. Adv Appl Prob 4:508–41. Image Anal Stereol 2014;33:147-155 Matheron G (1975). Random sets and integral geometry. New York: Wiley. Mecke J, Stoyan D (1980). Formulas for stationary planar .bre processes I – general theory. Math Operationsforsch Statist Ser Statist 11:267–79. Miles RE (1974). On the elimination of edge-effects in planar sampling. In: Harding EF, Kendall DG, Eds. Stochastic geometry. London: Wiley. 228–47. Miles RE (1978). The sampling, by quadrats, of planar aggregates. J Microsc 113:257–67. Niilo-Rämä M, Kärkkäinen S (2011). Inclusion ratio based estimator for the mean length of the Boolean line segment model. In: Proc 13th Int Congress Stereol (ICS-13), Oct 19-23, Beijing, China. pp 520–3. Nisslert R, Kvanström M, Lorén N, Nydén M, Rudemo M (2007). Identi.cation of the three-dimensional gel microstructure from transmission electron micrographs. J Microsc 225:10–21. Pöhler T, Lappalainen T, Tammelin T, Eronen P, Hiekkataipale P, Vehniäinen A, Koskinen T (2010). In.uence of .brillation method on the character of nano.brillated cellulose (NFC). In: Proc 2010 TAPPI Int Conf Nanotech Forest. Sept 27-29, Espoo, Finland. pp 437–58. R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project. org/. Ridler TW, Calvard S (1978). Picture thresholding using an iterative selection method. IEEE T Syst Man Cyb 8:630–2. Stoyan D, Kendall WS, Mecke J (1995). Stochastic geometry and its applications, 2nd Ed. Chichester, Wiley.