Stochastic Texture Analysis for Monitoring Stochastic Processes in Industry
Jacob Scharcanski a a UFRGS–Universidade
Federal do Rio Grande do Sul
Instituto de Inform´ atica Av. Bento Gon¸calves, 9500. Porto Alegre, RS, Brasil 91501-970
[email protected] Abstract Several continuous manufacturing processes use stochastic texture images for quality control and monitoring. Large amounts of pictorial data are acquired, providing both important information about the materials produced and about the manufacturing processes involved. However, it is often difficult to measure objectively the similarity among such images, or to discriminate between texture images of materials with distinct properties. The degree of discrimination required by industrial processes sometimes goes beyond the limits of human visual perception. This work presents a new method for multi-resolution stochastic texture analysis, interpretation and discrimination based on the Wavelet Transform. A multi-resolution distance measure for stochastic textures is proposed, and applications of our method in the non-woven textiles industry are reported. The conclusions include ideas for future work.
Key words: stochastic textures, wavelets, anisotropy, nonwoven textiles
1
The author thanks CNPq (Brazilian Research Council) for financial support, and
Mr. Osmar Machado (Riocell, Brasil) for providing experimental data; thanks are
Preprint submitted to Elsevier Science
27 September 2004
1
Introduction
In several continuous processes, static and dynamic stochastic texture images are acquired and used in quality control [1]. Often, industrial machine operators try to visually interpret stochastic texture images, and estimate the manufacturing process condition using their experience in the field. This empirical approach is subjective, and prone to failure, mainly because the human vision is limited in terms of its ability to distinguish between stochastic textures [2]. Despite advances in texture representation and classification over the past three decades [3] [4], the problem of stochastic texture feature interpretation and classification remains a challenge for researchers [5], and for a large segment of the industry [2]. There have been a variety of methods of extracting texture features from textured images, e.g. geometric, random field, fractal, and signal processing models for textures. A significant part of the recent works on textures concentrate on statistical modelling [3], which characterizes textures as probability distributions, and uses statistical theories to formulate and solve texture processing problems mathematically. Wavelet-based texture characterization has attracted attention recently because of its usefulness in several important applications, such as texture classification [4] [6], and texture segmentation [7]. Several approaches have been proposed to extract features in the wavelet domain with application in texture analysis, such as: (a) wavelet energy signatures, which were found useful for texture classification [4]; (b) second-order statistics of the wavelet transform were also used to improve the accuracy of texture characterization [8]; and (c) higher order dependencies of wavelet coefficients were studied for texture analysis [3] [7] [9]. These wavelet-based apdue also to Professor Roberto da Silva (Instituto de Informatica, UFRGS, Brasil) and to Professor Robin T. Clarke (Instituto de Pesquisas Hidraulicas, UFRGS, Brasil) for advice and useful discussions.
2
proaches have been found more effective than other methods based on secondorder statistics or random fields, which analyze textures in a single resolution without considering the human visual perception of textures [3] [2]. Most of the work on texture in the wavelet domain has concentrated on the analysis of visual textures, and were not designed for stochastic texture analysis. For example, often feature extraction is carried out assuming sub-band independence at each resolution (e.g. [4] [6]), which is not verified experimentally. Also, the methods based on higher order statistics generally do not make explicit relevant stochastic texture features, that are important for industrial applications, where process conditions are estimated based on specific texture parameters (e.g. [3] [7] [9] [2]). In this work, a multi-resolution scheme for stochastic texture representation and analysis is proposed. We begin by describing how we measure the image gradients in multiple resolutions. Based on this technique, the local grayscale variability and texture anisotropy are measured in multiple resolutions. Next, a multi-resolution distance measure for stochastic textures is introduced. Finally, we present some applications, experimental results and conclusions.
2
Our Proposed Texture Representation
In this work, we emphasize specific stochastic texture features that could be used to facilitate texture interpretation, and the discrimination between distinct process conditions
2
.
Our method relies on multiple resolution texture gradients and their magnitudes. To estimate the local gradients in multiple resolutions, we apply the redundant two-dimensional WT proposed by Mallat and Zhong [10]. The co2
One, or more samples, can be used as reference images, indicating a particular
process condition.
3
efficients W21j f (x, y) and W22j f (x, y) represent the details in the x and y directions, respectively, and approximate the image gradient at the resolution 2j . Since we are dealing with digital images f [n, m], we use the discrete version of the WT [10], and the discrete Wavelet coefficients are denoted in this work by W2ij f [n, m], for i = 1, 2. The gradient magnitudes at the resolution 2j are computed from q
M2j f [n, m] =
(W21j f [n, m])2 + (W22j f [n, m])2 .
(1)
2.1 Texture Directionality in Multiple Resolutions
Texture directionality, i.e. anisotropy, is an important parameter in the manufacture of foil-like materials. It correlates well with several mechanical and transport properties of such materials, as well as with commonly-monitored manufacturing variables. In this work, the occurrences of the coefficients W21j [n, m] and W22j [n, m] are approximated by Gaussian distributions. The normal plots in Figures 2(a) and (b) show that the Gaussian model can represent a wide range of wavelet coefficient values in stochastic textures, with some deviation from Gaussian occurring in the tails of the distribution, which was confirmed in a large-sample test set. It was also verified that, in general, wavelet coefficients associated with the same image position [m, n] in different sub-bands (i.e. W21j [n, m] and W22j [n, m]), at the same resolution 2j , are correlated (i.e. are not independent). Therefore, we represent the joint distribution of wavelet coefficients by 1 2 12 a bivariate Gaussian G12 2j (W2j [n, m], W2j [n, m]), denoted simply by Gj . The
iso-probability curves of the bivariate Gaussian G12 j are typically elliptic for anisotropic samples, and tend to be circular for samples that are isotropic. The joint coefficient distribution G12 j determines two orthogonal axes of ex4
tremal variance, that coincide with the directions of the eigenvectors vmax and vmin of the covariance matrix. In order to estimate the covariance matrix at resolution 2j , let us denote the means of the coefficients W21j and W22j by µ1 and µ2 , respectively. The covariance matrix is then calculated as follows: h
C2j = E (W21j − µ1
´³
2 σ1 W22j − µ2 )T = i
ρ12 σ1 σ2
ρ12 σ1 σ2 σ22
.
(2)
where, ρ12 is the correlation coefficient of W21j and W22j , and σ1 and σ2 are the standard deviations of W21j and W22j , respectively. The shape and orientation of the coefficient distribution, at resolution 2j , is described by a Gaussian ellipse with covariance matrix C2j , mean [µ1 µ2 ]T . 0
The orientation θ of its main semi-axis can be obtained from [11]: Ã 0
tan(2θ ) =
! 2ρ12 σ1 σ2 ) σ12 −σ22
−1
and θ = tan
(3)
!
à 0
,
0
tan(2θ ) tan(2θ0 )+2
. Figures 1(c) and (d) show Gaussian ellipses fit-
ting the actual angular distribution for isotropic and anisotropic texture images. Both ellipses were calculated at the finest resolution 21 , and provide a visual indication of the degree of anisotropy of each texture. It shall be noticed that often the joint distribution of coefficient values do not have zero means (i.e. µ1 6= 0 and µ2 6= 0), as illustrated in Figure 1(d). In order to measure quantitatively the distribution eccentricity, we calculate the eigenvectors vmax and vmin from (2), as well as the corresponding eigenvalues λmax and λmin . The eigenvalues λmax and λmin define the semi-axes of a Gaussian ellipse aligned with the eigenvectors directions. The distribution eccentricity e is given by the ratio of the eigenvalues: e=
λmax , λmin
(4)
5
which provides an estimate for the texture anisotropy at scale 2j . For example, the measured eccentricity e for the isotropic sample in Figure 1(a) is 1.008, and for the anisotropic sample in Figure 1(b) is 1.227. It should be noticed that the eccentricity can be calculated in multiple resolutions.
2.2 Texture Local Graylevel Variability in Multiple Resolutions
The texture local graylevel variability in multiple resolutions encodes important information about local density variability, and consequently about the material homogeneity. As discussed in the previous section, the stochastic texture image coefficients W21j f [n, m] and W22j f [n, m], for j = 1, 2, ..., J, or simply, W21j f and W22j f , may be considered as Gaussian distributed (See Figure 1). As a consequence, the corresponding distribution of magnitudes M2j f =
q
(W21j f )2 + (W22j f )2 , de-
noted simply by M2j , may be approximated by a Rayleigh probability density function [12]: −M 2
2j M2j Rj (M2j ) = j 2 e [2bj ]2 , [b ]
(5)
where, q j
b =
2 2σM 2j
4−π
,
(6)
and σM2j is the standard deviation of the gradient magnitudes at scale 2j . Figures 1(e) and (f) illustrate the Rayleigh model fitting the actual gradient magnitudes distribution. In fact, stochastic textures may vary in anisotropy, local variability, or both, at different resolutions, and a wide range of stochastic textures can be represented in terms of these features. Typically, samples of similar textures also have similar bj parameters, and this information may be 6
used for texture discrimination and similarity matching, as discussed later. A distance measure to compare these texture representations is introduced next.
3
A Stochastic Texture Distance Measure
In our approach, at each resolution, a texture image Ii is represented by a histogram with k disjoint intervals of equal length, {S1 , S2 , ..., Sk }. Under these circumstances, Do and Vetterli [6] showed that the Kullback-Leibler distance ranks, or classifies, a texture I consistently with Maximum Likelihood rule. The empirical distributions (i.e. data histograms), usually require large representation overheads (i.e. storing large numbers of histogram bins). Therefore, as detailed before, we model these distributions, at each resolution, by bivariate Gaussian G12 j and Rayleigh Rj probability density functions. The Kullback-Leibler distance between two bivariate Gaussian distributions G(C1 , M1 ) and G(C2 , M2 ) is given by [13] : 1 1 DGauss2 = log (det(C2 )/det(C1 )) + trace(C1 C2−1 ) + (M1 − M2 )T C2−1 (M1 − M2 ),(7) 2 2 where, C1 and C2 are covariance matrices, M1 and M2 are the means of the bivariate distributions, and det(C) denotes the matrix determinant. On the other hand, the Kullback-Leibler distance between two Rayleigh distributions R(b1 ) and R(b2 ) is : Ã
DRayl
1 b2 b22 − +1 = − log π b1 b21
!
(8)
where, b1 and b2 are the Rayleigh parameters of the textures, respectively. The Kullback-Leibler distance has some nice properties. Its convexity guarantees that a minimum exists, and in order to calculate the Kullback-Leibler distance from multiple data sets, such as from different sub-bands or feature 7
sets, we can use the chain rule [6]. This rule states that the Kullback-Leibler distance DKL between two pairs of conditional probability density functions is simply the sum of the distances from corresponding probability density function marginals. Therefore, the proposed multi-resolution stochastic texture distance measure is : Dj =
J X
j
j
(1 − α)DGauss2 + αDRayl
(9)
j=1
j
j
where, DGauss2 (²[0, 1]) and DRayl (²[0, 1]) are the symmetric and normalized Kullback-Leibler distances between bivariate Gaussian and univariate Rayleigh distributions at each resolution j, and α(²[0, 1]) is a parameter that j
j
controls the weights attributed to DGauss2 and DRayl .
4
Experimental Results
Stochastic texture images are widely used in manufacture of foil-like materials such as non-woven textiles, paper, polymer membranes, conductor and semiconductor coatings. Paper and non-woven textiles have a definite ”grain” caused by the greater orientation of fibers in the machine direction and by the stress/strain imposed during pressing and drying. The directionality of such materials affect substantially their physical properties. Due to fluctuations and irregularities in systems operation, there are corresponding variations in the makeup of fibers and orientation of fibers across the machine (i.e. web) [14]. Consequently, the physical properties may vary in both the MD (machine direction) and CD (cross direction), and often is necessary to sample and test at several locations across the machine. The orientation of paper and non-woven textiles can be determined by testing the mechanical resistance to tearing in the MD and CD. During the past two decades, automated non8
woven testing became a new trend in industry, and a popular orientation test is the ”TSI/TSO” test [14], which measures orientation by measuring the sound propagation speed in the samples, in different orientations. All the orientation tests mentioned above are destructive, or require contact with the sample. Non-destructive and non-contact testing has been a challenge for researchers [2]. To illustrate the performance of our approach in stochastic industrial texture classification, we used 315 distinct β-radiographic images of non-woven textile and paper samples, obtained from a standard industrial image database for non-woven textiles and paper [15]. We chose 9 samples of each type, from a variety of forming machines, with different furnish and grammage
3
(e.g. among
the 315 samples we have samples of headbox handsheets, repulped machine sheets, standard handsheets, board handsheets, gap formers, Fourdrinier formers, speedformers, tissue papers, and glass fiber mats). All these images have a resolution 140 × 140 pixels, with a spatial resolution of 0.2 × 0.2mm2 /pixel [15]. Each type of texture indicates a particular production condition (i.e. in terms of operating parameters and furnish). The sample textures representing different operating conditions were classified, and the system current condition could be identified. Potentially, this approach could help non-woven textile industry operators identify the current system condition, who usually rely on “ad-hoc” texture interpretation methods. In our classification experiments we compared our approach (i.e. Equation 9) with the Spatial Gray Level Dependence Method SGLDM
4
and Gabor filters
5
, in both cases the Euclidean distance was used as a metric. We also com-
3
Grammage is defined as mass density in g/m2 Each sample is represented by gray-level co-occurrence matrices for the 0o ,45o
4
and 90o orientations; matrices are described by 5 features: contrast, homogeneity, angular second moment, entropy and variance. 5 The Gabor features are the mean magnitudes of the 24 sub-bands (i.e. 4 scales
9
pared our approach with the performance of its two individual components, i.e. a bivariate Gaussian model for wavelet coefficients, and the Rayleigh model for gradient magnitudes, using as a metric the symmetric and normalized Kullback-Leibler distance for bivariate Gaussians and Rayleighs, respectively. Table 1 shows the correct classification rates obtained (as percentages), indicating that with our approach we obtain the highest correct classification rate (with α = 0.15). With the other texture representation approaches we obtained lower correct classification rates, at a similar, or even higher, computational costs, as shown in Table 1. Computational costs were evaluated in terms of the running times, i.e. seconds, of the Matlab routines implementing the above mentioned methods (on a PIII 870MHz notebook with 128MBytes of RAM). We also illustrate our approach using 5 sets of 8 samples equally spaced across the manufacturing web (i.e. 15cm2 each, spaced 50cm from each other). Three sample sets represent standard operating conditions, and two sample sets represent deviations from the standard operating conditions (i.e. were obtained before stopping production for maintenance). The central part of these samples were scanned using a scanner with transparency unit, at 600dpi, obtaining images with 1000 × 1100 pixels. Three consecutive dyadic scales were used (2j , for j = 1, 2, 3) in the texture analysis. These experiments were conducted to evaluate web uniformity and to rank samples based on their anisotropy, as illustrated in Figures 2(c) and (d). Our results correlate with the tensile test (ratio maximum/minimum tensile) in anisotropy measurements across the web (i.e. coefficient of correlation = 0.9059), and in sample ranking based on their anisotropy (i.e. coefficient of correlation = 0.8844). We also verified that our texture analysis results correlate less with TSI/TSO measurements across the web (i.e. coefficient of correlation = 0.7825). However, it should be reported that we found variability in our results, which we attribute to several and 6 orientations were used).
10
factors, such as: noisy laboratory mechanical measurements, light dispersion introduced by optical sample scanning, sample mass density and thickness variability. It shall be noticed that TSI/TSO focuses on sample sound propagation properties, which may not correlated well the optical properties of the sample. Therefore, it was necessary to take averages of the measurements to reduce their variability. In order to quantify fluctuations and irregularities in systems operation, it is often necessary to sample and test at several locations across the machine (i.e. web), and estimate the web homogeneity [14]. Given a set S of k samples across the web, the sum of pairwise distances between each sample s (² S) and the remaining k − 1 samples in S, namely Djacc (s)), can be used as an web homogeneity indicator: Djacc (s) =
k X J X
j
j
(1 − α)DGauss2 (s, p) + αDRayl (s, p)
(10)
p=1 j=1
An homogeneous web is characterized by small mean and small variability of across the web. Figures 2(e) and (f) show that we obtain different Djacc (s) profiles for different operating conditions. Under regular operating conditions, we obtain smaller Djacc (s) values as well as smaller variability across the web - even at web extremities, where higher variability is expected - (See Figure 2(f)); on the other hand, for non-standard operating conditions, we obtain higher Djacc (s) values and variability, as illustrated in Figure 2(f). It shall be noticed that this web homogeneity test compares stochastic texture images, and only requires that imaging conditions remain the same.
5
Concluding Remarks
In conclusion, we may say that stochastic texture images are acquired in large quantities in continuous industrial processes, and encode important quality 11
and process information. Consequently, methods for objective stochastic texture interpretation and discrimination are important for a large segment of the industry. This work presented a new multi-resolution method for stochastic texture interpretation and discrimination based on the Wavelet Transform. Also, a multi-resolution distance measure for stochastic textures was proposed, and applications of our method in the non-woven textiles industry were reported. The preliminary experimental results obtained by our approach were encouraging, and our texture features tend to correlate well with industrial procedures. However, more research is needed to estimate formation anisotropy using texture analysis methods that are statistically robust, and correlate better with the physical properties of stochastic materials. As future work, we intend to use our approach as a support for data mining operations in stochastic data repositories, with application in preventative maintenance and personnel training.
References
[1] X. Z. Wang, Data Mining and Knowledge Discovery for Process Monitoring and Control, Springer-Verlag, London, 1999. [2] J. Scharcanski, C. T. J. Dodson, Local spatial anisotropy and its variability, IEEE Transactions on Instrumentation and Measurement 49 (5) (2000) 971– 979. [3] G. Fan, X. Xia, Wavelet-based texture analysis and synthesis using hidden markov models., IEEE Transactions on Circuits and Systems - I: Fundamental Theory and Applications 50 (1) (2003) 106–120. [4] S. Arivazhagan, L. Ganesan, Texture classification using wavelet transform., Pattern Recognition Letters 24 (2003) 1513–1521.
12
[5] S. C. Zhu, Y. Wu, D. Mumford, Filters, random fields and maximum entropy., Int. J. Comput. Vision 27 (2) (1998) 1–20. [6] M. N. Do, M. Vetterli, Wavelet-based texture retrieval using generalized gaussian density and kullback-leibler distance, IEEE Transactions on Image Processing 11 (2) (2002) 146–158. [7] H. Choi, R. G. Baraniuk, Multiscale image segmentation using wavelet-domain hidden markov models, IEEE Transactions on Image Processing 10 (2001) 1309– 1321. [8] G. V. Wouver, P. Sheunders, D. V. Dyck, Statistical texture characterization from wavelet representations., IEEE Transactions on Image Processing 8 (1999) 592–598. [9] J. K. Romberg, H. Choi, R. G. Baraniuk, Bayesian tree structured image modeling using wavelet-domain hidden markov models, IEEE Transactions on Image Processing 10 (2001) 1056–1068. [10] S. G. Mallat, S. Zhong, Characterization of signals from multiscale edges, IEEE Transactions on Pattern Analysis and machine Intelligence 14 (7) (1992) 710– 732. [11] D. F. Mix, Random signal processing, Prentice Hall, New York, 1999. [12] H. J. Larson, B. O. Shubert, Probabilistic models in engineering sciences, Vol. I, John Wiley & Sons, New York, 1979. [13] S. Yoshizawa, K. Tanabe, Dual differential geometry associated with the kullback-leibler information on the gaussian distributions and its 2-parameter deformations, SUT Journal of Mathematics 35 (1) (1999) 113–137. [14] G. A. Smook, Handbook for pulp & paper technologists, 2nd. Ed., Angus Wilde Publications, Vancouver, 1994. [15] C. Dodson, W.K.Ng, R.R.Singh, Paper stochastic structure analysis - archive 2 (CD-ROM), University of Toronto, Toronto, Canada, 1995.
13
Classif. Approach
Correct Classif. Rate (%)
Classif. Time (sec.)
Feature Extraction Time (sec.)
GRD, KL
0.8127
0.71
1.37
2G, KL
0.7746
0.66
0.66
Gabor, Euclid.
0.6508
0.33
10.65
SGLDM, Euclid.
0.6476
0.10
15.99
Rayl., KL
0.6476
0.10
15.99
Table 1 Measured correct classification rate and computational costs: GRD, KL is our approach, i.e. bivariate Gaussian and Rayleigh models, and Kullback-Leibler distance (alpha = 0.15); 2G, KL is the approach of bivariate Gaussian model for the wavelet coefficients, and Kullback-Leibler distance; Gabor, Euclid. is the approach of Gabor filters, and Euclidean distance; SGLDM, Euclid. is the approach of SGLDM, and Euclidean distance; and Rayl., KL is the approach of Rayleigh model for the gradient magnitudes, and Kullback-Leibler distance.
14
(a)
(b)
90 10 120
90
60
8
120
8
60 6
6 150
30 150
4
30
4
2
2
180
0 180
210
0
330 210
240
330
300 240
270
300 270
(c)
(d)
4000 5000 3500
4500
4000
3000
3500
3000 Frequency
Frequency
2500
2000
1500
2500
2000
1500
1000
1000 500 500 0
0
5
10
15 20 Gradient Magnitudes
25
30
35
0
(e)
0
5
10
15 20 Gradient Magnitudes
25
30
35
(f)
Fig. 1. Comparative results for isotropic and anisotropic test samples. Stochastic texture images: (a) nearly isotropic and (b) anisotropic; Gaussian ellipses: (c) nearly isotropic and (d) anisotropic; Histograms of local gradient magnitudes (solid-Rayleigh model): (e) nearly isotropic and (f) anisotropic.
15
Normal Probability Plot Normal Probability Plot
0.999 0.997
0.999 0.997
0.99 0.98 0.95
0.99 0.98
0.90
0.95
Probability
0.50 0.25
0.75 0.50 0.25
0.10 0.05
0.10
0.02 0.01
0.05 0.02 0.01
0.003 0.001
0.003 0.001
−30
−20
−10
0 Data
10
20
30 −30
−20
−10
0
10
20
Data
(a)
(b)
Anisotropy Testing: Mechanical x Texture Measurements 1.4 1.25
Anisotropy: texture analysis (solid) , lab. tests (dotted)
Anisotropy (solid−texture, dashed−mechanical)
1.2
1
0.8
0.6
0.4
0.2
0
1
2
3
4
5
6
7
1.2
1.15
1.1
1.05
8
1
Sample
1
2
3
4
5
6
7
8
5
6
7
8
Samples
(c)
(d)
10 10 9 9 8 8 Accumulated Distance to Other Samples
Accumulated Distance to Other Samples
Probability
0.90 0.75
7
6
5
4
3
2
7
6
5
4
3
2 1 1 0
1
2
3
4
5
6
7
8
0
Sample
1
2
3
4 Sample
(e)
(f)
Fig. 2. Normal plots: (a) nearly isotropic and (b) anisotropic. Web uniformity testing results: (c) Sample anisotropy ranking (our approach: solid; tensile test: dotted); (d) Anisotropy measurements
16
across the web (our approach: solid; tensile test: dotted); (e) Distance Djacc profile (standard operating conditions); (f) Distance Djacc profile (non-standard operating conditions).