Histogram Matching, Hypothesis Testing, and Statistical Control-Chart

Sep 8, 2010 - Hypothesis testing using the Kolmogorov−Smirnof (KS) statistics was also ... This material is available free of charge via the Interne...
0 downloads 0 Views 3MB Size
9932

Ind. Eng. Chem. Res. 2010, 49, 9932–9944

Histogram Matching, Hypothesis Testing, and Statistical Control-Chart-Assisted Nucleation Detection Using Bulk Video Imaging for Optimal Switching between Nucleation and Seed Conditioning Steps Levente L. Simon,*,† Kaoutar Abbou Oucherif,⊥,‡ Zoltan K. Nagy,§ and Konrad Hungerbuhler† Institute of Chemical and Bioengineering, ETH Zurich, W. Pauli Street, 10, Zurich 8093, Switzerland, Department of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47906, Chemical Engineering Department, Loughborough UniVersity, Loughborough, LE11 3TU, United Kingdom

This article investigates the systematic nucleation detection performance of hypothesis testing, histogram matching, and control charts based on bulk video imaging. The control charts are built on the time series of the bin-by-bin distances between the reference and online acquired image histograms. In addition to the binby-bin histogram measures, the cross-bin earth mover’s distance measure has been evaluated. Hypothesis testing using the Kolmogorov-Smirnof (KS) statistics was also investigated as an additional method to detect nucleation. This work is relevant to the crystallization process control with in situ seed generation. Following the nucleation step in unseeded crystallization processes, it is common to keep the temperature constant or to increase it in order to allow the stabilization of the suspension, also called the digestion step. The nucleation detection methods discussed in this work allow the automated and robust switching between the nucleation and digestion steps. It is concluded that histogram matching can be successfully used for systematic nucleation detection, while the KS hypothesis-testing-based detection is strongly dependent on the histogram resolution. Furthermore, it is suggested that histogram matching is performed on the first principal components of color images in order to decrease the autocorrelation of histogram distance time series. 1. Introduction In industrial processes crystallization is a key unit operation, where it is used as a purification and separation step. Within the pharmaceutical industry, there usually exists at least one crystallization step during the drug development stages; thus, emphasizing the importance of understanding and controlling this unit operation. Metastable zone width (MSZW) identification is a prerequisite in the optimum design of crystallization systems where its accurate determination can allow operation within the desired supersaturation range of the material thus ensuring maximal productivity and product quality.1,2 To determine the MSZW several sensors can be considered such as turbidity,3,4 focused beam reflectance measurement,5,6 spectroscopy,5-8 ultrasonic velocity measurements,9 density10 and electrical conductivity monitoring,11 hot stage microscopy,12 quartz-crystal-based monitoring,13,14 and by visual inspection of the crystallizer content.15 MSZW detection by the naked eye still remains one of the most widely used approaches since most of the aforementioned techniques provide information based on local measurements and may be susceptible to errors due to artifacts or nonuniform mixing conditions. Imaging-based nucleation onset and crystallization state monitoring for protein crystallization was evaluated by Kawabata et al.16 The need for robust nucleation onset monitoring and prediction has been recognized by Pollanen et al.8 who proposed multivariate statistical process monitoring charts designed on the basis of IR spectra. * To whom correspondence should be addressed. Tel: +41 44 6334486. Fax: +41 44 6321189. E-mail address: levente.simon@ chem.ethz.ch. † ETH Zurich. ‡ Purdue University. § Loughborough University. ⊥ On research stay at ETH Zurich, Safety and Environmental Technology Group.

Recently bulk video imaging (BVI) was introduced to automate nucleation detection by visual inspection with the aim to minimize the time lag between nucleation of the first crystals and the ability to detect them. The particular feature of the BVI method is that it monitors the formation/dissolution of the solid phase using information from the bulk suspension providing a noninvasive17 or in situ,18 low-cost, automated method for MSZW determination, and making it a well-suited methodology in the pharmaceutical and food industries. Previously, it was shown that BVI can be used for nucleation onset monitoring and metastable zone identification and has exhibited similar performance to the existing focused beam reflectance measurement (FBRM), ultraviolet-visible (UV-vis), infrared (IR) spectroscopy, and calorimetric signals. In the previous works,17,18 one of the suggested image processing methods was to convert the colored images into grayscale and to calculate a mean gray intensity value. The nucleation onset was determined when the signal value rose above a threshold value, which was arbitrarily chosen. The other proposed alternative, which is more complex from a design and an implementation point of view, suggests employing image processing techniques to detect the first crystals as objects.18 Several other works related to image processing with applications in the field of particulate technology focus on image segmentation for particle identification,19-24 while recent advances for the 3D characterization are also reported.25-28 Simon et al.29 proposed a multivariate image analysis, control chart, and acoustic-signal-based method for nucleation detection. Several video data sampling and image processing technologies were evaluated to provide reliable solutions for nucleation detection based on external bulk video imaging. Multivariate image analysis based on statistical concepts was performed on the BVI images and proved to be one of the fastest methods for nucleation onset detection. Image feature descriptors were also evaluated because of their computational simplicity.

10.1021/ie100586p  2010 American Chemical Society Published on Web 09/08/2010

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Moreover acoustic-signal-based monitoring was investigated where gray scale images were converted into a 2 channel stereo sound. The latter method was the least robust at detecting nucleation due to its high degree of subjectivity. Since the BVI method relies on a video signal, the aim of this work is to investigate histogram-matching-based monitoring to provide systematic, sensitive, and robust solutions for nucleation detection. Histogram matching is a widespread technology for image comparison and video data retrieval. It is proposed that whenever feasible, the time series trends generated by the histogram distances are monitored using Shewhart and exponentially weighted moving average (EWMA) statistical control charts to provide systematic nucleation onset detection. To the authors’ knowledge, it is for the first time that the problem of systematic nucleation onset detection is tackled using histogram matching techniques combined with process control charts. Similar to the previous nucleation detection results,29 this work is relevant in the context of unseeded crystallization process control with in situ seed generation. According to this type of operation the seed crystals are formed in situ by generating a high degree of supersaturation in the first stage of the process which causes the appearance of nuclei via homogeneous or heterogeneous primary nucleation. After this first step, it is common to keep the temperature constant in order to allow the stabilization of the suspension30 or to increase it, also called the digestion step.31 In previous works, it was shown that the nucleation onset is a function of the vessel size,32 cooling rate, and stirrer speed; therefore, an automated sensory system is needed in order to detect nucleation and trigger the transition to the next processing step, for example, stopping the cooling and keeping the temperature constant or increase it. Such an automated nucleation detection system is justified for the case of cooling crystallization where the primary nucleation takes place via a heterogeneous mechanism mainly on the subcooled crystallizer walls and the nucleation in the bulk proceeds via contact secondary nucleation33,34 due to the collision of particles with the wall, stirrer, or other particles. For such a process, the nucleation onset is random and has a certain variability, which is also influenced by the presence of foreign solid particles in the liquid phase. To deal with these uncertainties without considering the intervention of operators, in this paper we evaluate several bulk video-imaging-based systematic nucleation detection solutions, which can be integrated in the supersaturation control strategy of batch crystallizers during the nucleation and suspension digestion stage. Another advantage of this automated nucleation detection system is that it is not necessary to perform metastable zone identification experiments in large-scale vessels and the variability of the MSZW is directly handled by the closed-loop switching to the suspension aging or digestion stage. As a result of this nucleation detection method, it is expected that process start-up time is saved, MSZW uncertainties are directly handled via feedback information, and the final product size variability is minimized. Other applications of this systematic nucleation detection method are automated metastable zone or induction time determination using imaging-based laboratory equipment.35 The article is structured as follows: after the description of the motivation and experimental setup, the histogram matching techniques are introduced and discussed. The results and discussions section presents the nucleation detection results, and final conclusions are presented in the last section.

9933

2. Experimental Section The external BVI (eBVI) method is based on recording the crystallization process with a video camera, which captures 25 frames per second. The video camera is placed externally, and the crystallizer vessel is covered with black folia in order to ensure a dark background.17 For large scale applications the camera can be used through an observation window. Anhydrous caffeine from Sigma-Aldrich with a minimum purity of 98% in deionized water was used as the model system. The experimental work which involves endoscopy-based insitu BVI (iBVI) has been carried out using a small scale crystallization calorimeter, and potash alum in water was used as the model substance. The most important feature of the CRCv4 calorimetric vessel is its capability to ensure conditions even in the case of highly exothermic reactions.36-38 According to the power compensation method the liquid temperature is controlled by an electric compensation heater, which ensures good control dynamics. For crystallization purposes which imply decreasing temperature trends, the calorimeter operation is modified so that the advantages offered by the power compensation technique are preserved.18,39 3. Bulk Video-Imaging-Based Nucleation Monitoring After selecting the image sampling rate the first principal component (T1) or the gray image is calculated, which is a twodimensional matrix. The histogram matching is performed by comparing an image showing the clear liquid against new images collected during the experiment. When the calculated distance between the two images generates a time series with low autocorrelation and normal distribution, Shewhart and EWMA process control charts can be employed to perform systematic nucleation detection. 3.1 Video Data Sampling. The video signal sampling may significantly influence the nucleation onset monitoring performance; therefore, this work evaluates the 25 Hz, 1.66 Hz, and 25 Hz sampling strategies with 15 frames overlapped. The 1.66 Hz sampling is used to remove dominating correlation patterns as described in a previous work.29 A RGB color image of M × N × 3 dimensions is defined by the stacked R, G, and B color channels: RGB ) {R, G, B}

(1)

where M is the image width [pixels], N is the height [pixels], R, G, and B are sets of red, green, and blue color components belonging to pixel q. R R ) {xR1 , xR2 , ..., xRq , ..., xMN },

G G ) {xG1 , xG2 , ..., xGq , ..., xMN }, B B B B B ) {x1 , x2 , ..., xq , ..., xMN} (2)

where q is the number of pixels in an image and xRi , xGi , xBi are the intensity values of color pixels in the image. The conversion of RGB images to grayscale (Gr) is performed by the weighted summation of the R, G, and B components: Grn,m ) 0.2989Rn,m + 0.5870Gn,m + 0.1140Bn,m

(3)

The obtained gray scale represents the luminance of the original image and the weights in eq 3 are derived from the NTSC-standard (National Television System Committee) used for analog television signals. The concept of multivariate image analysis and the calculation of T1 images are presented in Appendix A.

9934

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Previously it was observed that the time series trends of the average of gray scale images show oscillatory patterns.17,18 Detailed investigation has shown that these oscillations of the image luminosity are due to the light reflected by the stirrer. Recently, it was shown that working with T1 images decreases the autocorrelation of the time series;29 therefore, the histogram matching is applied on both T1 and gray images. After performing a principal component analysis on a RGB image three principal components are obtained; however, only the first component is kept. The T1 image has the size of the original RGB image (M × N). The R, G, B channels are weighed by the loadings of the first principal component to form the T1 image. The 25 Hz with frame overlapping method is used to concentrate the crystals found in different frames into one frame. The result of stacking or overlapping is an M × N matrix which is obtained by retaining the maximum intensity value along the third dimension and is described as cn,m ) max (T1n,m,f) T1 f)1,...,25

n ) 1, ..., N;

m ) 1, ..., M

(4) c 1c is a two-dimensional matrix. The result of overlapping where T is an image which contains all particles encountered in the 15 frames. The overlapping method is similar to the addition operation often used in image processing. Unlike the addition operation along the pixels where the intensities saturate above 255 and renormalization is needed, when the overlapping method is used the pixel values are always in the 0 to 255 range. This method is also expected to decrease the autocorrelation of the time series data since it performs an operation along the frames. 3.2. Histogram Matching. The nucleation onset detection using histogram matching is justified by the fact that these constitute one of the most popular methods for video data comparison. The information content of colored and gray scale images can be summarized in the form of histograms, which allow for fast image retrieval from image databases and video data. The image histogram is a collection of bins with each bin counting the intensity values of pixels that match the bin value. For image comparison applications, a distance measure is calculated from a reference to a current histogram and the candidate images are ranked according to the minimum distance. The histogram dissimilarity measures are grouped in three classes:40 bin-by-bin, cross-bin, and parameter-based dissimilarity measures. 3.2.1. Nucleation Detection Using Statistical Hypothesis Testing. For nucleation detection purposes the KS statisticsbased hypothesis testing is evaluated,41 and the problem of finding an optimal detection threshold is shifted in the statistical framework by specifying a null-hypothesis and an acceptance confidence interval as follows: H0 ) clear liquid, H1 ) nucleation onset. The two-sample Kolmogorov-Smirnov test compares the cumulative distributions of h1,cum and h2,cum with the null hypothesis that h1,cum and h2,cum are from the same continuous distribution. The test statistic is calculated on cumulative histograms and is given by the KS statistic: dKS )

max |hb1,cum - hb2,cum |

b)0,...,255

(5)

where b is the bin number, h1,cum is the reference, and h2,cum is the cumulative histogram of the observed frame, respectively. The KS test is used to perform statistical hypothesis testing to

Figure 1. Schematic representation of the EMD method for three bin histograms.

check whether the current distribution has been drawn from the reference distribution. The Kolmogorov-Smirnov test is a nonparametric test as it does not require the estimation of distribution parameters such as the mean or standard deviation. 3.2.2. Nucleation Detection Using Bin-by-Bin Histogram Distance Measures. The bin-by-bin dissimilarity measures compare contents of corresponding histogram bins; however, they discard the information about neighboring bins. Bin-bybin measures are sensitive to the histogram position. Some of the most widely used distance measures40 are the KolmogorovSmirnof (KS) distance, histogram intersection, Jeffrey divergence, and χ2 statistics. In addition to hypothesis testing fortythree histogram distance measures as presented by Cha et al.42 are calculated; however, instead of applying hypothesis testing, SPC charts are designed on the obtained time series of the histogram distances. 3.2.3. Nucleation Detection Using Cross-Bin Histogram Distance Measure: The Earth Movers Distance. To overcome the drawbacks presented by the bin-by-bin distances, cross-bin measures have been proposed, which contain additional information about neighboring bins and can be used for partial histogram matching. Some distances belonging to this family are the quadratic-form distance and the earth mover’s distance (EMD) introduced by Rubner et al.40,43,44 The EMD can be applied to histograms or signatures and is the solution of the transportation problem, which is a special case of linear programming (LP). Histograms can be viewed as a special case of signatures in the sense that each bin corresponds to an element in the signature. The histogram values are treated as the weights in a signature, and the bin locations are the cluster positions.40 For a one-dimensional histogram matching problem, the EMD formulation is described below. Given the index set Φ for bins where i is used to denote a bin Φ ) {i : 1 e i e b} and the set of flows is defined as FL ) {fli,j : ( i,j) ∈ Φ} where fli,j denotes a flow from bin i to bin j. P ) {pi : i ∈ Φ} (supply vector) and Q ) {qj : j ∈ Φ} (demand vector) are two histograms to be compared with equal bin numbers. In the context of the EMD the flow is the number of pixels in a certain bin. The EMD objective function is formulated in eq 6 and schematically presented in Figure 1:

( ) ∑ ∑ fl

i,jdi,j

EMD(P, Q, FL) ) min

i

j

∑ ∑ fl

(6)

i,j

i

j

where the ground distance di,j in this work is defined by the L2 norm, and ΣiΣjfli,j is the total flow. The following constraints apply when calculating the EMD: fli,j g 0

∀(i, j) ∈ Φ

(7)

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

9935

Figure 2. Kolmogorov-Smirnof test based probability values for the eBVI and iBVI data using histograms discretized in 64, 128, and 256 bins and 25 Hz sampling; horizontal lines are confidence levels at 1σ (- · - · - · ), 2σ (---), 3σ (s). First visually detected crystal is at 784 s.

∑ fl

i,j

e pi

∀i ∈ Φ

(8)

i,j

e qj

∀j ∈ Φ

(9)

j

∑ fl i

∑ ∑ fl

i,j

i

i

) min(

∑p ∑q) i

i

j

∀(i, j) ∈ Φ

j

(10) The first constraint allows moving “supplies” from P to Q and not vice versa. The next two constraints limit the amount of “supplies” that can be sent and express the constraint that Q cannot receive more “supplies” than the “demand”. The last constraint ensures that the maximum amount of “supplies” possible is moved. The drawback of this method is that for each matching a LP problem must be solved, which may hinder its real-time application in some cases. Note that the number of variables to be optimized (mcnc) is related to the size of the histograms or signatures. The EMD implementation considered in this work uses a histogram size of 50 bins for the eBVI and iBVI images. The solution of the optimization problem is calculated by the large-scale linear programming Matlab solver which implements a linear interior point solution strategy. 3.3. Statistical Process Control Using Univariate Shewhart and EWMA Charts. Statistical process control charts are efficient and established tools for process monitoring, and a detailed discussion with references is found in Simon et al.29 Univariate process control charts have been proposed in the forms of Shewhart, Cumulative Sum (CUSUM), and EWMA charts. These charts differ from each other by applying different weightings to past data samples: for example, the Shewhart chart weights only the last observation, the CUSUM chart assigns equal weights to all observations in a moving window, and the EWMA applies an exponentially decaying weighting.45 The fundamental assumption behind the traditional control charts is that the time series contains independent random variables. Deviation from the assumption of low autocorrelation yields the decrease of the in-control average run length leading to higher false alarm rates. To investigate the autocorrelation

patterns in the frequency domain, frequency plots or periodograms are used and the autocorrelation extent is evaluated by using autocorrelation charts. The second requirement to build a control chart is that the data should be normally distributed around the mean. In this work, normality is assessed by using data histograms and is quantified by the Lillefors test.46 The latter is a two-sided goodness-of-fit test suitable when a fully specified null distribution is unknown and its parameters must be estimated. The default null hypothesis is that the sample comes from a distribution in the normal family, against the alternative that it does not come from a normal distribution. The monitoring robustness of Shewhart and EWMA charts is increased by introducing the following rules: 5 or 10 chart points should be consecutively out of control. All the time series in this work are subject to autocorrelation and normality tests with a 3σ confidence interval. The following Matlab47 toolboxes were used within this project: Optimization Toolbox, Image Processing Toolbox, Statistical Toolbox, and the Time Series Analysis Tool. 4. Results and Discussion For plotting purposes the experimental data are presented about 10-15 min prior to nucleation. To provide an accurate evaluation of the monitoring performance, a visual investigation of the eBVI data was performed and the frame-by-frame analysis has shown that the first crystal appears after 784 s. Note that this crystal was not visible when the video was played at 25 frames/s. For the case of the iBVI data, the nucleation onset analysis based on visual inspection has shown that the first crystal appears after 627 s. In this work the Shewhart control charts have been represented with gray, while the EWMA charts are represented with black. 4.1. Hypothesis-Testing-Based Nucleation Detection. The results of the nucleation onset detection for the eBVI and iBVI data sampled at 25 Hz using the KS statistical hypothesis test are shown in Figure 2. It is observed that the null hypothesis is not rejected for data which monitor the clear liquid, for example, before 800 s. Upon nucleation, the detection performance is a function of the threshold confidence level and the number of

9936

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Table 1. Autocorrelation and Normality Tests Results for the eBVI and iBVI Time Series Obtained from Gray and T1 Images Sampled at 25 Hz time series source

autocorrelation test

normality test

Video Recording Method: eBVI mean of T1 images mean of Gray images

18, 20-23, 31, EMD 31, 36

mean of T1 images mean of Gray images

4, 9, 10, 12, 30, 31, 34, 35, 36, EMD EMD

18, 19, 24, 25, 28, 32-34, 39, 41, 42 9, 35

Video Recording Method: iBVI

bins used to represent the histogram. It is concluded that the higher the histogram resolution is, the faster and larger is the change of the p value. The nucleation detection results are as follows for a 25% confidence level: 957, 957, and 844 s for a 64, 128, and 256 bins discretization, respectively. While this confidence level might seem low for certain applications and for classical hypothesis testing, for the eBVI data it is concluded that low confidence intervals can be used to reject the nullhypothesis since no false alarms are detected. Note that the detection at 844 s is very similar to previous results.29 It is also observed that the larger the histogram resolution is, the lower are the probability values of images that show nucleation, which indicates that the images are considered more different compared to the reference image. It is concluded that for the case when the null hypothesis (the liquid does not contain crystals) has a very high probability value and is valid, the histogram resolution should be the largest possible. In contrast to the observations made for the eBVI data, the iBVI image histograms discretized in 256 bins reject the null hypothesis below a 3σ confidence level. This means that the images collected online are, from a statistical point of view, significantly different compared to the reference image. However, this does not mean that a significant change in the process occurred. This significant difference can also be due to background noise, which in this case is the change of image brightness due to light reflected from the stirrer blades. Note that by decreasing the histogram resolution it is possible to accept the null hypothesis at 2σ and 3σ for a discretization of 64 and 128 bins, respectively. According to this observation, it is concluded that by decreasing the histogram resolution it is possible to increase the histogram similarity probability for images which should be similar, for example, none of them contain crystals, and to perform hypothesis testing around 2σ or 3σ levels. Another observation is that many of the images after nucleation are considered to be more similar to the reference image as they have a high p value; some of the high p values before nucleation could be due to air bubbles in the liquid. The nucleation detection results for the 64 and 128 bin discretization are 754 s (2σ) and 723 s (3σ). 4.2. Histogram Distance and SPCs-Based Nucleation Detection. The second part of the results related to histogram matching techniques refers to the evaluation of 43 distances presented by Cha,42 which were applied to the eBVI and iBVI time series for both the means of the T1 and gray images at different sampling frequencies. Distances 26 and 27 were not included in the evaluation because of their mathematical nature, which may yield complex numbers. The autocorrelation and normality tests were evaluated for all distances with the aim of applying statistical control charts to those that fulfill both tests and without the need of performing any hypothesis testing. The autocorrelation and normality test results are tabulated for those distances which pass the tests. The power density plots, autocorrelation charts, and histogram plots for all time series are found as Supporting Information. Selected control charts are presented for those histogram distances, which pass both

2, 5-7, 9,10, 11, 12 -14, 16, 17, 34, 35, 45 2, 5-7, 10-14, 16-18, 24, 25, 28, 34, 35, 39, 44, 45 Table 2. Nucleation Detection Results Using the Mean of T1 and Gray Images Sampled at 25 Hz. First visually Detected Crystal is at 784 and 627 s for the eBVI and iBVI Data, Respectively histogram distance number Shewhart-5 (s) EWMA-5 (s) EWMA-10 (s) Mean of T1 Images Time Series 18 EMD 9 10 12 35 EMD

eBVI 1098 842 iBVI 643 664 664 643 647

844 842

989 842

638 640 640 639 640

639 640 640 640 642

Mean of Gray Images Time Series EMD

iBVI 647

642

642

tests; an exception being the EMD time series, which did not pass the normality test but for which a control chart was nevertheless designed. 4.2.1. The 25 Hz Sampling. Bin-by-Bin Histogram Matching Results. The results of the autocorrelation and normality tests for T1 and gray images sampled at 25 Hz for the eBVI and iBVI data are shown in Table 1. Numbers in bold show the histogram distances which passed both tests. The histogram distances for which the detection results are presented in figures are given in Appendix B. It is observed that more distances calculated on T1 images pass the autocorrelation test than the distances calculated on gray images, especially for the iBVI case. The nucleation detection results are given in Table 2. It is shown that the nucleation detection performance of distance 18 (Figure 3) is similar to the results obtained in a previous work.29 The time series of distance 18 for the T1 images shows a decreasing trend upon nucleation. To investigate this observation a detailed inspection of the gray and T1 images histograms before and after nucleation is shown in Figure 4. It is observed that upon nucleation the gray image histogram shifts to the right which shows that the overall image brightness increases since the crystals are bright and the background does not change. A different behavior is observed for the T1 images where the overall picture brightness decreases upon nucleation. Also note that the histogram of T1 images is broader compared to that of the gray image. The decreasing trend of distance 18 can be explained by the mathematical nature of the distance, which is the inner product or sum of the dot product of two

Figure 3. SPCs of the mean of T1 eBVI images sampled at 25 Hz for histogram distance 18.

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

9937

Figure 4. Histograms of gray (a) and T1 (b) images before and after nucleation for the eBVI data.

it is not important whether the time series are out of control at the upper or lower control limit.

Figure 5. SPCs of the mean of T1 iBVI images sampled at 25 Hz for histogram distance 9.

histograms. The decreasing distance upon nucleation for the T1 images is due to the shift of the histogram toward lower values which yield darker images. Note that this distance does not contain a difference term and it provides information about the location of the online histogram (the position of the reference histogram is constant). The detection results for the iBVI data are similar for several distances and also for all the control charts with different rules. This is due to the fact that in the iBVI experiments the dynamics of nucleation onset are fast and the image brightness changes quickly (Figure 5). In contrast to the iBVI results, the eBVI images show that initially only few crystals are present in the solution. For the purpose of nucleation detection using SPCs,

Cross-Bin EMD Matching Results. Apart from the case of the eBVI gray images, the EMD distance has also passed the autocorrelation test. The power density plots, autocorrelation charts, and EMD work histograms for the T1 and gray images for the eBVI and iBVI experimental data are presented in Figure 6 and Figure 7. In both cases it is observed that the T1 images, compared to the gray images, yield less autocorrelated time series. However, none of them pass the normality test as shown in the histogram plots. The EMD trend for the eBVI data calculated for T1 images is shown in Figure 8. For an EMD matching task with two histograms of 50 bins each, the optimization problem has 2500 parameters. The average computational time on a 3 GHz Core 2 Duo CPU is 3 and 1.6 s for the T1 and gray eBVI images, respectively. For the iBVI data, the EMD computational time is 2.7 and 2.1 s for the T1 and gray images, respectively. In all cases the tolerance of the objective function was 1 × 10-3. It should be noted that the EMD calculation time on T1 images is larger compared to the EMD on gray images, keeping in mind that T1 histograms are wider compared to gray images. Therefore, it is concluded

Figure 6. Power density plots, autocorrelation charts, and histograms of the EMD of T1 (a, c, e) and gray (b, d, f) eBVI images sampled at 25 Hz EMD.

9938

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Figure 7. Power density plots, autocorrelation charts, and histograms of the EMD of T1 (a, c, e) and gray (b, d, f) iBVI images sampled at 25 Hz EMD.

Figure 8. SPC on the EMD distance of the mean of T1 eBVI images sampled at 25 Hz.

that the wider the histogram is, the more nonzero flows have to be calculated which results in increased computational time. 4.2.2. The 1.66 Hz Sampling. In the previous work29 it was observed that due to light reflections from the stirrer the luminosity of the images changes significantly and a dominating frequency is detected at 1.66 Hz which is believed to be due to the rotation and oscillation of the stirrer shaft. The stirrer speed was set to 100 rpm. To remove this oscillation from the data, a sampling rate of 1.66 Hz was implemented, which means that only every 15th frame was processed. The autocorrelation and normality tests of the T1 and gray time series obtained using the eBVI and iBVI experiments are presented in Table 3. It is observed that all the eBVI and iBVI T1 images time series pass the autocorrelation test. Furthermore, it is observed that the

normality criterion is fulfilled by several histogram distance time series. The power density plots and autocorrelation charts for all histogram distances are found in Figure 9 and Figure 10, where it shown that there are no particular dominating frequencies and no particular patterns in the autocorrelation charts. The normality of the time series can be visually assessed by using histogram plots, as presented in Figure 11, where it is shown that several time series histogram distances fail to exhibit a normal distribution. The fastest detection for the T1 images is shown by histogram 8 and EMD using the EWMA-5 chart, as summarized in Table 4. The SPC charts of distance 8 and 25 built on T1 images of eBVI data are shown in Figure 12. Generally, it is observed that the EWMA-5 chart is the fastest method to detect the nucleation point for all histograms; this is due to the memory effect of the EWMA method. The SPCs built on the basis of gray images show false alarms (FA) and only the Shewhart-5 chart can be used for monitoring. The control chart for the EMD is shown in Figure 13, where it is observed that this method can also be used for nucleation detection and the EWMA-5 chart already signals nucleation at 846 s. The monitoring performance of the distances which pass both tests for the iBVI experimental data are presented in Table 5, while Figure 14 shows the SPC for distance 35 calculated for

Table 3. Autocorrelation and Normality Tests Results for the eBVI and iBVI Time Series Obtained from Gray and T1 Images Sampled at 1.66 Hz time series source

autocorrelation test

normality test

Video Recording Method: eBVI mean of T1 images mean of Gray images

1-45, EMD 9, 12, 18, 30, 31, 34-36, 40

mean of T1 images mean of Gray images

1-45, EMD EMD

8, 19, 20, 23-25, 28, 29, 31-33, 38-43 1-9, 11, 12-14, 16-18, 40, 45, EMD

Video Recording Method: iBVI 1-3, 5-14, 16,17, 19, 24, 28, 29, 32-35, 39-42, 45 1-25, 28, 29, 31-36, 38-45, EMD

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

9939

Figure 9. External bulk video imaging (eBVI) results: power density plots of histogram distances calculated on T1 images sampled at 1.66 Hz; x axis, frequency [Hz], y axis is power of histogram distance; the units are function of the distance measure.

Figure 10. External bulk video imaging (eBVI) results: autocorrelation charts of histogram distances calculated on T1 images. Sampling, 1.66 Hz; x axis, Lags [-]; y axis, autocorrelation coefficients.

the mean of T1 images. Similar to the 25 Hz sampling the nucleation is detected in a small time interval due to the fast change in image brightness. 4.2.3. The 25 Hz Sampling with 15 Frame Overlapping (1.66 Hz). The sampling at 25 Hz with frames overlapping was evaluated with the aim to increase the particle density in the frames. The autocorrelation results presented in Table 6 show that, with the exception of certain distances for the eBVI T1 data, none of the trends passed the autocorrelation test at 3σ. It was observed that the time series contained low-frequency oscillations, which hindered the development of control charts. The normality requirement was fulfilled by some of the distances calculated on gray images. 4.2.4. Final Discussion. It is concluded that for the eBVI data the fastest nucleation detection is performed using a 25 Hz sampling rate. Similar conclusions were drawn in a previous work using the mean of T1 images, which are essentially the

means of the image histogram, as shown in Appendix C. Apart from distance 8 and EMD, the 1.66 Hz sampling shows a detection delay of 106 s. According to the results presented in this work it is found that there is no noticeable trend in the distances or similarity measures that always obey both the normality and autocorrelation tests despite the type of series analyzed or the sampling frequency. One reason for the nucleation detection delay compared to visual observation of the BVI methods based on time series trends of a mean property of an image29 (e.g., mean gray) is due to the averaging operation which was proven to filter small particles.17,18 The method presented in this work uses histograms which are also reduced representations of images and may not be sensitive enough to show the appearance of small crystals which in the initial phase of the process might be only few in an image. Additionally, the signal-to-noise ratio is another important issue which needs consideration. This is the case when

9940

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

Figure 11. External bulk video imaging (eBVI) results: histograms of histogram distances calculated on T1 images sampled at 1.66 Hz. Table 4. Nucleation Detection Results Using the Mean of T1 Images Sampled at 1.66 Hz for the eBVI Data histogram distance number Shewhart-5 (s) EWMA-5 (s) EWMA-10 (s) Mean of T1 Images Time Series 8 19 20 23 24 25 28 29 31 32 33 38 39 40 41 42 43 EMD

1055 1069 1106 1106 1069 1055 1069 1076 1200 < 1069 1069 1076 1106 1078 1069 1069 1200 < 1055

847 954 979 979 954 920 979 954 1046 954 954 955 1022 1002 979 979 1050 848

1000 1002 1052 1052 1049 1007 1049 1049 1112 1049 1056 1052 1052 1007 1049 1054 1056 1002

Mean of Gray Images Time Series 9 12 18 40

996 896 751, FA 895

747, FA 475, FA 450, FA 469, FA

Figure 12. SPCs of the mean of T1 eBVI images sampled at 1.66 Hz for histogram distance 8 (a) and histogram distance 25 (b).

754, FA 482, FA 456, FA 472, FA

the crystals appear as white objects in the images and periodic light reflections from the stirrer blades or from the stirrer blades-crystallizer wall hinder the early nucleation detection. Note that the BVI method aims to gather information from a large volume and the light penetration in the liquid phase should be as large as possible while only allowing reflections from crystals. Another reason for delayed detection using time series (mean gray or histogram distance trends) can be attributed to the control charts which plot the mean of several measurements; therefore, SPCs attenuate the brightness change of frames. Nevertheless, compared to direct nucleation detection methods which identify the crystals as objects in the image18 or in the principal component score space,29 time-series- and SPCs-based methods are much easier to implement and do not require parameter tuning. Thus SPCs are more flexible and do not require

Figure 13. SPC on the EMD distance of the mean of T1 eBVI images sampled at 1.66 Hz.

algorithm retuning in case the substance to be crystallized, the vessel size, and illumination conditions change. Other histogram matching methods could also be evaluated, for example, parameter-based dissimilarity measures, which compute a set of parameters specific to the histogram, such as first three moments, peaks, means, and standard deviations. An alternative method for nucleation detection based on digital images is the employment of artificial intelligence techniques for pattern matching.48 Further extensions of the nucleation

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010 Table 5. Nucleation Detection Results Using the Mean of T1 and Gray Images Sampled at 1.66 Hz for the iBVI Data histogram distance number Shewhart-5 (s) EWMA-5 (s) EWMA-10 (s) Mean of T1 Images Time Series 1 2 3 5 6 7 8 9 10 11 12 13 14 16 17 19 24 28 29 32 33 34 35 39 40 41 42 45 EMD

686 649 686 650 649 648 650 648 664 649 649 649 653 649 652 655 660 660 689 655 655 646 647 651 668 655 655 651 653

676 649 673 648 645 645 645 649 644 649 645 645 645 645 647 649 655 656 677 649 648 644 644 646 661 649 649 646 648

682 679 654 651 653 655 650 653 656 652 652 651 652 652 653 655 661 662 682 655 654 650 650 651 667 655 655 651 654

detection methods discussed in this paper could find applicability in the context of direct nucleation control49 (DNC) or adaptive supersaturation control50 strategies. 5. Conclusions The work presented in this article investigates the systematic nucleation detection performance of hypothesis testing, histogram matching, and control charts based on bulk video imaging (BVI). The control charts are built on the time series of the bin-by-bin distances between the reference and online acquired image histograms. In addition to the bin-by-bin histogram measures, the cross-bin earth mover’s distance measure has been

Figure 14. SPC on the distance 35 of the mean of T1 iBVI images sampled at 1.66 Hz.

9941

evaluated. A further method to detect nucleation is to perform hypothesis testing using the Kolmogorov-Smirnof statistics. Histogram matching techniques can be used for nucleation detection when applied on the principal component of a RGB image. The best detection performance was achieved at the maximum sampling rate which is 25 Hz. It is suggested that histogram matching is performed on the first principal components (T1) of color images in order to decrease the autocorrelation of histogram distances time series. It should be noted that computing the EMD on a principal component image is computationally demanding making its application impractical at high data sampling rates and the detection capabilities are similar to certain bin-by-bin distance measures. Furthermore, it was found that the KS hypothesis-testing-based detection is strongly dependent on the histogram resolution. Finally, it is concluded that using the EWMA chart with the rule of five consecutive chart points out of control provides the fastest control chart setting. Supporting Information Available: Autocorrelation charts, power density plots, and histograms of the distance time series. This material is available free of charge via the Internet at http:// pubs.acs.org. Appendix A Low-cost video cameras use a single light sensitive sensor and a color filter array (CFA) with each pixel unit recording the intensity value of each color channel typically red, green, and blue. The most widespread CFA is the Bayer filter.51 In this scheme the green filters are in a quincunx, interlaced grid with the red and blue completing the rest of the empty cells. Note that the CFA captures only one-third of the needed color intensities and it provides a 2D color matrix. Furthermore, the green color filters make up half of the cells while the rest is shared by the blue and red filters. This is because the luminance response curve of the eye peaks at around the frequency of green light. As a next step, this 2D structure is converted into the well-known 3D or RGB color format by an interpolation method. This process is called debayering or demosaicking. Some of the interpolation methods are nearest neighborhood, bilinear interpolation, smooth hue transition, adaptive color plane interpolation, and gradient-based interpolation.52 One source of correlation might be due to the fact that only every third pixel value is a true measurement and the rest of the pixel values are obtained by interpolation. Furthermore, image compression algorithms may introduce additional correlation among the color channels; therefore, multivariate image analysis has been proposed to compress and extract uncorrelated information using multiway principal component analysis.53-55 MPCA decomposes the unfolded multivariate image array U into principal components made of score matrices ta, loading vectors pa, and into a residual matrix E: A

Table 6. Autocorrelation and Normality Tests Results for the eBVI and iBVI Time Series Obtained from Gray and T1 Images Sampled at 25 Hz and Overlapping of 15 Frames (1.66 Hz) time series source

autocorrelation test

normality test

Video Recording Method: eBVI 19, 20, 23-25, 28-33, 36-45 mean of Gray images 1-14, 16-18, 37, 40, 43, 45, EMD mean of T1 images

Video Recording Method: iBVI mean of T1 images 34, 35 mean of Gray images 4, 20, 44

U)

∑t p

T a a

+E

(A.1)

a)1

where A is the number of principal components. The multiway PCA implementation proposed for the nucleation onset monitoring follows the methodology described in previous works,54,56 and the loading vectors are computed using the efficient kernel algorithm.54 The kernel matrix (UTU) is first constructed and a singular value decomposition (SVD) is performed on the matrix to extract the loading vectors pa, and the corresponding score vectors ta are calculated according to

9942

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

ta ) Upa

(A.2)

The calculation results have been validated against the results provided by the MACCMIA Matlab tool57 developed at McMaster University, Canada. Several MIA-based applications in the context of process engineering have been reported: online monitoring of snack food processes,58,59 flames monitoring,60,61 flotation froth monitoring,62 anode cover alumina content63 and mineral ore composition estimation,64 counterfeit drug detection,65 emulsion monitoring,66,67 quality control of polymer blends,68 semiconductor characteristics monitoring,69 table coating uniformity,70 and recently crystal24 and fiber size distribution determination.71

and the sum of all pixel values in an image is 255

∑N b b

(C.3)

b)0

By using the notation introduced above, the mean gray intensity is calculated as follows: 255

Nb

∑ ∑p

b n

b)0 n)1 255 b

(C.4)

∑N b)0

Appendix B Distance measures for which SPCs are shown in the manuscript:

while the histogram mean is defined as follows: 255

∑N b

Kulczynski (8):

b

d

∑ |h

b 1

b)0 255

- hb2 |

b)1 d

dkul )

∑N

(B.1)

∑ min(h , h ) b 1

(C.5) b

b)0

b 2

Using eqs C.2 and C.3 it is shown that eqs C.4 and C.5 are equivalent.

b)1

Canberra (9): |hb1 - hb2 |

d

dCan )



Nomenclature

(B.2)

hb1 + hb2

b)1

Inner Product (18): d

∑h h

dIP )

b b 1 2

(B.3)

b)1

Bhattacharyya (25): d

dBhat ) -ln

∑ √h h

b b 1 2

(B.4)

b)1

Clarck (35): dClk )

∑( d

|hb1 - hb2 |

b)1

hb1 + hb2

)

2

(B.5)

where d is the number of bins in histograms. In this work d ) 256.

A ) number of principal components b ) intensity value or histogram bin number d ) histogram distance E ) residuals matrix fl ) flow FL ) set of flows Gr ) gray image h1 ) reference histogram h2 ) observed histogram M ) image width [pixels] N ) image height [pixels] or number of pixles in a histogram bin p ) loading vector or histogram bin value P ) histogram t ) score matrix T1 ) first principal component ps ) position in a signature q ) histogram bin value Q ) histogram x ) pixel intensity value along a color channel U ) multivariate image

Appendix C

Greek Symbols

The sum of pixel values of a gray image can be calculated by the addition of the number of pixels for all intensity values b:

λ ) EWMA memory factor [unitless] Φ ) Index set χ2 ) chi-square distance

p01 + ... + p0n + ... + pN0 b + ... + pb1 + ... + pbn + ... +

Indices

255

pNb b + ... + pb1 + ... + pbn + ... + pN255 b )

Nb

∑ ∑p

b n

(C.1)

b)0 n)1

b

N pnb where Nb is the number of pixels of intensity b, the Σn)1 term is the sum of pixel values corresponding to intensity value b. The multiple summation of the same value can also be expressed in the form of a multiplication:

N

b

∑p

b n

n)1

) Nbb

(C.2)

a ) principal component f ) image i ) pixel along image height j ) pixel along image width q ) pixels in an image b ) bin number AbbreViations eBVI ) external bulk video imaging EMD ) earth movers distance EWMA ) exponentially weighted moving average

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010 FBRM ) focused beam reflectance measurement iBVI ) in-situ bulk video imaging KS ) Kolmogorov-Smirnov statistic MSZW ) metastable zone width RGB ) color image defined by the red, green, and blue color channels SPC ) statistical process control UV-vis ) ultraviolet-visible SVD ) singular value decomposition

Literature Cited (1) Braatz, R. D. Advanced control of crystallization processes. Annu. ReV. Control 2002, 26, 87. (2) Nagy, Z. K.; Chew, J. W.; Fujiwara, M.; Braatz, R. D. Comparative performance of concentration and temperature controlled batch crystallizations. J Process Control 2008, 18 (3-4), 399. (3) Parsons, A. R.; Black, S. N.; Colling, R. Automated measurement of metastable zones for pharmaceutical compounds. Chem. Eng. Res. Des. 2003, 81 (A6), 700. (4) De Anda, J. C.; Wang, X. Z.; Lai, X.; Roberts, K. J.; Jennings, K. H.; Wilkinson, M. J.; Watson, D.; Roberts, D. Real-time product morphology monitoring in crystallization using imaging technique. AIChE J. 2005, 51 (5), 1406. (5) O’Grady, D.; Barrett, M.; Casey, E.; Glennon, B. The effect of mixing on the metastable zone width and nucleation kinetics in the antisolvent crystallization of benzoic acid. Chem. Eng. Res. Des. 2007, 85 (A7), 945. (6) Fujiwara, M.; Chow, P. S.; Ma, D. L.; Braatz, R. D. Paracetamol crystallization using laser backscattering and ATR-FTIR spectroscopy: Metastability, agglomeration, and control. Cryst. Growth Des. 2002, 2 (5), 363. (7) Anderson, J. E.; Moore, S.; Tarczynski, F.; Walker, D. Determination of the onset of crystallization of N1-2-(thiazolyl)sulfanilamide (sulfathiazole) by UV-vis and calorimetry using an automated reaction platform; subsequent characterization of polymorphic forms using dispersive Raman spectroscopy. Spectrochim. Acta, Part A 2001, 57 (9), 1793. (8) Pollanen, K.; Hakkinen, A.; Reinikainen, S. P.; Rantanen, J.; Minkkinen, P. Dynamic PCA-based MSPC charts for nucleation prediction in batch cooling crystallization processes. Chemom. Intell. Lab. Syst. 2006, 84 (1-2), 126. (9) Gurbuz, H.; Ozdemir, B. Experimental determination of the metastable zone width of borax decahydrate by ultrasonic velocity measurement. J. Cryst. Growth 2003, 252 (1-3), 343. (10) Marciniak, B. Density and ultrasonic velocity of undersaturated and supersaturated solutions of fluoranthene in trichloroethylene, and study of their metastable zone width. J. Cryst. Growth 2002, 236, 347. (11) Lyczko, N.; Espitalier, F.; Louisnard, O.; Schwartzentruber, J. Effect of ultrasound on the induction time and the metastable zone widths of potassium sulphate. Chem. Eng. J. 2002, 86 (3), 233. (12) Kumar, F. J.; Moorthy, S. G.; Jayaraman, D.; Subramanian, C. Estimation of metastable zone width, interfacial energy and growth rates of KTiOPO4 crystallizing from K6P4O13 flux by hot stage microscopy. J. Cryst. Growth 1996, 160 (1-2), 129. (13) Loffelmann, M.; Mersmann, A. How to measure supersaturation. Chem. Eng. Sci. 2002, 57 (20), 4301. (14) Joung, O. J.; Kim, Y. H.; Fukui, K. Determination of metastable zone width in cooling crystallization with a quartz crystal sensor. Sens. Actuators, B 2005, 105 (2), 464. (15) Sohnel, O.; Mullin, J. W. The role of time in metastable zone width determinations. Chem. Eng. Res. Des. 1988, 66 (6), 537. (16) Kawabata, K.; Takahashi, M.; Saitoh, K.; Asama, H.; Mishima, T.; Sugahara, M.; Miyano, M. Evaluation of crystalline objects in crystallizing protein droplets based on line-segment information in greyscale images. Acta Crystallogr., Sect. D 2006, 62, 239. (17) Simon, L. L.; Nagy, Z. K.; Hungerbuhler, K. Comparison of external bulk video imaging with focused beam reflectance and ultraviolet-visible spectroscopy for crystallization nucleation detection and metastable zone identification. Chem. Eng. Sci. 2009, 64, 3344. (18) Simon, L. L.; Nagy, Z. K.; Hungerbuhler, K. Endoscopy-based in situ bulk video imaging of batch crystallization processes. Org. Process Res. DeV. 2009, 13 (6), 1254. (19) Calderon De Anda, J.; Wang, X. Z.; Roberts, K. J. Multiscale segmentation image analysis for the in-process monitoring of particle shape with batch crystallisers. Chem. Eng. Sci. 2005, 60 (4), 1053.

9943

(20) Eggers, J.; Kempkes, M.; Mazzotti, M. Measurement of size and shape distributions of particles through image analysis. Chem. Eng. Sci. 2008, 63 (22), 5513. (21) Larsen, P. A.; Rawlings, J. B. The potential of current highresolution imaging-based particle size distribution measurements for crystallization monitoring. AIChE J. 2009, 55 (4), 896. (22) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. An algorithm for analyzing noisy, in situ images of high-aspect-ratio crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61 (16), 5236. (23) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. Model-based object recognition to measure crystal size and shape distributions from in situ video images. Chem. Eng. Sci. 2007, 62 (5), 1430. (24) Sarkar, D.; Doan, X.-T.; Ying, Z.; Srinivasan, R. In situ particle size estimation for crystallization processes by multivariate image analysis. Chem. Eng. Sci. 2009, 64 (1), 9. (25) Darakis, E.; Khanam, T.; Rajendran, A.; Kariwala, V.; Naughton, T. J.; Asundi, A. K. Microparticle characterization using digital holography. Chem. Eng. Sci. 2010, 65 (2), 1037. (26) Khanam, T.; Darakis, E.; Rajendran, A.; Kariwala, V.; Asundi, A. K.; Naughton, T. J. On-line digital holographic measurement of size and shape of microparticles for crystallization processes. Proc. SPIE 2008, 7155, art. no. 71551K. (27) Kempkes, M.; Darakis, E.; Khanam, T.; Rajendran, A.; Kariwala, V.; Mazzotti, M.; Naughton, T. J.; Asundi, A. K. Three dimensional digital holographic profiling of microfibers. Opt. Express 2009, 17 (4), 2938. (28) Kempkes, M.; Vetter, T.; Mazzotti, M. Measurement of 3D particle size distributions by stereoscopic imaging. Chem. Eng. Sci. 2010, 65 (4), 1362. (29) Simon, L. L.; Abbou Oucherif, K.; Nagy, Z. K.; Hungerbuhler, K. Bulk video imaging based multivariate image analysis, process control chart, and acoustic signal assisted nucleation detection. Chem. Eng. Sci. 2010, 65 (17), 4983. (30) Kim, H.; Bang, Y.; Lee, K. S.; Yang, D. R. Use of calorimetry model and batch control technique for scale-up of unseeded batch cooling crystallization of poly(hydroxybenzophenone). Ind. Eng. Chem. Res. 2009, 48 (14), 6776. (31) Harner, R. S.; Ressler, R. J.; Briggs, R. L.; Hitt, J. E.; Larsen, P. A.; Frank, T. C. Use of a fiber-optic turbidity probe to monitor and control commercial-scale unseeded batch crystallizations. Org. Process Res. DeV. 2009, 13 (1), 114. (32) Yi, Y. J.; Myerson, A. S. Laboratory scale batch crystallization and the role of vessel size. Chem. Eng. Res. Des. 2006, 84 (8), 721. (33) Bauer, L. G.; Rousseau, R. W.; McCabe, W. L. Influence of crystal size on rate of contact nucleation in stirred-tank crystallizers. AIChE J. 1974, 20 (4), 653. (34) Garside, J.; Davey, R. J. Secondary contact nucleationskinetics, growth and scale-up. Chem. Eng. Commun. 1980, 4 (4-5), 393. (35) Pino-Garcia, O.; Rasmuson, A. C. Primary nucleation of vanillin explored by a novel multicell device. Ind. Eng. Chem. Res. 2003, 42 (20), 4899. (36) Visentin, F.; Gianoli, S. I.; Zogg, A.; Kut, O. M.; Hungerbuhler, K. Pressure-resistant small-scale reaction calorimeter that combines the principles of power compensation and heat balance (CRC.v4). Org. Process Res. DeV. 2004, 8 (5), 725. (37) Billeter, J.; Neuhold, Y. M.; Simon, L.; Puxty, G.; Hungerbuhler, K. Uncertainties and error propagation in kinetic hard-modelling of spectroscopic data. Chemom. Intell. Lab. Syst. 2008, 93 (2), 120. (38) Richner, G.; Neuhold, Y. M.; Papadokonstantakis, S.; Hungerbuhler, K. Temperature oscillation calorimetry for the determination of the heat capacity in a small-scale reactor. Chem. Eng. Sci. 2008, 63 (14), 3755. (39) Richner, G.; Neuhold, Y.-M.; Hungerbuhler, K. Nonisothermal calorimetry for fast thermokinetic reaction analysis: solvent-free esterification of n-butanol by acetic anhydride. Org. Process Res. DeV. 2010, 14 (3), 524. (40) Rubner, Y.; Tomasi, C.; Guibas, L. J. The earth mover’s distance as a metric for image retrieval. Int. J. Comput. Vision. 2000, 40 (2), 99. (41) Massey, F. J. The Kolmogorov-Smirnov test for goodness of fit. J. Am. Soc. Stat. Assoc. 1951, 46 (253), 68. (42) Cha, S.-H. Comprehensive survey on distance/similarity measures between probability density functions. Int. J. Math. Models Meth. Appl. Sci. 2007, 1 (4), 300. (43) Rubner, Y.; Puzicha, J.; Tomasi, C.; Buhmann, J. M. Empirical evaluation of dissimilarity measures for color and texture. Comput. Vision Image Understand. 2001, 84 (1), 25. (44) Rubner, Y.; Tomasi, C.; Guibas, L. J. A metric for distributions with applications to image databases. 6th Int. Conf. Comput. Vision 1998, 59. (45) Hunter, J. S. The exponentially weighted moving average. J. Qual. Tech. 1986, 18 (4), 203.

9944

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

(46) Lilliefors, H. On Kolmogorov-Smirnov test for normality with mean and variance unknown. J. Am. Soc. Stat. Assoc. 1967, 62 (318), 399. (47) Matlab, version 7.6.0.324(R2008a); The Mathworks, Inc: Natick, MA, 2008, www.mathworks.com. (48) Simon, L. L.; Hungerbuhler, K. Industrial batch dryer data mining using intelligent pattern classifiers: Neural network, neuro-fuzzy and Takagi-Sugeno fuzzy models. Chem. Eng. J. 2010, 157 (2-3), 568. (49) Abu Bakar, M. R.; Nagy, Z. K.; Saleemi, A. N.; Rielly, C. D. The impact of direct nucleation control on crystal size distribution in pharmaceutical crystallization processes. Cryst. Growth Des. 2009, 9 (3), 1378. (50) Woo, X. Y.; Nagy, Z. K.; Tan, R. B. H.; Braatz, R. D. Adaptive concentration control of cooling and antisolvent crystallization with laser backscattering measurement. Cryst. Growth Des. 2009, 9 (1), 182. (51) Bayer, B. E. Color imaging array. U.S. Patent 3,971,065, 1976. (52) Ramanath, R.; Snyder, W. E.; Bilbro, G. L.; Sander, W. A. Demosaicking methods for Bayer color arrays. J. Electron. Imaging 2002, 11 (3), 306. (53) Esbensen, K.; Geladi, P. Strategy of multivariate image-analysis (MIA). Chemom. Intell. Lab. Syst. 1989, 7 (1-2), 67. (54) Geladi, P.; Isaksson, H.; Lindqvist, L.; Wold, S.; Esbensen, K. Principal component analysis of multivariate images. Chemom. Intell. Lab. Syst. 1989, 5 (3), 209. (55) Geladi, P.; Grahn, H., MultiVariate Image Analysis; John Wiley & Sons: Chichester, U.K., 1996. (56) Bharati, M. H.; MacGregor, J. F. Multivariate image analysis for real-time process monitoring and control. Ind. Eng. Chem. Res. 1998, 37 (12), 4715. (57) Dunn, K. MACCMIA1.81; McMaster University: Hamilton, Ontario, Canada, 2007. (58) Yu, H. L.; MacGregor, J. F.; Haarsma, G.; Bourg, W. Digital imaging for online monitoring and control of industrial snack food processes. Ind. Eng. Chem. Res. 2003, 42 (13), 3036. (59) Yu, H. L.; MacGregor, J. F. Multivariate image analysis and regression for prediction of coating content and distribution in the production of snack foods. Chemom. Intell. Lab. Syst. 2003, 67 (2), 125. (60) Yu, H. L.; MacGregor, J. F. Monitoring flames in an industrial boiler using multivariate image analysis. AIChE J. 2004, 50 (7), 1474. (61) Szatvanyi, G.; Duchesne, C.; Bartolacci, G. Multivariate image analysis of flames for product quality and combustion control in rotary kilns. Ind. Eng. Chem. Res. 2006, 45 (13), 4706.

(62) Liu, J. J.; MacGregor, J. F.; Duchesne, C.; Bartolacci, G. Flotation froth monitoring using multiresolutional multivariate image analysis. Miner. Eng. 2005, 18 (1), 65. (63) Tessier, J.; Duchesne, C.; Gauthier, C.; Dufour, G. Estimation of alumina content of anode cover materials using multivariate image analysis techniques. Chem. Eng. Sci. 2008, 63 (5), 1370. (64) Tessier, J.; Duchesne, C.; Bartolacci, G. A machine vision approach to on-line estimation of run-of-mine ore composition on conveyor belts. Miner. Eng. 2007, 20 (12), 1129. (65) Rodionova, O. Y.; Houmoller, L. P.; Pomerantsev, A. L.; Geladi, P.; Burger, J.; Dorofeyev, V. L.; Arzamastsev, A. P. NIR spectrometry for counterfeit drug detectionsA feasibility study. Anal. Chim. Acta 2005, 549 (1-2), 151. (66) de Juan, A.; Maeder, M.; Hancewicz, T.; Tauler, R. Use of local rank-based spatial information for resolution of spectroscopic images. J. Chemom. 2008, 22 (5-6), 291. (67) de Juan, A.; Maeder, M.; Hancewicz, T.; Tauler, R. Local rank analysis for exploratory spectroscopic image analysis. Fixed size image window-evolving factor analysis. Chemom. Intell. Lab. Syst. 2005, 77 (12), 64. (68) Gosselin, R.; Rodrigue, D.; Gonzalez-Nunez, R.; Duchesne, C. Potential of hyperspectral imaging for quality control of polymer blend films. Ind. Eng. Chem. Res. 2009, 48 (6), 3033. (69) Facco, P.; Bezzo, F.; Barolo, M.; Mukherjee, R.; Romagnoli, J., A. Monitoring roughness and edge shape on semiconductors through multiresolution and multivariate image analysis. AIChE J. 2009, 55 (5), 1147. (70) Garcı´a-Mun˜oz, S.; Gierer, D. S. Coating uniformity assessment for colored immediate release tablets using multivariate image analysis. Int. J. Pharm. 2010, 395 (1-2), 104. (71) Tomba, E.; Facco, P.; Roso, M.; Modesti, M.; Bezzo, F.; Barolo, M. Artificial vision system for the automatic measurement of interfiber pore characteristics and fiber diameter distribution in nanofiber assemblies. Ind. Eng. Chem. Res. 2010, 49 (6), 2957–2968.

ReceiVed for reView March 11, 2010 ReVised manuscript receiVed July 21, 2010 Accepted August 9, 2010 IE100586P