Characterization of fiber and vessel elements in pulp suspension images

(1)

Degree Program in Information Technology

Master’s Thesis

Tatiana Kurakina

CHARACTERIZATION OF FIBER AND VESSEL ELEMENTS IN PULP SUSPENSION IMAGES

Examiners: Professor Heikki Kälviäinen D. Sc. (Tech.) Tuomas Eerola Supervisors: M. Sc. Nataliya Strokina

Professor Heikki Kälviäinen D. Sc. (Tech.) Lasse Lensu D. Sc. (Tech.) Tuomas Eerola

(2)

ABSTRACT

Lappeenranta University of Technology Faculty of Technology Management

Degree Program in Information Technology Tatiana Kurakina

Characterization of fiber and vessel elements in pulp suspension images Master’s Thesis

2012

51 pages, 22 figures, 3 tables, and 1 appendix.

Examiners: Professor Heikki Kälviäinen D. Sc. (Tech.) Tuomas Eerola

Keywords: image processing, machine vision, image segmentation, tensor voting, curve extraction

The thesis is related to the topic of image-based characterization of fibers in pulp suspension during the papermaking process. Papermaking industry is focusing on process control optimization and automatization, which makes it possible to manufacture high- quality products in a resource-efficient way. Being a part of the process control, pulp suspension analysis allows to predict and modify properties of the end product. This work is a part of the tree species identification task and focuses on analysis of fiber parameters in the pulp suspension at the wet stage of paper production.

The existing machine vision methods for pulp characterization were investigated, and a method exploiting direction sensitive filtering, non-maximum suppression, hysteresis thresholding, tensor voting, and curve extraction from tensor maps was developed. Ap- plication of the method to the microscopic grayscale pulp images made it possible to detect curves corresponding to fibers in the pulp image and to compute their morphological characteristics. Performance of the method was evaluated based on the manually produced ground truth data. An accuracy of fiber characteristics estimation, including length, width, and curvature, for the acacia pulp images was found to be 84, 85, and 60%

correspondingly.

(3)

PREFACE

I wish to thank my practical supervisor Nataliya Strokina for her patience, discussions and answers to unlimited questions. Without her ideas and support this thesis would have never happened. I am also grateful to Heikki Kälviäinen, Lasse Lensu and Tuomas Eerola for their guidance.

I wish to thank my family, especially Mom and Dad, who were always eager to help and support me while writing the thesis and during the rest of the life. I believe, this thesis is our joint achievement.

I am grateful to my friends, who never ceased to believe in me and kindly motivated to work persistently.

Lappeenranta, November 7th, 2012

Tatiana Kurakina

(4)

ABBREVIATIONS AND SYMBOLS

FDR false discovery rate FP false positives FN false negatives GT ground truth

LoG Laplacian of Gaussian

MVPR Machine Vision and Pattern Recognition PPV positive predictive value or precision TP true positives

TPR true positive rate

CIacc accuracy of experimental curl index evaluation CI_exp experimental fiber curl index evaluation CI_GT GT fiber curl index estimation

F1 F-measure of a test accuracy

L_acc accuracy of experimental fiber length evaluation L_exp experimental fiber length evaluation

L_GT GT fiber length evaluation

N_exp number of experimentally extracted fibers NGT number of GT fibers

T time required for processing of one image T₁ time required for edge detection

T₂ time required for tensor voting T₃ time required for curve extraction

W_acc accuracy of experimental fiber width evaluation W_exp experimental fiber width evaluation

W_GT GT fiber width evaluation

(7)

1 INTRODUCTION

1.1 Background

Wood is a natural material and therefore it cannot be homogeneous. Its properties vary not only within wood types (hardwood or softwood), but also within one family of wood.

The variations of wood properties cause inhomogeneous properties of pulp, which is fed into the paper machine. The characteristics of pulp used in papermaking determine the properties of the produced paper. Since a suspension typically includes pulp of different tree species, analyzing the pulp properties at the early stages allows the properties of the produced paper to be predicted [1]. In papermaking industry, the analysis has been traditionally made in laboratory. Firstly, the samples are taken and then the tests are performed. Papermaking industry nowadays is interested in transferring of the laboratory measurements to the in-line measurements process. Therefore, pulp analysis methods suitable for industrial usage are needed.

In the Machine Vision and Pattern Recognition Laboratory (MVPR) of Lappeenranta Uni- versity of Technology [2] image analysis methods for machine vision applications for paper and printing industry are being developed. The PulpVision project [3] focuses on the development of image-based measurement and characterization methods related to the quality of pulp as a raw material for papermaking. The presented work is held within the PulpVision project and focuses on image-based characterization of pulp suspension.

According to the tree species identification task, pulp can be classified by the properties of its components. Characteristics of the pulp can be computed from the images of pulp suspension. Pulp suspension is inhomogeneous and consists of different elements. The main components of the pulp suspension are fibers and vessel cells [1]. Fibers in the suspension interact with each other, forming fiber bundles. In the production process some of the fibers are mechanically damaged, causing appearance of small segments of fibers, called fines. Fines enhance optical properties of paper, however, impairing the strength properties. Abrupt changes in fiber curvature (kinks) affect the formation of paper and its strength properties. Vessel segments are dead and hollow objects, perforated at the ends. The size, character, and number of vessel cells differ between tree species and makes it possible to classify pulp species. Pulp species can be also classified based on fiber length to width distribution [4].

(8)

1.2 Objectives and Restrictions

The master thesis is related to the task of pulp suspension characterization and classification. The scope of this thesis is to detect the fibers in the suspension and to compute their characteristics.

The research tasks can be formulated as follows:

1. to investigate existing systems for pulp suspension characterization and develop a method for fiber segmentation and characterization,

2. to apply the developed method to the real pulp images and compute the fiber parameters,

3. to create the ground truth data and perform result evaluation.

The real test data include microscopic pulp suspension images of acacia made with 2.5 magnification, captured in laboratory conditions. In the presented work the data is produced by imaging suspension flows and is therefore closer to in-line measurements than traditional laboratory measurements.

1.3 Structure of the Thesis

In Section 2, a literature review of existing methods is presented. In Section 3, the developed method for fiber detection and characterization is described. In Section 4, the test data is introduced, experiments with the developed method are presented, evaluation of the results is performed, and the results are discussed. Finally, in Section 5, the conclu- sions are drawn.

(9)

2 PAPER PULP CHARACTERIZATION

2.1 Pulp and papermaking

Papermaking is a complicated technological process that utilizes a large number of resources. It starts with removing bark from logs and cutting wood into smaller pieces.

That is called chipping. Wood chips are mechanically or chemically disintegrated as a result of which wood pulp is obtained. Cellulose fibers in the pulp are separated from each other. After that, at the wet stage, the cellulose fibers are diluted with water, chem- icals are added, and pulp suspension with necessary properties is prepared. Finally, this suspension is refined, cleaned from dirt particles, and fed into the papermaking machine, which produces paper [1]. Scheme of papermaking process is presented in Figure 1.

Figure 1.Papermaking process [5].

Raw materials, intermediates, and end products are analyzed at all stages of industrial production in order to control the papermaking process [6]. The thesis focuses on characterization of pulp suspension during the wet stage of paper production.

It is important to measure properties of pulp at the wet stage so that chemists could im- prove the technological process and if the process is controlled on-line, analysis gives faster response, therefore, paper is produced more economically by environmentally- friendly processes with optimization of resource consumption [6].

(10)

The characteristics of fibers have a significant effect on qualities of the pulp and the end product. For example, fiber length and coarseness affect the flocculation and the mo- bility of fibers [7]. Thus, measurements of pulp suspension properties are an important component of the efficient process and product quality control [6].

2.2 Fiber analyzers

Solutions for the pulp characterization problem using the image-processing methods have been developed since the 90s [8]. According to the review of Hirn and Bauer [7], there are several well-known commercial analyzers: Pulp Fibre Analyzer (Galai CIS-100) [9], FiberLab [10], FiberMaster [10], Fiber Quality Analyzer(FQA) [11], MorFi [12], FS200 [10].

All of them take pulp samples and analyze them in laboratory. FQA, for instance, incor- porates the cytometric flow cell that orients fibers for precise measurements [11]. This approach is not suitable for in-line measurements of the papermaking process, where pulp is analyzed without sample withdrawal from the technological process. Development of new methods for in-line measurements is motivated by an absence of industrial analyzers that could be embedded into the papermaking equipment.

2.3 Analysis of the pulp suspension images

Four fundamental steps of a machine vision system [13] in terms of pulp suspension image analysis are formulated as follows:

1. Image acquisition.

Since acquisition of the pulp images is complicated by insufficient light, high tem- perature, and continuous movement of the pulp, the images are captured with spe- cialized equipment [14]. The stroboscope is used to produce light pulses. While flash of light is discharged, a microscopic image of the pulp is taken. The amount of light has to be optimally balanced so that all objects could be clearly visible. Fiber analyzers, using flow cells, are able to investigate the pulp from several sides [8].

2. Image preprocessing.

Quality of the images, acquired in industrial conditions, is typically poor. Image enhancement, including, for example, illumination correction and noise removal, is

(11)

applied to make details of the image more distinguishable [13]. Redundant objects in the image can be segmented at the stage of preprocessing as well.

3. Fiber segmentation.

The term “segmentation“ usually denotes separation of objects from the background [13].

In the present thesis segmentation includes segmentation of fibers from the background and separation of fibers from each other (localization). Presence of fiber bundles in the pulp image complicates the task of fiber localization and parametrization. In Sections 2.4 and 2.5 fiber segmentation methods are discussed.

4. Fiber characterization.

The fibers are represented as parameterized curves. Morphological properties of fibers are computed, including length, width, and curvature [8]. In Section 2.6 fiber morphology is considered.

2.4 Fiber segmentation

At the segmentation stage, objects are separated from the image background. Classical approach for image segmentation is thresholding. The basic idea of the method is to determine the threshold value and classify pixels of the input image according to their intensities: if intensity of the pixel is greater than the threshold value, the pixel is classified as an ’object’ pixel, otherwise it is a background pixel. According to the survey of Sezgin and Sankur [15], histogram shape-based (Weszka and Rosenfeld [16]), clustering-based (Otsu [17]), entropy-based (Kapur [18]) and other methods have been developed in order to determine the threshold value. The thresholding methods are suitable for segmentation tasks, where objects can be separated by the intensity value. One could take a grayscale pulp image as an input and find distribution of intensity values by calculating histogram of the image. In ideal case peaks of the histogram correspond to background and object pixels. The intensity value between this peaks can be selected as the threshold value. All pixels of the image with intensity values smaller than this threshold value will be set to zero, all the rest - to one. As a result of this transformation input image is represented as a binary one, where objects and background pixels have different intensity values.

When the foreground pixels have been separated from the background and objects in the thresholded images are presented as areas of nonzero pixels, geometrical and topological characteristics of these thresholded objects can be represented by the skeletons of fiber image. Skeleton of a set of pixels X is defined as a collection of the centers of all the

(12)

maximal elements contained in it [19]. Medial axis is a special case of skeleton [20], that represents an object as an infinite union of balls [21]. It is important to represent each fiber by a thinner object for further fiber intersection processing and fiber analysis.

For this purpose distance transform [22] can be calculated by labeling each pixel with the shortest distance to the background. The medial axis can be defined as the local maximum of the distance transform labels [23]. However, a skeleton acquired using this method is disconnected and sometimes has more than one pixel in width. The Hilditch thinning algorithm, described in [24], performs multiple thinning passes on the image, so that the final skeleton is connected, centered, and stay as thin as possible.

Morphological operations are another approach for skeletonization: skeletons of the thresholded areas can be defined as a result of mathematical morphological operations with a structural element in a binary image [13].

Since the pulp images were captured in difficult light conditions, background pixels may have the same intensity as the badly illuminated fibers and some data can be lost. More- over, skeleton representation of fibers causes difficulties with separation of fiber intersection areas. Therefore, classical segmentation method based on thresholding is not applicable for the task of fiber segmentation.

2.4.1 Edge detection operators

Object edges are identified as points, where the image brightness sharply changes. Edges in the binary image are simply defined as a boundary between one and zero pixels. When the image is filtered with a predefined operator, the boundary between two regions becomes more distinguishable. For the edge detection approach different mask have been developed: the well-known method of Canny [25], Sobel, Prewitt, Roberts, Harris, Lapla- cian of Gaussian [13]. These operators are extremely sensible to noises, because they use the first or the second derivative of adjacent edges only. An edge detection method should have good anti-interference (resistance to noise) ability and a better quality of the edge testing [26]. In [27], Hui used Canny edge detector and local curvature function for pulp corner detection problem, though evaluation of the method performance on pulp images has not been published yet.

Ridge edge detection. Ridge edge detector consists of image presmoothing by a Gaus- sian kernel followed by non-maximum suppression. It is able to detect points, at which the gradient magnitude is a maximum in the gradient direction. In this way, edges can be

(13)

detected at any scale in scale-space [28]. The ridges (or a set of curves whose points are local maxima in one dimension) are obtained as a result [29].

Line detection. Line segments can be also detected by filtering an image with masks, corresponding to vertical, horizontal, 45 and -45 degree lines. Response of the mask is calculated at every point of the image. The points are associated with a line in the direction of mask. In order to detect all lines in the image in the direction defined by a mask, the filtered image is thresholded [13]. This method was investigated in the Master⁰s thesis of Laaksonen [30], where a new mask, giving accurate results on fibers, was described.

2.4.2 Wavelet transforms

The two-dimensional wavelet transform can he seen as a one-dimensional wavelet transform along thex- andy- axes. The image is decomposed using a high- and low-pass filters into four parts: vertical higher frequencies (horizontal edges), horizontal higher frequencies (vertical edges), higher frequencies in both directions, and lower frequencies [31].

Then both horizontal and vertical details are recomposed, and the recomposed images are fused again [32].

Directional wavelet transform is a form of wavelet transform that can be described by imposing the direction constraint on the wavelet function. Idea of the refined method is to calculate separable wavelet transformφ(x, y)along sets of different directionsθ_i [33]:

φ(x, y, θ_i) =φ(xcosθ_i, ysinθ_i), (1) Recognition of linear objects based on edge detection with directional wavelet transform was applied to the fiber detection problem by Quin [26]. This transform reflects the edge characteristics and texture of the fiber image better than the original wavelet. The experiments of Quin [26] suggest that this method is robust and reliable in dealing with pulp fiber images.

Curvelet transform. Introduced by Candes and Donoho [34], a curvelet transform is based on the Fourier transform, wavelet transform, and radon transform. The idea is to represent a curve as a superposition of functions of various lengths and widths, obeying the scaling lawwidth=length²[35]. This can be done by first image decomposition into subbands, i.e. separation of the object into a series of disjoint scales. Each scale is then analyzed by the means of a local ridgelet transform (a traditional radon transform [36]

(14)

combined with a wavelet transform, introduced by Candes in [37]). Curvelets are based on multiscale ridgelets combined with a spatial bandpass filtering operation to isolate different scales.

Curvelet transform is used in the edge processing of pulp fiber images for fiber edge feature extraction and in contrast enhancement for noise removal as well (example is presented in Figure 2) [38].

Figure 2. Using curvelet transform for processing the pulp fiber image [38].

The curvelet transform based edge detector outperforms the traditional Sobel operator, Canny edge detector and wavelet transform [38]. In [38], Bian and Qui presented an algorithm, which combined edge extraction method, based on the curvelet transform, and skeleton tracing method, based on the preferential orientation, for recognition of fibers and their intersections in pulp fiber images. However, in the proposed approach the edge information is combined with fiber skeleton that cannot be precisely calculated for all fiber intersections. Therefore, this is not the best solution for the fiber segmentation and localization task.

2.5 Fiber localization and parameterization

Input of fiber parametrization is a binary image with nonzero pixels corresponding to the fibers. The task is to decide, which pixel corresponds to which fiber, and to fit a curve to each fiber. In certain laboratory measurements, fibers are intentionally rotated to make the segmentation easier. This, however, is not the case with the in-line measurements. Fibers in the pulp images can intersect with each other or be located one under another. In order to process these cases, overlapping fiber edges have to be separated (possible intersections are shown in Figure 3).

(15)

Figure 3.Possible fiber intersections [26].

In order to extract the junction of fibers and separate the intersected fibers, local curvature function can be calculated. The curvature of a point (x_i, y_i) on the contour can be expressed as the difference of the angle between tangent of two adjacent points andx-axis as follows [27]:

f(i) =arctany_i−yi−w

x_i−xi−w

−arctanyi−1−yi−w−1

xi−1−xi−w−1

, f(i)∈(−π, π), (2) where(x_i, y_i)are coordinates of the points, andwis size of the window, which is used in tangent angle calculation.

By comparing adjacent values of local curvature function, intersecting fibers can be dis- tinguished (each fiber contour would have different angle withx-axis). A suitable curve point for the skeleton is defined so that it minimizes the angle difference between the curve orientation, prior to arriving at the junction, and the curve orientation just after the junction [38].

Since in practice fiber segments usually do not intersect orthogonally and therefore have more than one intersection point, skeleton tracing (traversing each fiber starting from an end point and then proceeding from one pixel to its adjacent until reaching another end point [39]) is not able to process such segments as one fiber. Therefore, more global methods for segment grouping are needed, i.e. cost function optimization and optimal model fitting or grouping of segments according to perceptual principles [38].

2.5.1 Model fitting

The term geometric model fitting means representation of data objects as geometrical models. Model is easier to operate with, than a set of data points, and the task of model fitting is to determine the parameters of the models. The task of fitting the fibers in the pulp suspension images with geometric objects was not discussed in scientific papers before.

Y. Boykov et al. have investigated the problem of energy-based geometric multi-model

(16)

fitting [40], and the tests demonstrated that their approach significantly outperformed the RANSAC (RANdom SAmple Consensus) [41] based geometric model fitting methods that can be considered standard methods. The idea of the algorithm is to assign different labels to each observation (data point, object segment et al.) and select the optimal labeling which would correspond to the minimum value of energy function. The proposed approach PEaRL (Propose, Expand and Re-estimate Labels) [42] combines model sam- pling from data points as in RANSAC with iterative re-estimation of model parameters based on energy function minimization for each set of labels. For multi-model fitting each label corresponds to a candidate model, defined by a set of parameter values. A larger set of label proposals is more likely to contain models that approximate the real objects.

Model fitting can be applied for the tasks where objects can be represented as simple geometric primitives. Bubble detection problem was solved by fitting circles to the objects by Strokina [43]. Since fibers in the paper pulp are curved they cannot be approximated by simple geometric objects and therefore optimal cost function cannot be defined.

2.5.2 Perceptual grouping

Though fiber localization problem can be hardly solved mathematically (by model fitting), for human vision system solution of this problem is obvious. Perceptual grouping method groups the edge information depending on surrounding context, contextual and global information, as a human would do. According to the state-of-the-art survey of Papari et al. [44], various approaches to contour detection that have been proposed in the last two decades. Global operators of contour detection, including contour saliency computation and grouping of pixels together into contours using Gestalt psychology principles, that can be applied for fiber detection, are presented below.

Contour saliency. The first approach is similar to human psychology where the result is influenced by the presence of other stimuli in the surroundings. Psychological studies have inspired several authors to perform contour detection in two steps: first, local edge strength is computated; second, edge map is inhibited or enhanced depending on the surrounding context. The result is a more informative quantity of edges that is called contour saliency [44].

Several approaches for computing contour saliency have been developed: tensor voting, label relaxation et al.

(17)

• Relaxation labeling. It is a powerful probabilistic graph-based framework for the computation of saliency. The problems are cast as the assignment of labels to the elements of the scene from a set of possible labels (fiber indices for the fiber localization problem). Labels that violate consistency according to the predefined criteria are iteratively removed from the tokens until convergence [45].

• Tensor voting.Guy and Medioni [45] have introduced the tensor voting framework for computing contour saliency. Their method is considered to be one of the best known approaches for facilitation [46]. It was presented in details in [45]. At first, their approach was applied to curves [47], later it was expanded to surfaces as well [48]. A saliency-enhancing operator capable of highlighting features (edges, junctions etc.), which are considered psychologically ’important’, is introduced in their work. A vector field called the extension field is associated with each pixel of the image. It can be either chosen so that it ensures a smooth curve continuation, or oriented orthogonally to the original pixel’s direction [49]. Therefore, each pixel of the image is characterized by the distribution of the values from all other pixels of the image (the so called votes) and a saliency map is obtained from the statistics of this distribution [44]. The covariance matrix is considered and both saliency values and local edge orientation are obtained from its eigenvalues and eigenvectors. This approach is particularly useful for gap filling, contour closure, and in the presence of a large amount of noise [47].

The novelty of the approach is that there is no objective function, explicitly defined and optimized according to the global criteria. Tensor voting is performed locally and the saliency of perceptual structures is estimated as a function of the support tokens received from their neighbors. Tokens with compatible orientations that can form salient structures reinforce each other. The support of a token for its neighbors is expressed by the votes, cast according to the Gestalt principles of proximity, co- linearity and co-curvilinearity [45]. The authors are still working on this topic, and a new modification of their approach, Closed-form solution to tensor voting, has been published recently [50].

Capability of grouping fiber segments with compatible orientation and local representation of fibers makes tensor voting a suitable approach for solution of the fiber localization problem.

Grouping pixels together into contours according to Gestalt principles. Principles of Gestalt psychology which guide perceptual organization in the Human Visual System include similarity, proximity, good continuation, simplicity, closure, collinearity, and co-

(18)

curvilinearity [51]. Graph theory is a natural mathematical framework for this purpose:

nodes of a graphGcorrespond to edge pixels and interconnections between nodes model potential links between contiguous edges. A contourCis represented by a subgraph ofG and it is detected by minimizing a cost functionC computed overG[44].

Since formulation of the cost function for fiber segment grouping is a challenge, tensor voting occurs to be the optimal method for solution of the current problem.

2.6 Characterization

When the fibers have been localized and parameterized, their characteristics can be evaluated. In [8], Hirn provided an overview of fiber properties, which are commonly summarized under the term "fiber morphology" and can be considered a starting point in fiber characterization. The morphology, i.e. the structural appearance of fibers, is commonly described by five parameters: length, width, coarseness, kink, and curl.

A set of parameters, characterizing the pulp suspension, is:

• Suspension consistency (fiber concentration).

• Proportion of fines.

• Fiber length distribution.

• Fiber width distribution.

• Fiber curvature distribution.

• Presence of kinks, cells.

The fiber length is defined as the fiber contour length,L, or as the end-to-end (projected) length, l [11]. Diagrammatic representation is shown in Figure 4. Lcan be found as a length of an approximating spline or as a distance between pixels in the skeleton, whilel corresponds to the distance between end points.

While analyzing fibers and their quantity, it is important to take not an average mean, but a length weighted mean (to take fines into account) [8]. The length weighted average fiber length [10] is calculated as:

Lˆ =

PL_i²

PL_i , (3)

(19)

Figure 4. A diagrammatic representation of a curled fiber with two definitions of length [4].

whereLˆ is the length-weighted average length,L_i is the length of the fiber.

Fiber kink and curl provide information about the curvature of the fiber. Th curl describes the degree of curvature (non-straightness) of a fiber, but it gives no information if a fiber curves gradually or abruptly [8]. Fiber kink denotes small regions of very high curvature, i.e. sharp bends, along the fiber (as shown in Figure 5). It indicates fiber deformations, mostly caused by a mechanical damage.

Figure 5.Fibers with different kink but the same curl [8].

A curl index of zero indicates that no curl is present [4]. The curl index is calculated for each individual fiber as:

CI = L

l −1. (4)

Kink is an abrupt change in the fiber curvature. The most widely used definition for a kink is the Kibblewhite⁰s equation [4]. Kibblewhite established that large kinks in fibers had more impact on paper properties than small kinks. Therefore, his equation places greater

(20)

importance on larger angled kinks:

KinkIndex = 2N(21−45)+ 3N(46−90)+ 4N(91−180)

L_total , (5)

where N_i,j is the number of fiber segments with an edge orientation change of i − j degrees,L_total is the sum of fiber lengthsL.

As it was investigated in [4], fiber length and coarseness have a large impact on paper properties such as tensile index, folding endurance, improved formation et al. Fiber curl and kink have been shown to affect tensile stiffness, tear index, porosity, absorbency et al.

(21)

3 DESCRIPTION OF THE DEVELOPED METHOD

The developed method of fiber characterization in the pulp suspension images consists of three steps. Firstly, foreground data is separated from the image background. Secondly, the segments, obtained at the previous step, are grouped so that distinct data segments corresponding to the same fiber are joint. Finally, when each fiber of the input image has been represented as a curve, characteristics of the fibers are estimated. Overall scheme of fiber image processing is given in Figure 6.

Figure 6. Fiber characterization scheme.

3.1 Segmentation

At the first stage of pulp image processing the fibers are segmented from the background using the ridge edge detection method [29]. The approach comprises direction sensitive filtering, selection of the dominant orientation, non-maximum suppression, and hysteresis thresholding. The segmentation method takes a grayscale image with dark fibers and returns an image with thin bright lines, corresponding to the fibers in the original image.

3.1.1 Direction sensitive filtering

In order to detect dark fibers on the light background appropriate filter needs to be designed. The filter should give maximum response at the center of a fiber. Laplacian of

(22)

Gaussian (LoG) filter suggested by Marr and Hildreth [52] was selected (see Figure 7).

30 40 50 60 70 80 90 100

−0.15

−0.1

−0.05 0 0.05 0.1 0.15 0.2 0.25

Figure 7.Laplacian of Gaussian filter.

The LoG filter, shown in Figure 7, gives the maximum response for bright objects, for detecting dark fibers on bright background it has to be inverted. The LoG operator, defined as [52]

∂²

∂²xG_σ(x,y) =−x²−σ² σ⁴ ·e⁻ ^x

2

2·σ2 (6)

is computed for different values of xpoints with constant scale σ. This vector of LoG values is repeated, so that a square matrix for the LoG kernel is obtained. σ determines the width of the center part of the LoG operator and depends on width of fibers that need to be segmented. Size of the filter kernel influences the sum of all elements of the kernel:

if the sum is equal to zero, than the result of convolution for homogeneous regions is zero.

In order to avoid abrupt changes on the edges of LOG filer response and to reduce false alarms, the Gaussian filter is applied to the obtained operator [52]. For this purpose a kernel of the same size as LoG is created, that contains values for the first-order derivative with the same sigma:

G_σ(x,y) =− 1

σ·√

2·π ·e⁻ ^x

2

2·σ2 (7)

This kernel is rotated in order to be orthogonal to the LoG kernel. Convolution of these two square kernels of the same size results in the horizontally oriented filter. Sum of all the elements of this combined filter remains zero.

As fibers in the pulp suspension are randomly oriented, filters for their detection on the images should be also oriented. Input grayscale image is filtered with eight filter masks, oriented atπ/8 degrees angle. Filtering gives 8 grayscale direction-filtered images as a result. Filter masks are demonstrated in Figure 8.

(23)

Figure 8.8 filter masks, rotated byπ/8 angle.

3.1.2 Selection of dominant orientation

For 8 images, filtered in different directions, orientation that gives a maximum response is selected. Filter responses for each pixel are investigated in all directions, the maximum value and the corresponding orientation are saved in 2 result matrices (see Fig- ure 9). Therefore, from the 8 grayscale filtered images one image with a maximum filter responses and orientation of pixels is obtained.

Figure 9.Example of the dominant orientation selection.

3.1.3 Non-maximum suppression

The values of the neighboring pixels are suppressed, and the edges of the result image are represented as thin lines. For each pixel two neighboring pixels are taken into account.

If the value of the current pixel is the maximum, values of the neighboring pixels are set

(24)

to zero (Figure 10). The neighboring pixels are chosen in the direction, orthogonal to the orientation of the pixel.

(a) (b)

Figure 10. Examples of non-maximum suppression: (a) No suppression (center pixel is not the maximum); (b) Suppression of non-maximum intensities.

Each pixel can have one of 8 orientations, therefore it is not always straightforward to define the neighboring pixels in the orthogonal direction. For the pixels with ⁽²ⁿ⁺¹⁾₈ π orientation, where n = 0,1,2,3, the neighboring pixels are located ’between’ the real pixels of the image. Bilinear interpolation is used to solve this problem and values of the neighboring pixels in ⁽²ⁿ⁺¹⁾₈ π direction are interpolated between two nearest real pixel values in the selected direction [53].

3.1.4 Hysteresis thresholding

In order to select the edges that belong to the fibers, thresholding is applied to the grayscale non-maximum suppression output. In hysteresis thresholding two threshold values are used [54]. At the first step pixels with values, that are greater than the upper threshold value are set to maximum value. At the next step the maximum value is assigned to the pixels with values, that are greater than the low threshold value and belong to the same connected components, as the pixels, selected on the previous step. All the other pixels are set to zero. Example is shown in Figure 11.

3.2 Fiber segment grouping

The result of segmentation is an oriented edge-map. Since the output of the non-maximum suppression may contain disjoint parts, these segments need to be grouped. Perceptual grouping with tensor voting framework makes it possible to continue the fibers and to

(25)

(a) (b) (c)

Figure 11.Examples of hysteresis thresholding: (a) Original image; (b) Result of the first step of thresholding; (c) Result of the second step of hysteresis thresholding.

bridge the gaps based on Gestalt psychology principles of proximity and good continuation. The tensor voting methodology consists of four steps [55]: input data encoding, tensor voting, tensor field decomposition, and feature extraction. Firstly, input data is represented in tensor notation (encoded) so that each input token (object with a value, initially - point position) becomes a tensor token (oriented object). Secondly, distribution of votes from sparse tokens and accumulation of votes at every location (tensor voting) is performed. Thirdly, the tensor field is analyzed and saliency maps are extracted. Finally, salient features (curves and junction points) are extracted from the saliency maps. The overall scheme of the tensor voting approach [55] is presented in Figure 12.

Encode Encode

Tensor Voting Tensor Voting

Decompose Decompose

Feature extraction Feature extraction

Input tokens (sparse) Input tokens

(sparse)

Tensor tokens Tensor tokens

Saliency tensor field (dense) Saliency tensor

field (dense)

Junction saliency map Junction saliency

map Curve saliency

map

Curves Curves

Figure 12.Overall scheme of tensor voting [55].

(26)

3.2.1 Tensor representation

In the 2-D case the salient features to be extracted are curves and junction points. The saliency (perceived importance) of each type of the perceptual structure is indicated with a symmetric non-negative definite tensor, which is geometrically equivalent to an ellipse.

Tensor is a generalization of a scalar and a vector, it can represent points with no associated direction in contrast to vector, which would have zero size in case of direction information absence. Shape of the tensor indicates the type of structure (stick or ball component) and its size - the saliency of this information [45]. Tensor (or ellipse) can be described by a 2×2 eigensystem, where the eigenvectorse₁ ande₂give the ellipse orientation and the eigenvaluesλ1 and λ2 (λ1 > λ2) represent the ellipse size. The tensor is internally represented as a matrix S:

S =λ₁ ·e₁e₁^T +λ₂·e₂e₂^T (8) An input token that represents a point on the curve element is encoded as a stick tensor, wheree₂represents the curve tangent ande₁the curve normal, whileλ₁ = 1andλ₂ = 0.

3.2.2 Voting with tensors

Tokens are communicating with each other through a voting process, where each token generates a vote at each site in its neighborhood. The token, which generates votes, is a voter, and the tokens at its neighborhood, which receive votes, are receivers. The size and shape of the neighborhood, the vote strength and orientation are stored in pre-calculated voting fields (kernels). In the current problem, when all the points have normals, stick voting field can be used. The voting field is generated based on a single parameter - a scale factor σ. Vote orientation corresponds to the best (smoothest) local curve continuation from the voter to recipient. Strength of the voting fieldV S(−→

d)decreases with distance

|−→

d|from voter to recipient, and with curvatureρ [49]:

V S(−→ d) =e⁻

|→− d|2

+ρ2

σ2 (9)

In Figure 13 (a) generation of votes for a 2-D stick field with Equation 9 is given.

A tensor in pointP where the curve information is locally known (depicted by the curve normal−→

N_P ) casts votes at its neighborQ. An orientation of the voting field that ensures a smooth curve continuation (through a circular arc) from voterP to recipient Q is chosen.

(27)

(a) (b)

Figure 13.Voting fields: (a) Vote generation; (b) 2D voting field [49].

Thus obtained curve normal−→

N needs to be propagated by the vote from P to Q. The vote V_stick(−→

d)received from P at Q is encoded as a tensor as follows:

V_stick(−→

d) = V S(−→ d)·−→

N−→

N^T (10)

where−→

d =Q−P.

It is important to note that the vote strength at points Q⁰ and Q⁰⁰ is smaller than at Q because Q’ is farther, andQ⁰⁰ demands a higher curvature thanQ[49]. In Figure 13 (b) the result stick voting field with its color-coded strength is shown.

3.2.3 Vote accumulation

The votes are collected through simple tensor addition (sum of matricesV S(−→

d)) at every point. Tokens that lie on a smooth curve reinforce each other after the voting process, and the tensors deform according to the dominant orientation of the result ellipse. Each tensor encodes the local orientation of the geometric features such as curves (given by the tensor orientation) and confidence of this knowledge (also calledsaliency, given by the tensor shape and size). For a general tensor, its curve saliency is given by (λ₁−λ₂) and the curve normal orientation bye₁, while its point saliency is given byλ₂ [49]. However, some noise can be obtained for tokens that receive very little support. Curve and point saliency are stored in corresponding saliency maps.

(28)

3.2.4 Curve extraction

Curves are extracted as the local maxima of curve saliency map, while junctions are in- ferred as the local maxima of the junction saliency map [55]. Since calculating votes for every location in the volume containing the data points is pointless and impractical, sur- face and curve extraction begins from seeds, location with highest saliency, and voting is performed only towards the directions indicated by curve tangents [56].

Salient points are arranged in descending order of their saliency. Each point, if it has not been already assigned to another curve, can be a start point. For each start point the most salient neighbor in the curve tangent direction is selected and recursive curve growing starts for this salient neighbor. The process stops when an end point is reached, or when the saliency is smaller than the threshold. The curve is grown from the other side of the start point in the curve tangent direction by the same recursive algorithm. At the junctions the curve continues growing with the same tangent direction as the previous point, that was added to the curve before the junction area. This assumption is made because of uncertainty of curve tangent in the junction area, where tensors have shape of a ball, not stick.

The overall process of tensor voting is summarized in Algorithms 1 and 2.

Algorithm 1: Tensor voting Input:oriented edge-map.

Output:saliency maps.

Initialize saliency tensor field.

for alloriented points of edge-mapdo

Generate an ellipse with current orientation at the corresponding location of sparse saliency tensor field.

end for

for alltokensdodistribute votes

Rotate the voting field corresponding to the token dominant orientation.

Add the voting field to the saliency tensor field (summarize S matrices at each point).

end for

Represent saliency tensor field with tensor maps.

for alltokensdo

Store (λ₁ −λ₂) saliency at the corresponding point of the curve map.

Storeλ2 saliency at the corresponding point of the junction map.

end for

Curve extractionalgorithm.

(29)

Algorithm 2: Curve extraction Input:

saliency tensor field;

curve saliency map;

junction saliency map;

binary image with curve endpoints.

Output:

extracted fiber curves.

Arrange curve saliency point in the descending order of their saliency.

Set pixels in an-neighborhood of each salient point to zero.

for allsalient pointsP dostart growing curves

while(P is not an end point) and (saliency ofP is greater than the threshold)do dir= dominant tensor orientation atP.

CURVE GROWING fromP in the most salientdirdirection.

CURVE GROWING fromP in the most salient−dirdirection.

end while end for

functionCURVE GROWING(saliency maps, seed pointP, tangent directiondir) current seed pointcurP =curP.

ifcurP is not at the fiber junctionthen

dir= dominant tensor orientation atcurP.

curP = the most salient point in thedirdirection.

AddcurP to the curve.

Savedir else

curP = the most salient point in thedirdirection.

AddcurP to the curve.

end if end function

(30)

3.3 Fiber characterization

The methods for fiber parameter estimation are based on curve representation (a set of points with coordinates x_i and y_i). Length of a fiber is sum of distances between the curve points:

L=

n

X

i=2

p(x_i−xi−1)²+ (y_i−yi−1)². (11)

Projected length of the fiber is estimated as a distance between the curve end points:

l =p

(xn−x1)²+ (yn−y1)². (12)

Curl index is calculated according to the equation (4):

CI = L l −1.

Width of the fiber (see Figure 14) is calculated using the coordinates of the curve points and a binary image (result of input image global thresholding with Otsu method). Fiber width is estimated at the equally distributed curve points, and the average width is considered to be the width of the current fiber. Width at each point is evaluated as the minimum width (in nonzero pixels) of the thresholded original image at the corresponding point.

Figure 14.Fiber width estimation.

(31)

4 EXPERIMENTS

4.1 Test data

The proposed approach to fiber detection and characterization was tested on a dataset of pulp suspension images, provided by the CEMIS-OULU laboratory of University of Oulu.

The images were captured with a microscopic CCD camera with 2.5x magnification. The algorithm was tested on 50 grayscale acacia pulp suspension images. Image dimensions are 600x800 pixels. Examples of the test images are presented in Figure 15.

Figure 15.Examples of the acacia images.

4.2 Ground truth data

The results of fiber detection and characterization are evaluated based on the manually created ground truth (GT) data. A non-expert indicated 7 key points for each fiber in the image. Coordinates of these points in most cases are not monotonic, and they may have similarxorycoordinates (for horizontal or vertical parts of fibers). This properties make curve model fitting quite a difficult task. Matlab Curve Fitting tool was used for this purpose, but no model suitable for precise fiber interpolation was found in the Carthesian coordinate system. In the current implementation fiber curves are represented as a set of line segments between the key points. Example of these key points and curve segments is demonstrated on a fragment of the original image in Figure 16. As an idea for future, the set of points can be possibly divided into monotonic subsets with further model fitting, so that a curve would be represented by several models. Another probable approach to curve representation is to use a specific coordinate system. A non-monotonic curve in a shape

(32)

of a circle can be represented in the polar coordinate system, the same solution could be possibly applied for the current task.

Figure 16.Example of the ground truth data.

GT fiber length is estimated according to Equation 11, curl index of fibers is evaluated according to Equation 4. Fiber width is evaluated as described in Section 3.3, with the only difference that GT key points instead of equally distributed points are considered for average width estimation.

4.3 Step by step input image processing

The original image is introduced in Figures 17(a) and 17(d). The result of filtering and the dominant orientation selection is presented in Figures 17(b) and 17(e). Shape and size of the filter kernel were determined depending on the size of the recognized objects and are presented in Table 1. In order to avoid the double responses of the filter, the highlighted parts of fibers are assigned the average background intensity. The hysteresis thresholding result in grayscale can be seen in Figures 17(c) and 17(f).

Table 1.Parameters of the developed method Name of the parameter value

Filter kernel size 15

Filter scale 1.9

Voting scale 10

The next stage of segmentation is tensor voting and curve extraction from the saliency

(33)

(a) (b) (c)

(d) (e) (f)

Figure 17.Examples of image segmentation: (a),(d) original images; (b), (e) filtered images; (c), (f) hysteresis thresholding results.

maps. A stick tensor field of a size, big enough for connecting discontinuities, but suitable for not connecting the distinct segments, was experimentally chosen and the value of the parameter is presented in Table 1. A curve saliency map is shown in Figures 18(a) and 18(c), while a junction saliency map is shown in Figures 18(b) and 18(d).

(34)

(a) (b)

(c) (d)

Figure 18. Examples of saliency maps: (a), (c) curve saliency maps; (b), (d) junction saliency maps.

The GT fibers are shown in Figures 19(a) and 19(c), the extracted curves are presented in Figures 19(b) and 19(d). Each curve is represented with a different color.

(35)

(a) (b)

(c) (d)

Figure 19. Examples of fibers (each fiber is marked with different color): (a), (c) GT fibers; (b), (d) extracted fibers.

At the stage of result evaluation the GT fibers (see Figure 20(a)) are compared with the computed curves 20(b). Correctly detected fibers (see Section 4.4 for detailed description) or true positives (TP) are colored in blue. False alarms (extracted fibers, which do not correspond to any GT fibers) or false positives (FP) are marked with yellow. GT fibers that were not detected by the method are false negatives (FN), and are colored in red. As not all the fines are marked in GT, the minimum fiber length was learnt from the GT, and the extracted fibers with length smaller than the minimum GT fiber length are not taken into consideration at the result evaluation stage.

An example of the fiber parameter evaluation is presented in Figures 21 and 22. In Fig- ure 21 (a) length of the ground truth fibers, calculated according to Equation 11, is shown.

In Figure 21 (b) length of the extracted fibers and length estimation error are presented.

(36)

(a) (b)

Figure 20. Examples of fiber images: (a) original image with GT fibers; (b) original image with extracted fibers (TP - blue, FP - yellow, FN - red).

The error of parameter evaluation is defined as follows:

P_err = P_GT −P_exp

P_GT , (13)

where P_GT is the GT value of the parameter, and P_exp is the experimental value of the parameter.

(a) (b)

Figure 21. Examples of fiber length: (a) GT fiber length; (b) fiber length evaluationLerrwith the developed method.

Since the GT fibers are marked manually by non-expert (see Figure 22 (a)), curl index can be calculated with an error (see Figure 22 (b)). This kind of errors is a drawback of current GT fiber representation with 7 key points.

(37)

(a) (b)

Figure 22. Examples of fiber curl index: (a) GT curl index; (b) curl index evaluationCI_err with the developed method.

4.4 Numerical evaluation of the results

The detected fibers and their characteristics are compared to the GT in order to evaluate the results. Firstly, for each GT fiber the corresponding extracted fiber is found. In order to do this a 2-pixel neighborhood area of each GT key point is looked through, and the indices of possible nearby extracted curves are stored. All the unique indices are sorted in such an order, that the first will be the most frequent index, and the last – the rarest.

Then extracted curve corresponding to each found index is compared with the GT fiber by carrying out the following binary comparison. Each GT fiber is defined by 7 points and can be represented as a binary mask, where each nonzero element corresponds to a key point. As the points were marked manually, their coordinates are not precise and cannot be directly compared with the extracted fibers. Each extracted fiber is a thin curve, so dilation is applied 10 times, and the curve becomes a wide strip. This strip can be also represented as a binary mask with nonzero elements, corresponding to the strip pixels.

Logical conjunction of these two masks indicates how many key points support the detected fiber (distance of the key point from the detected fiber is less than 10). If 6 or 7 key points support the detected fiber, then it is considered to be correctly detected and is marked as a TP. If the fiber does not fit to any strip, then it was not detected by the method, and it is a FN. Finally, the extracted curves that do not correspond to any GT fibers are considered to be FPs. The evaluation procedure is described in Algorithm 3.

(38)

Algorithm 3: Result evaluation Input:

set of manually marked GT fibers, each fiber is represented with 7 coordinate points;

set of extracted fibers, each fiber is represented with a set of points;

image with pixels corresponding to fiber labels (pixels of one fiber are marked with its index).

Output:

Number of TP, FP, TN.

for allGT fibersdo

Create a binary image with key points marked as nonzero pixels.

for allGT key pointsdo

Find values of pixels (indices of curves), located in a 2-pixel neighboring area of current point on the input image.

end for

Sort all the indices in order of decreasing frequency.

for allindicesdo

Create a binary image where pixels with current index are marked with ones, the other pixels are zero.

Apply morphological dilating to the binary image ten times (obtain wide fiber strip).

Create a binary image with key points marked as nonzero pixels.

Calculate logical conjunction of key point binary image and dilated fiber curve image.

iflogical conjunction has one pixel less than the key point binary imagethen Mark the fiber with current index as a TP

end if end for

ifcurrent GT fiber does not fit to any stripthen Mark the fiber as TN

end if end for

for allextracted fibersdo

iffiber is not marked as a TPthen Mark the fiber as FP

end if end for

(39)

Based on the classification results, performance of the method is evaluated using the following metrics:

• true positive rate (TPR) or sensitivity: T P R= _{T P}^{T P}_{+F N},

• false discovery rate (FDR):F DR = _{T P}^{F P}_{+F P},

• positive predicted value (PPV) or precision: P P V = _{T P}^{T P}_{+F P},

• statistical measure of test accuracy [57] (F-measure): F1 = _{2T P}_{+F P}^{2T P}_{+F N}.

The TPR measures the proportion of correctly detected fibers from the GT data. The FDR is designed to control the proportion of incorrectly detected fibers. The PPV specifies the proportion of correctly detected fibers from all the extracted fibers, while the F-score assesses the classification performance.

Estimation of the method performance, fiber parameters and their accuracy, calculated for the proposed method results on 50 images, are given in Tables 2 and 3. Accuracy of fiber parameters evaluation is computed as follows:

P_acc= (1−P_GT −P_exp

P_GT )·100%, (14)

whereP_GT is GT value of the parameter,P_exp is experimental value of the parameter. A full report on the results of the experiments is presented in Appendix 1.

Table 2.Evaluation of the method performance.

Characteristics Average value

per image

Number of GT fibers 50

Number of extracted fibers 63

Number of true positives (TP) 40

Number of false positives (FP) 25

Number of false negatives (FN) 9

True positive rate (TPR) or sensitivity 82%

Positive predictive value (PPV) or precision 62%

False discovery rate (FDR) 38%

F-measure of a test accuracy 70 %

Proportional execution time, required for

edge detection 10%

tensor voting 74%

curve extraction 16%

Average evaluation of processing time 251 seconds

(40)

Table 3.Fiber characteristics evaluation.

XXXX

Method

Parameter Average Average Average

length curl index width Proposed method 85 pixels 0.10 4.08 pixels

GT 99 pixels 0.16 4.68 pixels

Accuracy 84% 60% 85%

4.5 Discussion

The efficiency of the proposed method can be evaluated based on the experiments. 82%

of the GT fibers were detected correctly, though there were about 38% of false alarms, caused mostly by localizing one fiber as two separate parts at the intersection point. F- score that measures the effectiveness of fiber curve extraction [58] is considered to be 70%. For fiber length and width the error is about 15%, while curl index evaluation error is about 40%. Curl index is determined by the ratio of lengths, therefore, it is indirectly affected by the length evaluation error. The GT curl index cannot be evaluated precisely because of current GT fiber representation and non-expert GT data acquisition. These factors cause significant error of curl index estimation.

The weak point of the method is end and junction point inference: some end points are located close to the points of intersection, and therefore are detected incorrectly. That is the reason, why some curves were grown longer than necessary, or why the process of curve growing was stopped too early. Precision of the method can be increased by end point inference through first order tensor voting. The first order tensor (commonly known as vector) is used to encode the asymmetric orientation information at an end point.

Directional distribution of neighbors is computed at every point. First order moment of the received stick votes, which corresponds to the vector sum of the stick vote contribution, is collected instead of the second order moment collection (hence end-point votes are collected as vectors instead of tensors) [45].

As for the execution time, on 1.6 GHz Intel Celeron M processor it took approximately 4 minutes to process each image. In case it is necessary to decrease the processing time so that the method could characterize pulp during the process of paper production, more computational resources can be involved. Though the results, obtained with the developed method, can be used for average estimation of pulp consistency and stage of its prepara- tion, the method execution time is not suitable for in-line pulp analysis. Tensor voting is the slowest stage of fiber image processing and takes 74% of the processing time. Since tensor accumulation processes multiple receivers for each voter, a lot of tensor computa-

(41)

tions have to be done. Tensor size reduction would considerably reduce the required processing time, though, size of discontinuities, which are connected by the method, would be also decreased.

The proposed fiber detection method, based on direction sensitive filtering, non-maximum suppression, hysteresis thresholding, second order tensor voting, and curve extraction, is intended to detect elongated objects on homogeneous background, and can be applied not only for fiber detection in pulp images, but also for other objects with ridge edges, and for objects, that can be represented as geometric models.

(42)

5 CONCLUSION

The scope of this Master’s thesis was to detect elements in the pulp suspension and to characterize them according to the pulp properties. In the Section 2 the overall scheme of the fiber image processing was given, review of fiber morphology and existing methods of segmentation was made. The method based on direction sensitive filtering and tensor voting was proposed for image segmentation. This approach has not been used for fiber detection problem before.

The proposed method accomplishes the main goal of the research and shows that this method can give good results for segmentation of intersecting fibers. The developed method can be used for estimation of the pulp consistency as a part of in-line industrial measurements process. One possible way of implementation is an embedded system or a computer, that receives pulp images and processes them in order to get more information about current state of the pulp.

The developed method was successfully tested on 50 acacia pulp images. Ground truth data for result evaluation was created manually by non-expert. As a result, performance of the method implementation and suspension features were evaluated and statistical characteristics were estimated. The proposed method achieved considerably good results of fiber detection with 82% of sensitivity. One image is analyzed for about 4 minutes. How- ever, the precision is only about 62%, so the algorithm still needs improvement. There are several points where performance of the method can be improved. The images could be possibly captured without blinks. As a part of future work, collecting more extensive information about ground truth fibers can be considered, so that accurate ground truth fiber curl index could be estimated. As for the result of tensor voting, it can be improved by first order tensor voting for precise definition of curve end points. Algorithm of curve extraction can be also improved at the intersection points.

(43)

REFERENCES

[1] KnowPulp Learning Environment. www.knowpulp.com, 2012.

[2] Machine Vision and Pattern Recognition Laboratory (MVPR).

http://www2.it.lut.fi/mvpr/, 2012.

[3] Pulp Vision Project. http://www2.it.lut.fi/project/pulpvision, 2011.

[4] P. Sutton, C. Joss, and B. Crossely. Factors affecting fiber characteristics in pulp. In Tappi journal, 2000.

[5] The Paper Company Ltd (PaperCo). http://www.paperco.co.uk/, 2012.

[6] J. Levlin and L. Soderhjelm. Pulp and Paper Testing. Fapet Oy, 1999.

[7] G. Robertson, J. Olson, and R. Seth. Measurement of fiber length, coarseness, and shape with the fiber quality analyzer. Tappi journal, 82:93–98, 1999.

[8] U. Hirn and W. Bauer. A review of image analysis based methods to evaluate fiber properties. Lenzinger Berichte, 86:96–105, 2006.

[9] M. G. Carvalho, P. Ferreira, and A. Martins. A comparative study of two automated techniques for measuring fiber length. Tappi journal, 80:137–142, 1997.

[10] G. Don, N. Sutherland., and W. Rantanen. Comparison of fiber length analyzers. In TAPPI Practical Papermaking Conference, 2005.

[11] R. Trepanier. Automatic fiber length and shape measurement by image analysis.

Tappi journal, 81:152–154, 1998.

[12] A. Bartl and D. Pico. Characterization of short fibers. In The 9th International Conference on Chemical and Process Engineering, 2009.

[13] R. Woods and R. Gonzales. Digital Image Processing. Pearson Prentice Hall, 3rd edition, 2008.

[14] QImaging cameras. www.qimaging.com. 2012.

[15] M. Sezgin and B. L. Sankur. Survey over image thresholding techniques and quan- titative performance evaluation. Journal of Electronic Imaging, 13:146–165, 2004.

[16] J. Weszka and A. Rosenfeld. Threshold evaluation techniques. IEEE Transactions on Systems, Man and Cybernetics, 8:622–629, 1978.

(44)

[17] N. Otsu. A threshold selection method from gray-level histograms. IEEE Transac- tions on Systems, Man and Cybernetics, 9:62–66, 1979.

[18] N. Kapur, P. Sahoo, and A. Wong. A new method for gray-level picture thresholding using the entropy of the histogram. Computer Vision, Graphics, and Image Processing, 29:273–285, 1985.

[19] J. P. Serra and P. Soille, editors. Mathematical Morphology and its application to Image processing, volume 6671 of Lecture Notes in Computer Science. Kluwer Academic Publishers, 1994.

[20] J. Serra. Image Analysis and Mathematical Morphology. Academic Press, 1982.

[21] N. Amenta, S. Choi, and R. Kolluri. The power crust, unions of balls, and the medial axis transform. Computational Geometry: Theory and Applications, 19:127–153, 2000.

[22] H. Breu. Linear time euclidean distance transform algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17:529 – 533, 1995.

[23] A. Ferreira and S. Ubeda. Computing the medial axis transform in parallel with eight scan operations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 21:277 – 282, 1999.

[24] N. Naccache and R. Shinghal. An investigation into the skeletonization approach of Hildich. Pattern Recognition, 17:279–284, 1984.

[25] J. Canny. A computational approach to edge detection. IEEE Transactions on Pat- tern Analysis and Machine Intelligence, 8:679–698, 1986.

[26] L. Qin and J. Jian. The recognition of linear objects based on edge detection and skeleton tracking. Third International Symposium on Information Science and En- gineering, pages 379–381, 2010.

[27] G. Hui and Q. Shubo. Pulp fiber corner detection arithmetic based on local curvature function. InProceedings of the 27th Chinese Control Conference, pages 160–163, 2008.

[28] T. Lindeberg. Scale-space theory: A basic tool for analysing structures at different scales. Journal of Applied Statistics, 2:224–270, 1994.

[29] T. Lindeberg. Edge detection and ridge detection with automatic scale selection.

International Journal of Computer Vision, 30:117–154, 1998.

Characterization of fiber and vessel elements in pulp suspension images

CHARACTERIZATION OF FIBER AND VESSEL ELEMENTS IN PULP SUSPENSION IMAGES

ABSTRACT

PREFACE

CONTENTS

ABBREVIATIONS AND SYMBOLS

1 INTRODUCTION

1.1 Background

1.2 Objectives and Restrictions

1.3 Structure of the Thesis

2 PAPER PULP CHARACTERIZATION

2.1 Pulp and papermaking

2.2 Fiber analyzers

2.3 Analysis of the pulp suspension images

2.4 Fiber segmentation

2.5 Fiber localization and parameterization

2.6 Characterization

3 DESCRIPTION OF THE DEVELOPED METHOD

3.1 Segmentation

3.2 Fiber segment grouping

3.3 Fiber characterization

4 EXPERIMENTS

4.1 Test data

4.2 Ground truth data

4.3 Step by step input image processing

4.4 Numerical evaluation of the results

4.5 Discussion

5 CONCLUSION

REFERENCES