Image-Based Multiresolution-ANN Approach for Online Particle Size

Article pubs.acs.org/IECR

Image-Based Multiresolution-ANN Approach for Online Particle Size Characterization B. Zhang,† R. Willis,† J. A. Romagnoli,*,† C. Fois,‡ S. Tronci,‡ and R. Baratti‡ †

Department of Chemical Engineering, Louisiana State University, Baton Rouge, Louisiana 70803, United States Dipartimento di Ingegneria Meccanica, Chimica e dei Materiali, Università di Cagliari, Cagliari 09123, Italy

‡

ABSTRACT: In this work an image-based multiresolution sensor for online prediction of crystal size distribution (CSD) is proposed. The mean and standard deviation of log-normal probability density function as the CSD can be predicted through the online sensor. In the proposed approach, texture analysis of fractal dimension (FD) and energy signatures as characteristic parameters to follow the crystal growth is utilized. The methodology consists of a combination of thresholding and wavelettexture algorithms. The thresholding method is used to identify crystal clusters and substrate empty backgrounds. Wavelet-fractal and energy signatures are performed afterward to estimate texture on crystal clusters. Following the texture information extraction, a nonlinear mapping consisting of an artificial neural network (ANN) is incorporated using as inputs the texture information in conjunction with the available online process conditions (flow rate and temperature). A software framework developed in MATLAB enables the configuration of the image acquisition parameters as well as the processing of the online images. Validations against experimental data are presented for the NaCl−water−ethanol anti-solvent crystallization system.

1. INTRODUCTION Crystallization is a critical part of processing in the fine chemical and pharmaceutical industries. During manufacturing, transition substances and/or end-products can exist in crystal particulate form. Crystal size, as an important external crystal structure, can largely influence the textural and physical properties of the final commercial products. For example, the quality of the sugar or the dissolution and bioavailability of a drug with protein crystals depend on the crystal size.1,2 A narrow size distribution with desired mean size is preferred by industries, and therefore it is important to efficiently monitor and control the process for assuring that target end-product size can be obtained. The development of monitoring and feedback control tools requires the availability of real-time measurements of the crystal size. The available devices suited for this issue are commonly based on the technologies of laser diffraction,3 ultrasound attenuation,4 and laser reflectance.5 The main drawbacks of such instruments are high costs, assumptions of specific shapes of crystals, and requirements of crystal dispersion prior to size measurement. On the other hand, the recent progress in high speed imaging devices, powerful computers at reasonable costs, and the adaptability to real-time applications have pushed the researcher toward the development of image-analysis-based methods for the direct measure of crystal size, size distribution, and shape characterization. In the crystallization field, Argaw et al. proposed an approach consisting of foreground and background markers, ultimate erosion, distance transform, and watershed segmentation to automatically measure the crystal size distribution using images.2 The accuracy of the results, according to the paper, depends on the number of crystals counted and the extent of crystal overlapping. Wavelet decomposition of crystal images coupled with artificial neural network modeling had been applied to estimate size measurements during sugar crystallization.6 This proposed approach © 2014 American Chemical Society

used the neural model to approximate the relationship between coefficients of wavelet decomposition and mean grain size, obtaining good estimation when clear images were available, but failing in the case of overlapping, noise, and nonuniform illumination. Papers by Larsen and Rawlings7,8 address the development of algorithms using artificial images to measure the size distribution in conditions of low, moderate, and high solid concentrations. They showed that these algorithms work well in low and moderate concentrations but do not perform well at high solid concentrations. The possibility of applying the image-based technique to industrial practice is attractive because of the acceptable cost of such an approach, but it is mandatory to ensure reliability of the measurements. A robust measurement method suitable for industrial application should be able to extract relevant information from real images even when they are not clear because of illumination issues, in the presence of noise, and in the case of agglomeration and overlapping phenomena. The previously proposed methods2,6−8 will perform well in measuring size or size distribution if every individual crystal in the image can be identified, but when dealing with real processes that assumption does not always hold. Separating overlapping crystals or avoiding agglomeration can be achieved by physical dispersion prior to capturing the images; however, this operation cannot completely eliminate these phenomena, because dispersed crystals will flow together, especially in high solid concentration conditions. An algorithm based on a multiresolution fuzzy clustering approach was proposed by Zhang et al.,9 with the purpose of separating touching and overlapping regions. The authors performed texture analysis on Received: Revised: Accepted: Published: 7008

June 17, 2013 March 25, 2014 April 3, 2014 April 3, 2014 dx.doi.org/10.1021/ie4019098 | Ind. Eng. Chem. Res. 2014, 53, 7008−7018

Industrial & Engineering Chemistry Research

Article

If only one feature is of interest in an image, separating the image into two classes, a foreground and a background, by applying a single threshold value is sufficient during the thresholding process. The foreground is those pixels with intensity values higher than the threshold value, and the upper boundary for thresholding implementation is the highest pixel value in the image; otherwise is the background. In the thresholding algorithm, the threshold value is the significant parameter for successfully dividing segments, and thus appropriately selecting threshold value becomes the major concern. Sezgin and Sankur24 conducted a survey of methods of thresholding value selection. They categorized the methods according to the information they are exploiting into histogram shape-based methods, clustering-based methods, entropy-based methods, object attribute-based methods, spatial methods, and local methods. They concluded that the use of a certain algorithm depends on the characteristics of systems. In our case, attention was paid to the complexity of the crystal edge appearance and the degree of light illumination uniformity. Due to the unique structure of crystals, some edges can be easily recognized (with distinct different intensities compared with their neighborhoods) while others (slight edges) can be hardly identified (with similar intensity levels around) since crystal flats can face with different orientations. Taking these two characteristics into consideration, in our application methods that work directly on the variance of intensity itself are better, as well as a correction factor to reduce the influence of nonuniform illumination if necessary for specific image setup systems. Otsu’s global thresholding has been widely used and proved to be one of the most successful techniques in image thresholding to find the optimal threshold value when only two texture classes in an image are to be identified, and this method is utilized in our approach.25 The algorithm is performed by comparing the variance within class as the summation of weighted variance in each class (eq 2) with the variance between classes shown in eq 3. The optimum value is the one that can generate the lowest within-class variance and highest between-classes variance.

crystal clusters aiming at accurate crystal image processing and subsequently correct control action. The other disturbances that can influence extracting relevant information from images are noise, in spatial and frequency domain and in high and low frequency domain. Nonuniform illumination, as low frequency noise, can happen during manufacturing processes due to aging filaments, faulty reference voltages, or contaminated apertures.10 Techniques have been developed to remove those noises according to their existing domains, such as linear/ nonlinear spatial filtering reducing spatial noise and Fourier transform deleting high frequency noise.11 Although the nonuniform illumination problem can be solved by mechanical operations, its occurrence can necessitate additional diagnostic effort for the monitoring and control systems. Because also the repair time can delay manufacturing, an algorithm that can deal with or consider noise in all of the above-mentioned domains without identifying the noise type will be attractive and welcomed for broad application in industry. In this paper, an image-based multiresolution sensor for online prediction of crystal size distribution (CSD) is proposed, where the mean and standard deviation of log-normal probability density function as the CSD can be predicted through the online sensor. A novel approach is proposed for crystal image texture analysis based on combing thresholding and wavelet-fractal-energy algorithms. The thresholding method can quickly detect crystal clusters and remove the background on the basis of the intensity value and is a good choice for real-time application. Fractal dimension (FD) as the texture analysis parameter,12−15 estimated through a waveletfractal approach,13,16,17 is a useful and quantitative analytical parameter to characterize many kinds of complicated selfsimilar substances in nature. Wavelet power spectra have been widely adopted for analysis in chemical fields.9,18−22 FD provides a noninteger value to describe crystal growth from crystal images from the point of view of surface roughness. Generally, the higher the value of the fractal dimension, the more rough the surface is, while the preprocessing step of subtracting backgrounds eliminates their contribution to surface roughness. Neural networks are then used to identify the inputoutput relationship between FD and crystal size distribution, which is the main issue in the crystallization control framework. In particular the neural model predicts mean and standard deviation for a given image from the calculated FD along with the energy signature obtained from wavelets and process conditions. Validations against experimental data are presented for the NaCl−water−ethanol anti-solvent crystallization system.

σ W2(t ) = wb(t )σb2(t ) + wf (t )σf2(t )

where w and σ2 are the weight and variance. The subscripts W, b, and f mean within class, background, and foreground. t represents a threshold value. σB2(t ) = σ 2(t ) − σ W2(t )

n = 1, 2, ..., N

(3)

The subscript B means between classes. Each image in our approach is preprocessed to retain objects and remove backgrounds by segmentation. The threshold method is applied here to identify the crystal edges generating a binary image (crystal cluster binary image) that shapes the crystal edges and provides the location information on crystal clusters. This information could be extracted by setting the x-y axes coordinates on the crystal cluster binary image. The location for each crystal cluster, in the crystal cluster binary image, is the same as in the original one since the size of the former image is the same as the latter. Subtracting backgrounds from the original image allows us to identify crystal clusters through their locations and store them in an empty image crystal (crystal clusters image). The intensities of the crystal clusters are then rearranged in a sequence of the same column. The rearrangement transfers the

2. IMAGE SEGMENTATION When an image contains different texture objects, thresholding can basically identify them. It is a segmentation approach operating in a way of comparing objects’ intensity values with a certain threshold value.23 Given an upper or lower boundary value of a specific segment, pixels whose intensity values are less or greater than the threshold value are recognized as the elements of that segment. Supposing an image containing N segments, the mathematic model of image segmentation by the thresholding method can be established as the following set: Gn = {Tn − 1 ≤ f (x , y) ≤ Tn}

(2)

(1)

where Gn represents nth segment and Tn,Tn−1 are its upper and lower boundaries, respectively. f(x,y) is the value of the pixel located in the xth row and yth column. 7009

dx.doi.org/10.1021/ie4019098 | Ind. Eng. Chem. Res. 2014, 53, 7008−7018


Article +∞

2-D image problem into a 1-D vector, containing less but meaningful data. The intensity vector is then used as an input data for the wavelet transformation.

Sm , n =

1 a0m

0

+∞

∑ ∑

Tm , nψm , n(x)

m =−∞ n =−∞

(7)

3.2. Wavelet Energy Signature. The basic idea of wavelet texture analysis is to extract texture features using detail coefficients. High frequency components, represented as detail coefficients, carry the characteristic information on the signal such as edges of objects distinguishing its correspondence signal from others.20 Low frequency components, approximation coefficients at higher decomposition level, on the other hand, usually reflect the foregrounds’ outline providing no sizerelated information. Wavelet energy signature,29 the square root of the variances of detail coefficients at several scales, is adopted as texture features under the assumption that each texture display its unique energy distribution at all scales. After wavelet decomposition, the energy Em stored in the detail coefficients at scale index m is calculated through the following equation: 2M − m − 1

Em =

∑

(Tm , n)2

n=0

(8)

Dividing by 2M−m on both sides, we get Em 2M − m

2M − m − 1

=

∑n = 0

(Tm , n)2

2M − m

(9)

where the left-hand is the coefficient variance at scale m labeled ⟨Tm,n2⟩m· According to this expression, the wavelet coefficient variance can be explained as the average energy wrapped up per coefficient at each scale m. In this context, the square root form of the coefficient variance is used. 3.3. Fractal-Wavelet Feature. Fractal dimension (FD) is a statistical real number that measures how complicated a fractal is. Fractals are objects that possess self-similar property, where fractals appear identical at different scales or numerical or statistical measures of fractals are consistent across scales. The reason why FD is selected as the texture feature is that the dynamics of crystal growth follows a fractal process14 and crystals in images exhibit statistical self-similarity in reality. Several approaches are available for calculating FD and in our context we use a fractal-wavelet method that requires not only wavelet decomposition but also a fractional Brownian motion (fBm) model. A fBm is a continuous-time Gaussian process defined as a nonstationary and zero-mean Gaussian random function shown in the following equation:30

ψ (a0−mx − nb0)

f (x)ψm , n(x) dx

Sm0 , nϕm , n(x) +

n =−∞

(4)

+∞

∫−∞

(6)

m0

+∞

∑

f (x ) =

where a0 is a fixed dilation parameter and is greater than 1, and b0 is the location parameter and is greater than zero. The discrete wavelet transformation of function f(x) is thus a function of m and n: Tm , n =

f (x)ϕm , n(x) dx

Taking multiresolution into consideration, we can write the function f(x) as a sum of approximation of the function at arbitrary scale index m0 and detail function from scale m0 to −∞, using the approximation coefficients and detail coefficients:

3. WAVELET TEXTURE ANALYSIS On the examination of crystal images, it can be found that crystal size information reflects on the intensity variation in a given small neighborhood: images for small crystals contain more variations compared with those for large crystals. Capturing crystal size information equals describing these spatial variations in pixel intensities. A good descriptor of this characteristic is the image texture, according to its loose definition by Russ.26 Implementing texture analysis can extract meaningful information and present them as textural features. Among the diverse texture analysis approaches such as statistical, structural, model-based, and transform-based methods,27 the wavelet transform method is used in our approach not only to extract textural characteristics but also as a tool for noise removal. 3.1. Wavelet Transformation. Wavelet transform technique converts a signal under investigation into another form that presents the signal information in a more succinct and/or useful representation. The transformation does not change the signal information contents but decomposes them into approximation and detail. Not only can the original signal be decomposed, but also the obtained approximation coefficients can acquire decomposition in the next resolution, thereby forming a multiresolution decomposition. Theoretically, the maximum wavelet decomposition scale is M, which should satisfy the constraint of 2M ≤ N < 2M+1, where N is the total sample data of the signal. The transformation can be done in continuous or discrete manner, termed continuous wavelet transformation (CWT) and discrete wavelet transformation (DWT), respectively. CWT requires calculating the wavelet coefficients at every possible scale and generates plenty of data, which leads to arduous work and decreases analysis efficiency. DWT, on the other hand, as a more efficient yet still accurate method, calculates the coefficients at a selected scale a and location b parameter.28 They are discretized in such a manner that a and b are linked. The scale a is generally discretized in a logarithmic way, am0 where m is an integer. Each location b can be reached in discrete steps n (an integer) from an origin. It is also proportional to the scale am0 . Thus, it can be represented as nb0am0 . In wavelet transformation, a wavelike function called wavelet is used to transform a signal or function f(x) in space or time into another form by convolution. The wavelet function ψm,n in the discrete form is ψm , n(x) =

∫−∞

(5)

The discrete wavelet transformation Tm,n is known as the detail coefficients. The approximation coefficients Sm,n, like the detail coefficients, are given by convoluting the function f(x) and another set of functions called scaling function, ϕm,n: 7010



Article

BH (0) = 0 1 Γ(H + 0.5)

BH (t ) − BH (s) =

0

∫−∞ [(t − s)H−0.5 − s H−0.5] dBH (s) ∞ + ∫ (t − s)H − 0.5 dBH (s)} 0 {

(10)

where BH represents the Gaussian process with observation time at 0, t and s. Γ is the Gamma function. The Hurst exponent, H, has a value between 0 and 1 characterizing a fBm. It is an indicator to reflect the smoothness of the fBm function: the higher the Hurst exponent is, the smoother the fBm function is. Hurst exponent is directly related to FD via the relationship FD = DE + 1 − H

Figure 1. A topology structure of a three-layer feedforward ANN. n2 + 1

z 3(i) = fe ( ∑ w2(j , i)z 2(j)), with z 2(n2 + 1) = 1

(11)

j=1

Here DE is the Euclidean dimension. For a fBm, it also scales to wavelet power spectra Pw(f m) and frequency f m such as28 Pw(fm ) ∝ fm−(2H + 1)

The proposed neural network considers a bias term in the input and hidden layers, which have been indicated with z1(n1 + 1) and z2(n2 + 1) in eqs 16 and 17. For the case at hand, two ANNs are developed using the Neural Network Toolbox of Matlab, to predict the mean and variance of the crystal size distribution, using the fractal and the wavelet energy signatures along with process conditions as model inputs. Nonlinear transfer functions between layers allow the network capturing the nonlinear relationships between input and output. In this work, a tan-sigmoidal and linear activation functions are used, respectively, for input eq 18 and hidden neurons eq 19:

(12)

Meanwhile, wavelet power spectra Pw(f m) is related to the variance of discrete wavelet coefficients as Pw(fm ) ∝ T ⟨m2 , n ⟩m

(13)

Combining eq 12 with 13, and the fact that frequency f m is inversely proportional to the wavelet scale am, the scaling relationship becomes

⟨Tm2 , n⟩m ∝ am(2H + 1)

(14)

f0 =

Computing the Hurst exponent (and consequently FD) can then be achieved by taking the base 2 logarithms of both sides of eq 14, leading to eq 15. The logarithmic plot of the variance of the discrete wavelet coefficients provides a method for calculating the Hurst exponent from the slope of the plot: log 2(⟨Tm2 , n⟩m )

= (2H + 1)m + constant

n1+ 1 w1(j , i)z1(j))

1 + e−2(∑ j=1

−1 (18)

∑ w2(j , i)z 2(j) (19)

j=1

The estimation of the network parameters (weights) is generally referred to as training of the network and has been developed by means of the Levenberg−Marquardt algorithm. The objective function to be minimized is the mean square error:

(15)

4. ARTIFICIAL NEURAL NETWORK The information extracted from the crystal images through the wavelet texture analysis needs to be related to CSD mean and variance, in order to obtain a monitoring tool suited for online applications. In this work, a feedforward artificial neural network (ANN) with three layers has been used to address the underlying identification problem. The general ANN scheme is shown in Figure 1, for n1 inputs, represented by the vector z1 and one output, z3.31 In more details, the input signals z1 are multiplied by the adjustable parameters, called weights, w1(i,j), and then all the contributions are summed and processed by the activations function f 0 (eq 16). The resulting signal, z2, is multiplied by the weight w2(i,j), and its components are summed and mapped into the vector z3 by the activation function fe (eq 17):

n

MSE =

∑ i=1

(y − t )2 n

(20)

where n is the number of data (targets) used for the training, y indicates the mean (or variance) predicted by the network, that is, y = z3(1) by considering the scheme in Figure 1, and t the experimental (target) value. The training of the network terminates when the neural model has learnt the task, and therefore some criteria are necessary for measuring this condition. One stop criteria can be related to MSE, by setting an acceptable error or applying a threshold to the incremental change in MSE. However, this criterion does not allow obtaining a neural network with the best generalization. This task can be accomplished by dividing training data into two sets: training and validation sets. The error calculated on the validation set, which is not used for training, decreases with the number of training iterations but then starts to increase: this means that the neural model is becoming too specialized on the

n1+ 1 j=1

2

n2 + 1

fe =

where the constant depends both on the wavelet function and the Hurst exponent.

z 2(i) = f0 ( ∑ w1(j , i)z1(j)), with z1(n1 + 1) = 1

(17)

(16) 7011



Article

Figure 2. Overall architecture of CSD prediction.

Figure 3. Curves for changed flow rate (a) and temperature (b) versus crystallization time.

training set (overtraining),31 and therefore test set performance starts deteriorating. The number of iterations leading to the minimum validation errors should allow the best generalization.

images are then stored in a database in the computer and put through image analysis. The sequence of steps for implementing the proposed methodology for CSD prediction from crystal images are shown in rectangular boxes in Figure 2. An input image is treated as a 2D array of pixel intensities. A thresholding algorithm is applied for extracting features of interest, which are the crystal clusters in this work. This is accomplished through a series of three substeps: (a) detection of crystal edges by a threshold value differentiating them from the background, (b) detection of the locations of the crystal clusters with the help of x-y coordinates on the binary image of crystal clusters, and (c) determining and extracting the intensity values belonging to the crystal clusters. The information data of the clusters is then restored into a vector, being processed to

5. METHODOLOGY OF THE IMAGE-BASED MULTIRESOLUTION SENSOR The overall architecture for online crystal size distribution prediction by the image-based multiresolution sensor is given in Figure 2. The crystallization is carried out in a reactor with suitable data acquisition and control system for adjusting both temperature and anti-solvent addition. Crystal suspension is circulated through a flow cell where the crystal images are taken by a camera connecting to a microscope and a computer. The 7012



Article

generate the texture features by means of wavelet-fractal-energy algorithm. In this regard, it is first decomposed by wavelet transformation at several levels into details and an approximation. The detail from lower decomposition level and the approximation are considered as the high and low frequency noise to be removed. The remaining details are then used for finding their variance and the 2-based logarithm of the variance. The 2-based logarithm of the variance at each scale and decomposition scale are then used to calculate the Hurst exponent and consequently the FD. The texture features as well as crystallization conditions (temperature and feed flow rate) are then used to train an ANN (during the training process) for the prediction of crystal mean size or standard deviation. When a new set of images arrive (during online application), they go through the same texture analysis procedures and generate a set of input variables to be used in conjunction with the trained ANN to predict the CSD.

Table 1. Database for Different Crystallization Conditions and Stages condition temp (°C)

6. CASE STUDY 6.1. Experiment. A set of anti-solvent crystallization experiments at both constant and varying (temperature and feed flow) conditions was carried out with anti-solvent flow rate and temperature as the process parameters. For constant conditions, three different values were chosen: 0.7, 1.5, and 3 mL/min as flow rate and 10, 20, and 30 °C as temperature (training data set). For the varying conditions (test data set), the temperature and flow rate were changed according to the profiles given in Figure 3, respectively. At the startup conditions, the crystallizer was filled with the saline solution composed of 100 g of water and 34 g of NaCl. The solvent solution was made up of 95% ethanol and was added by the time to the initial solution using a peristaltic pump. The liter batch reactor used was continuously stirred and immersed in a bath, and the solution temperature was controlled through a resistance temperature sensor connected to a cooling/heating control system. As for the crystal growth examination, this was performed choosing to work in an online fashion by using a peristaltic pump to direct the solution into a glass cell. A laboratory scale software/hardware framework for capturing crystal images for this case study was set up. The experimental setting utilizes a USB microscope camera (model MD900) with a resolution of 1280 × 960 pixels, which fits into the side tube of the microscope with one of the supplied adapters and connects to a computer. The AMSCOPE software is utilized to capture images and for manual measurement (individual particle analysis). The magnification used is 25x, which corresponds to 0.775 μm/pixel. This conversion factor was used to manually measure individual crystal sizes on each image. Images were taken at different crystallization stages (listed in Table 1) for each experimental condition. At each crystal growth stage, a set of at least 10 images capturing different amounts of crystals were utilized. 6.2. Texture Analysis on Crystal Images. A sample image from crystallization at 10 °C and 0.7 mL/min at 30 min is used to illustrate the method of FD estimation from crystal images. The original gray crystal image in TIF format was imported for image analysis. The original gray crystal image was compressed and cropped with Microsoft Office Picture Manager into a size of 960 × 1280 in JPG format as shown in Figure 4a to improve execution time for further analysis. Then the cropped image was imported into Matlab for image analysis, wavelet analysis, and ANN modeling. To make calculations simpler, the image intensities have been normalized

anti-solvent flow rate (mL/min)

10

0.7

10

1.5

10

3.0

20

0.7

20

1.5

20

3.0

30

0.7

30

1.5

30

3.0

varying

varying

crystallization stage (min) 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 10, 15, 20, 30, 60, 90, 120, 180, 240, 300 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 10, 15, 20, 30, 60, 90, 120, 180, 240, 300 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 30, 60, 90, 120, 180, 240, 300, 360, 420, 480 10, 15, 20, 30, 60, 90, 120, 180, 240, 300 20, 30, 60, 90, 120, 180, 240, 300, 360, 420, 480

with intensity value in a range of 0 to 1. The normalized image then underwent thresholding and morphological operations to detect crystal clusters. Otsu’s method multiplied by a parameter of 1.1 (correction factor to reduce the influence of nonuniform illumination) generated a threshold value of 0.7182 for the normalized image, implying that pixels with intensity value below 0.7182 belong to edges. The threshold method produced a binary image as shown in Figure 4b. The holes presented in the edge detection binary image correspond to crystal bodies in the original image. A morphological operation was used to fill the holes and generated the image shown in Figure 4c (crystal cluster binary image), which provides the location information on crystal clusters. The locations of each crystal cluster in the crystal cluster binary image are the same as those in the original image since the size of the former image is the same as the latter. Identifying crystal clusters and subtracting backgrounds was achieved through determining the locations of crystal clusters and storing them in an empty image. The detected crystal clusters are depicted in Figure 4d. The intensities of the crystal clusters were then rearranged in a sequence of the same column order connecting with the next row as a vector with a length of 224203. The rearrangement transfers the 2-D image problem into a 1-D problem. The intensity vector was the input data for the wavelet-fractal analysis. The wavelet function and decomposition level are the two parameters for wavelet transformation. In this study, wavelet function ‘db3′ was chosen, and decompositions at six levels were performed. The transformation results can be written in the form W = {cd1, cd 2, ..., cd6, ca6}

(23)

where W represents the signal data, and cd and ca are details and approximation, respectively. The subscript number indicates the decomposition level. We select details cd3−cd6 for further analysis because cd1 and cd2 contain the high frequency noise and ca6 captures the low frequency noise as nonuniform illumination. The variance of the detail coefficients and their 2-based logarithm were calculated as listed in Table 2. 7013



Article

Figure 4. (a) Preprocessed, (b) edge detection binary image, (c) crystal clusters binary image, and (d) crystal clusters image.

flow-rate and temperature, those variables have also been selected as inputs. Furthermore, the time dependence (sampling time) of crystal growth can be used to give further information to the model and improve its prediction capability. This seems to be the case when blurred or unclear images are present, which are possible when dealing with practical applications. As the inputs have been selected, it is necessary to set the number of hidden neurons that leads to the best prediction capability of the ANN. This parameter cannot be determined from the knowledge of the process, because hidden units elaborate signals that have lost the physical meaning of the inputs. The choice of the number of the hidden neurons has been based on the results obtained during the learning phase with different network structures. In particular, the number of hidden units has been selected considering the most parsimonious model showing the best generalization capability (low training and validation error). The available experimental data used as targets for ANN weight estimation are the manually measured mean size and CSD standard deviation (obtained from the best global model using nonlinear Fokker−Planck Equation32). They have been divided into two sets: (i) training set used for estimating the model weight and (ii) test set used to assess the prediction capability of the developed neural network (unseen data). Part of the training data set (unused images) has been devoted to model validation (30%), in order to ensure best generalization of the network. Images from the constant conditions comprise the training/validation set for the ANN model, while those obtained from the varied conditions are used for testing. Different images have been taken at the same stage from a specific condition, and all of them are used in the learning phase. To improve training efficiency the inputs and outputs data have been normalized through the following function:

Table 2. Detail Variance and Its log2 at Four Decomposition Levels decomposition level detail variance log2 of detail variance

3

4

5

6

0.0112 −6.480 376

0.082 725 −3.595 54

0.256 18 −1.964 77

0.489 008 −1.032 07

Performing a least-squares linear fit for 2-based logarithm of the detail variance and its corresponding decomposition level, as shown in Figure 5, we found the slope to be 1.7976. Using eqs

Figure 5. log2(detail variance) against decomposition level plots with the linear fitted line. The Hurst coefficient can be obtained through the fitted equation.

15 and the previously calculated slope, the Hurst exponent and FD were evaluated to be 0.3988 and 1.6012, respectively, when 1 was adopted as the Euclidean dimension for the 1-D situation. 6.3. ANN Models and Prediction. Two three-layer neural models were developed for estimating mean crystal size and variance, meaning that each ANN has one output. The selection of the ANN inputs and number of hidden neurons is not straightforward and necessitates a closer analysis of the problem being studied. Inputs to the network have been selected on the basis of the available information extracted by the wavelet texture analysis and the knowledge on the crystallization process. The first input of the network is the fractal dimension, it being related to the regularity of the images. Further information on CSD can be obtained from the wavelet energy signature, and it has been found that decomposition level from 3 to 6 contains relevant information on the crystal distribution and improves the network prediction capability. The crystal growth being dependent on anti-solvent

xn = 2

x − xmin −1 xmax − xmin

(24)

where x is a specific element of the original input/target data, xmin and xmax represent the minimum and maximum values of each row of the input/target matrix, and xn is the corresponding normalized value. The learning task of the ANN starts by randomly initializing the W vector value and it stops when the following conditions are satisfied: (i) the gradient of the objective function eq 20 reach a value of 1.0 × 10−5, and (ii) the error on validation set does not improve for six successive iterations. The prediction capability of the developed neural model has been assessed by comparison with experimental data, considering the following performance index: root-mean-square error (RMSE) and the error standard deviation (STDV): 7014



Article

Figure 6. Comparison of manually measured mean size and predicted value for test crystallization run (a) and expected standard deviation and predicted value (b) for the ANNs without sampling time as input.

RMSE =

∑ (e)2 n

STDV =

∑ (e − e ̅ )2 n

Table 4. Statistical Parameters of ANN-for-STD without Sampling Time for STD Prediction for Both Training and Testing Sets

(25)

predicted STD by ANN

(26) 10 °C, 0.7 mL/min 10 °C, 1.5 mL/min 10 °C, 3.0 mL/min 20 °C, 0.7 mL/min 20 °C, 1.5 mL/min 20 °C, 3.0 mL/min 30 °C, 0.7 mL/min 30 °C, 1.5 mL/min 30 °C, 3.0 mL/min varying conditions

where e indicates the difference between experimental and predicted values. Figure 6a represents the comparison between manually measured mean crystal size and the predicted one obtained with a network with 6 inputs (FD, energy signatures for details level from 3 to 6, temperature, and anti-solvent flow rate) and 8 hidden neurons for test data set (obtained under varying conditions). The RMSE and STDV performance indexes are summarized in Table 3 for both training and test data sets. Table 3. Statistical Parameters of ANN-for-Size without Sampling Time for Size Prediction for Both Training and Testing Data Sets RMSE

STDV

4.18 9.28 6.66 4.62 8.76 6.23 7.5 11.3 10.7 5.37

4.31 9.40 6.24 4.80 9.17 5.85 7.76 11.87 10.06 5.54

STDV

2.78 2.75 3.78 1.45 2.93 1.23 1.32 0.89 3.80 2.50

2.81 2.79 4.00 1.52 3.07 1.30 1.47 0.94 4.02 2.63

data set, as shown by the results, is very good. However, the performance may be affected in some instances due to the presence of blurred images, in the training set, as detected during our investigations. Another test was performed by augmenting the input vector of the ANN and adding a variable that takes into account the batch time (sampling time). The results obtained with the new structure of the two neural models (7 inputs and 8 hidden neurons for mean size and 10 for STD) are reported in terms of performance indexes in Table 5 (mean crystal size) and Table 6 (standard deviation) and Figure 7. We can appreciate that information related to the duration of the crystallization increases the robustness of the strategy against image quality deterioration and effectively improves the prediction capabilities (slightly smaller values of RMSE and STDV). 6.4. CSD Prediction and Comparison. The capabilities of the developed monitoring system have also been verified in terms of the crystal size distribution. First, the log-normal distribution parameters μ and σ are calculated from the mean m and the standard deviation s of a normal distribution (eqs 27 and 28), considering the manually measured mean crystal size and standard deviation of counted crystals and the predicted ones:

predicted size by ANN 10 °C, 0.7 mL/min 10 °C, 1.5 mL/min 10 °C, 3.0 mL/min 20 °C, 0.7 mL/min 20 °C, 1.5 mL/min 20 °C, 3.0 mL/min 30 °C, 0.7 mL/min 30 °C, 1.5 mL/min 30 °C, 3.0 mL/min varying conditions

RMSE

Figure 6b, on the other hand, represents the results obtained with the second ANN for the prediction of CSD variance. The ANN for the second moment prediction has, again, 6 inputs (FD, energy signatures for details level from 3 to 6, temperature, and anti-solvent flow rate) and 20 hidden neurons, and the corresponding performance index are summarized in Table 4 for both training and test data set. In general the prediction capability of the two ANNs for the test 7015



Article

The calculated variables are reported in Table 7, considering the prediction obtained by the two ANN structures described in the previous section.

Table 5. Statistical Parameters of ANN-for-Size with Sampling Time for Size Prediction for Both Training and Testing Sets predicted size by ANN 10 °C, 0.7 mL/min 10 °C, l.5 mL/min 10 °C, 3.0 mL/min 20 °C, 0.7 mL/min 20 °C, l.5 mL/min 20 °C, 3.0 mL/min 30 °C, 0.7 mL/min 30 °C, 1.5 mL/min 30 °C, 3.0 mL/min varying conditions

RMSE

STDV

1.29 0.81 1.52 2.47 1.62 0.86 1.78 1.03 0.87 3.57

1.28 0.82 1.61 2.60 1.70 0.86 1.88 1.08 0.91 3.31

Table 7. Manually Measured and Predicted Log-normal Distribution Parameters at Each Sampling Time during the Batch Run (Varying Conditions) manually measured

Table 6. Statistical Parameters of ANN-for-STD with Sampling Time for STD Prediction for Both Training and Testing Sets predicted STD by ANN 10 °C, 0.7 mL/min 10 °C, 1.5 mL/min 10 °C, 3.0 mL/min 20 °C, 0.7 mL/min 20 °C, 1.5 mL/min 20 °C, 3.0 mL/min 30 °C, 0.7 mL/min 30 °C, l.5 mL/min 30 °C, 3.0 mL/min varying conditions

μ = ln(m2 / s 2 + m2)

σ=

⎞ ⎛ s2 ln⎜ 2 + 1⎟ ⎠ ⎝m

RMSE

STDV

0.37 1.33 0.57 0.27 1.18 0.87 1.16 0.60 0.63 3.28

0.39 1.36 0.59 0.29 1.24 0.93 1.20 0.63 0.67 2.17

predicted (with sampling time)

predicted (without sampling time)

time (min)

μ

σ

μ

σ

μ

σ

20 30 60 90 120 180 240 300 360 420 480

4.45 4.43 4.59 4.69 4.72 4.93 5.02 5.06 5.05 5.06 5.06

0.43 0.55 0.471 0.47 0.47 0.40 0.37 0.35 0.36 0.34 0.34

4.47 4.51 4.57 4.62 4.77 4.90 4.99 5.03 5.04 5.05 5.05

0.39 0.41 0.45 0.44 0.40 0.40 0.38 0.38 0.38 0.38 0.39

4.50 4.55 4.57 4.67 4.77 5.00 5.05 5.05 5.02 5.02 5.00

0.44 0.43 0.43 0.42 0.44 0.38 0.39 0.39 0.40 0.40 0.41

The CSD model prediction obtained with the mean and standard deviations as predicted by the proposed sensor, in the case of varied conditions (test data set), are reported in Figure 8 along with the raw histogram of the experimental data and the estimated probability density function of the raw data. Different samples obtained at 20, 30, 60, 90, 120, 180, 240, 300, 360, 420, and 480 min during the crystallization experiment were collected for a comprehensive representation. It is shown, through the comparison between experimental crystal size distribution (obtained from raw data) and the estimated one, that the proposed sensor predictions are again very good and the differences are within the experimental error for both with and without inclusion of sampling time as input variable (performance is slightly better with the inclusion of the sampling time).

(27)

(28)

Figure 7. Comparison of manually measured mean size and predicted value for test crystallization run (a) and expected standard deviation and predicted value (b) for the ANNs with sampling time as input. 7016



Article

Figure 8. Comparison among predicted CSD, raw histogram, and the smoothed approximation with the manually measured log-normal distribution parameters at each sampling time from the crystallization batch run with varying conditions. (The legends are the same as those at 20 min.)

with the results that higher FD can be obtain at higher flow rate for a given temperature and at lower temperature for a given flow rate. The pattern of FD and crystal mean size for different crystallization conditions were compared, showing their similarity. The relationship of FD and crystal mean size was extracted and built as an ANN model for predicting crystal mean size as well as for standard deviation. The predicted CSD in log-normal probability density function was plotted and compared with experimental data. These results attest to the

7. CONCLUSIONS An image-based approach of texture analysis combining thresholding and wavelet-fractal with ANN for the prediction of CSD was proposed and applied on a case study. This method could successfully and automatically identify crystal clusters and estimate the texture by means of FD and wavelet energy signature. The FD transformation tendency for different crystallization operation conditions had been investigated 7017



Article

(19) Jetter, K.; Depczynski, U.; Molt, K.; Niemoller, A. Principles and applications of wavelet transformation to chemometrics. Anal. Chem. Acta 2000, 420 (2), 169−180. (20) Reis, M. S.; Bauer, A. Wavelet texture analysis of on-line acquired images for paper formation assessment and monitoring. Chemom. Intell. Lab. Syst. 2009, 95 (2), 129−137. (21) Facco, P.; Masiero, A.; Bezzo, F.; Barolo, M.; Beghi, A. Improved multivariate image analysis for product quality monitoring. Chemom. Intell. Lab. Syst. 2011, 109 (1), 42−50. (22) Sergey, K. Extracting useful information from images. Chemom. Intell. Lab. Syst. 2011, 108 (1), 2−12. (23) Shapiro, L. G.; Stockman, G. C. Computer Vision, 1st ed.; Prentice Hall: Upper Saddle River, 2001. (24) Sezgin, M.; Sankur, B. Survery over image thresholding techniques and quantitative performance evaluation. J. Electron. Imaging 2004, 13 (1), 146−165. (25) Otsu, N. Threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern. 1979, 9 (1), 62−66. (26) Russ, J. C. The Image Processing Handbook. 5th ed.; Taylor & Francis, Inc.: Abington, 2006. (27) Bharati, M. H.; Liu, J. J.; MacGregor, J. F. Image texture analysis: methods and comparisons. Chemom. Intellig. Lab. Syst. 2004, 72 (1), 57−71. (28) Addison, P. S. The Illustrated Wavelet Transform Handbook; Institute of Physics Publishing Bristol and Philadelphia, Napier University: Edinbrugh, U.K., 2002. (29) Barner, K. E.; Barner, B. E.; Arce, G. R. Nonlinear Signal and Image Processing: Theory, Methods, and Applications; CRC Press: Boca Raton, 2003. (30) Biagini, F.; Hu, Y.; Oksendal, B.; Zhang, T. Stochastic Calculus for Fractional Brownian Motion and Applications; Springer-Verlag: New York, 2008. (31) Principe, J. C.; Euliano, N. R.; Lefebvre, W. C. Neural and Adaptive Systems: Fundamental through Simulations; John Wiley & Sons: New York, 1999. (32) Grosso, M.; Cogoni, G.; Baratti, R.; Romagnoli, J. A. Stochastic approach for the prediction of PSD in crystallization processes: Formulation and comparative assessment of different stochastic models. Ind. Eng. Chem. Res. 2011, 50 (4), 2133−2143.

potential application of the proposed method for crystal production process monitoring and control.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (225) 578-1377. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS J.A.R. kindly acknowledges financial support by NSF through the Award no. 1132324 REFERENCES

(1) Wang, X. Z.; Roberts, K. J.; Ma, C. Crystal growth measurement using 2D and 3D imaging and the perspectives for shape control. Chem. Eng. Sci. 2008, 63 (5), 1173−1184. (2) Argaw, G. A.; Alport, M. J.; Malinga, S. B. Automatic measurement of crystal size distribution using image processing. Proc. Annu. Congr. - S. Afr. Sugar Technol. Assoc. 2006, 399−411. (3) Myerson, A. S. Handbook of Industrial Crystallization, 2nd ed.; Butterworth-Heinemann: Oxford, 2002. (4) Raj, B.; Palanichamy, P.; Rajendran, V. Science and Technology of Ultrasonics, 1st ed.; Alpha Science International, Ltd.: Oxford, 2004. (5) Bakeev, K. Process Analytical Technology: Spectroscopic Tools and Implementation Strategies for the Chemical and Pharmaceutical Industries, 2nd ed.; Wiley John & Sons, Incorporated: New York, 2010. (6) Mhlongo, A. Z.; Alport, M. J. Application of artificial neural network techniques for measuring grain sizes during sugar crystallization. Proc. Annu. Congr. - S. Afr. Sugar Technol. Assoc. 2002, 460−468. (7) Larsen, P. A.; Rawlings, J. B. The potential of current highresolution imaging-based particle size distribution measurements for crystallization monitoring. AlChE J. 2009, 55 (4), 896−905. (8) Larsen, P. A.; Rawlings, J. B. Assessing the reliability of particle number density measurements obtained by image analysis. Part. Part. Syst. Charact. 2009, 25 (5−6), 420−433. (9) Zhang, B.; Abbas, A.; Romagnoli, J. A. Multi-resolution fuzzy clustering approach for image-based particle characterization for particle systems. Chemom. Intellig. Lab. Syst. 2011, 107 (1), 155−164. (10) Meek, G. A. Practical Electron Microscopy for Biologists; Wiley & Sons: London, 1976. (11) Gonzalez, R. C.; Woods, R. E.; Eddins, S. L. Digital Image Processing Using Matlab; Dorling Kindersley: London, 2004. (12) Hagiwara, T.; Wang, H. L.; Suzuki, T.; Takai, R. Fractal analysis of ice crystals in frozen food. J. Agric. Food. Chem. 2002, 50 (11), 3085−3089. (13) Iftekharuddin, K. M.; Zheng, J.; Islam, M. A.; Ogg, R. J. Fractalbased brain tumor detection in multimodal MRI. Appl. Math. Comput. 2009, 207 (1), 23−41. (14) Velazquez-Camilo, O.; Bolanos-Reynoso, E.; Rodriguez, E.; Alvarez-Ramirez, J. Fractal analysis of crystallization slurry images. J. Cryst. Growth 2010, 312 (6), 842−850. (15) Barrera, E.; Gonzalez, F.; Rodriguez, E.; Alvarez-Ramirez, J. Correlation of optical properties with the fractal microstructure of black molybdenum coatings. Appl. Surf. Sci. 2010, 256 (6), 1756− 1763. (16) Kang, T. J.; Kim, S. C.; Sul, I. H.; Youn, J. R.; Chung, K. Fabric surface roughness evaluation using wavelet-fractal method - Part I: Wrinkle, smoothness and seam pucker. Text. Res. J. 2005, 75 (11), 751−760. (17) Kim, S. C.; Kang, T. J. Fabric surface roughness evaluation using wavelet-fractal method - Part II: Fabric pilling evaluation. Text. Res. J. 2005, 75 (11), 761−770. (18) Klinzing, G. E.; Rizk, F.; Marcus, R.; Leung, L. S. Pneumatic Conveying of Solids: A Theoretical and Practical Approach; Springer: New York, 1997. 7018


Image-Based Multiresolution-ANN Approach for Online Particle Size

Recommend Documents