Size measurement and prediction for L-glutamic acid crystal growth during stirred crystallization based on imaging analysis

1.   Introduction

Crystallization is widely utilized as a separation and purification technology in fine chemical fermentation and biochemical engineering industries to obtain high purity solid products. It is also the first step in a series of sub-processes to develop final product ^[1]. In order to guarantee the particle quality, it is necessary to monitor crystal evolution processes for real-time optimization ^[2]. The optimal control of crystal growth during the crystallization stage is critical for the production of crystals. Size feature is very important for crystal production, which should meet specific requirements. Imaging technology facilitating the detection of crystal characteristics is developing rapidly ^[3,4,5]. As an important process analysis technology, image measurement method has been widely applied to the observation of crystal morphology and size ^[2,6]. In order to measure crystal size in the crystallization process, several recent studies were carried out using image analysis methods based on in-situ microscopic imaging systems ^[3]. In-situ imaging systems for detecting crystallization can be divided into the two types: invasive imaging systems ^[7] and non-invasive imaging systems ^[8].

Recently developing on-line image analysis strategies are becoming important to produce crystal products of high quality with desired crystal size distribution (CSD) in an in-reactor environment ^[9,10]. It is well known that the first step of image analysis is to extract the crystal objects, so image segmentation has been one of the most crucial tasks in many researches ^[3,11,12]. Larsen et al. ^[13] studied a robust and efficient segmentation method for measuring CSD information from on-line images of high aspect-ratio crystals. Zhang et al. ^[14] presented a novel image-based approach combining wavelet transform and Fuzzy C-means Clustering (FCM) for particle image segmentation. For the measurement of multi-dimensional sizes, a real-time image analysis method ^[15] was proposed to measure the two-dimension (2D) sizes including length and width of crystals during crystallization. Gao et al. ^[16] utilized a deep-learning detection model to segment LGA crystals and measure the two-dimensional sizes of particles. Ma et al. ^[17] proposed an integrated particle dispersion-imaging system to characterize the growth rate of crystal mean sizes in a cooling crystallization. However, it is worth noting that there are less studies on the prediction for crystal sizes with unchanged operation conditions by imaging system. In addition, most methods did not fully take into account solution turbulence, uneven illumination, and noise in a stirred reactor for crystallization measurement.

In this work, considering the fact that crystals are growing in a stirred reactor under continuous motion, a systematic method based on image analysis is presented for crystal size distribution and size prediction. The method includes image processing (i.e., image restoration and image segmentation), size measurement and size prediction. Figure 1 summarizes all the concrete steps in the developed image analysis model that is applied to crystal images, under complex conditions (e.g., uneven illumination, solution movement etc.). Firstly, for eliminating the influence from solution turbulence, uneven illumination background, and noise, the crystal images are restored from the noised and blurred images by an effective restoration method. Thus, a threshold segmentation is used to extract the particles from the image background. Then the 2D sizes are measured and the corresponding CSD is obtained with probability density estimator. Subsequently, the crystal sizes in the short term are predicted. Experiments on the cooling crystallization process of β-form LGA are conducted to demonstrate the effectiveness of the proposed image analysis method.

The paper is organized as follows. Section 2 shows the experimental set-up. The methodology of image processing for complex problems is introduced in detail in Section 3. In Section 4, size measurement and prediction methods are demonstrated. Section 5 presents the experimental study and the effectiveness of the proposed approach. Finally, conclusions are presented in the last section.

2.   Experimental setup

Experiments were performed with a 1 L glass jacketed reactor, a four-paddle agitator and a Pt 100 temperature probe. For controlling temperature, a Julabo-CF41 thermostatic circulator was employed. The non-invasive imaging system with a high-speed camera was used to record the crystal images during crystallization. The camera was installed outside the glass reactor wall and a LED light source was used to provide illumination. The scheme of the experimental setup is shown in Figure 2.

The β-form LGA crystals were selected for demonstrating 2D sizes (i.e., the length and the width) by the imaging analysis technique. The β-form LGA crystals with a purity of 99% were taken as the solute, and the distilled water was used as the solvent. Seed crystals were produced with β-form LGA crystals by the preparation process including milling, sieving, dissolving, drying etc. The β-form LGA crystals were taken from the solution in a stirred reactor. The agitator stirred at 200 rpm. The solution was heated to 75 ℃ to dissolve all LGA solutes. The solution was cooled down to the temperature 45 ℃ at a speed of 0.5 ℃/min. When the solution reached 45 ℃, seed crystals were added into the reactor. Meanwhile, in-situ crystal images began to be recorded during the crystallization process.

3.   Crystal image processing

The purpose of the crystal image processing method is to obtain the crystals from the recorded image, including two important steps, such as image restoration and segmentation, which are listed in the following subsections.

3.1.   Image restoration

In image acquisition, recorded images are always influenced by Gaussian degradation and noise pollution etc. The images are also blurred for crystal and solution movement in the stirred reactor. It is easy to affect the subsequent particle extraction for crystal images. Therefore, restoration method is considered to solve this problem. The module of plug-and-play alternating direction method of multipliers ^[18] is used for removing image noise and deblurring. Suppose $p(x, y)$ is the observed image, $q(x, y)$ is the unknown clean image, $G(\cdot)$ is a regularization function, $H$ is a blur matrix. The optimization problem is

where $\lambda$ is the regularization parameter.

The inversion step for image deblurring is as

Consequently, the closed form solution is

The blur matrix $H$ is diagonalizable by using the discrete Fourier transform. This model can effectively remove the image fuzziness and noises.

3.2.   Image segmentation

Image segmentation is a key step for crystal size measurement. There were some useful segmentation methods used in crystal image segmentation ^[3,11], but these methods did not focus on the problem of strong non-uniform illumination imaging with low-contrast. In this study, a threshold segmentation algorithm based on morphology is considered due to its speedability.

With regard to the non-uniform illumination of the crystal image $q(x, y)$, the image background is obtained by the opening operation of image morphology ^[19] in this study. Accordingly, the image $T(x, y)$ without uneven background by using the top-hat transformation is computed as

where $b$ is the structure element, $\circ$ represents the opening operation.

The threshold ${t^*}$ is computed from the images by applying Otsu algorithm ^[19]. However, classical Ostu algorithm uses exhaustion method for determining threshold, so the computation amount of Otsu algorithm increases sharply with the increase of image size. Then particle swarm optimization (PSO) algorithm ^[20] is used to search the optimal threshold for improving the efficiency of image segmentation. The fitness function is

where ${\sigma _1}$ is between-class variance, ${\sigma _2}$ is within-class variance in Ostu algorithm.

Thus, the binary image is obtained by the following equation,

where ${\eta _1}$ and ${\eta _2}$ are respectively the background intensity and crystal intensity in the binary image.

After the segmentation, morphological region-filling operation is the next step to fill the holes which are closed due to the reflective and so on. In crystallization processes, due to crystal collision and supersaturating degree in a stirred reactor, there will be prone to breakage, crushing and adhesion for particle movement. Therefore, the tiny pieces and overlapped particles should be excluded from the extracted particles before size measurement. Referring to the particle sieving method ^[15], area threshold and convex parameter are then used to remove meaningless particles.

4.   Crystal size measurement and prediction

4.1.   Size measurement

For estimating the two-dimensional (2D) sizes of β-form LGA crystals, it has been proposed to measure the length and width of the minimum enclosing rectangle of a crystal contour, which denote the length and width of the crystal. The pixel numbers of the length and width of a β-form LGA contour are measured with the best-fit rectangle method ^[21], respectively. Denote by ${L_a}$ the pixel number of measured length and by ${S_a}$ the pixel number of measured width of a crystal contour as defined. The calibrated pixel equivalent ${P_e}$ is computed by the calibration method ^[15]. The real 2D sizes of the crystal, i.e., the real physical length denoted by $L$ and the real physical width denoted by $W$ as shown in Figure 3, are computed by

Based on the image analysis, the crystal sizes can be measured. Because it is impossible to measure the same crystals at different times in a stirring reactor. Therefore, a histogram of CSD is generally plotted from the computed sizes of crystals. However, since the histogram may not be adequate to reflect a continuously differentiable distribution, Gaussian function is used to smooth CSD ^[4,22,23], as to obtain the distribution curve representing the current population size model. For the sizes of length ${l_i}, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = 1, 2, ...$, the probability density estimator of Gaussian function for length variable $l$ is given by

where $l$ conforms to the Gaussian distribution ${N_l}{\rm{(}}{\mu _l}, {\sigma _l}{\rm{)}}$ with mean ${\mu _l}$ and standard deviation ${\sigma _l}$.

For representing the 2D crystal size distribution, the probability density function ${P_{L, W}}$ of 2D Gaussian distribution is considered as

where $(l, {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} w)$ conform to the 2D Gaussian distribution ${N_{l, w}}{\rm{(}}{\mu _l}, {\mu _w}, {\sigma _l}, {\sigma _w}, \rho {\rm{)}}$.

4.2.   Size prediction

According to the current values of crystal sizes, the following size values can be approximately forecasted for predictive control and operation optimization of crystallization processes. Exponential smoothing method ^[24] is a time-series analysis method developed from moving average method. By calculating smoothing values, it can forecast future values with a certain forecasting model. The predicted value is a weighted sum of previous observations with greater weight of new data and smaller weight of old data. In this work, triple exponential smoothing is employed for the short-term trend prediction of crystal size growth. Suppose the smoothed statistic ${S_t}$ is a weighted average of the current observation $\mu$ and the previous smoothed statistic ${S_{t{\rm{ - }}1}}$ for time $t$. For instance, $S_t^{(3)}$ denotes a triple exponential smoothing statistic at time $t$.

These exponential smoothing statistics are given by

where $\alpha$($0 < \alpha < 1$) is the smoothing factor.

The output of the prediction model for the future $M$ values is

The parameters of the above equation are given by

In order to avoid the possible disadvantage of static model parameter, the smoothing factor $\alpha$ is defined dynamically in the study. According to the minimum of the error square sum for historical data, the optimization function is given by

In addition, PSO is used to find the optimal factor $\alpha$ of the previous data which is used for the size prediction.

5.   Experimental results

5.1.   Crystal image processing

The experiment on the crystallization process of β-form LGA was conducted. A typical in-situ image is shown in Figure 4(a). The result of image restoration for deblurring and denoising is shown in Figure 4(b). Figure 5(a) shows the 3D intensity of the image with non-uniform background, which can lead to the difficulty of image segmentation. Figure 5(b) illustrates the the 3D intensity of the image after the removing step of non-uniform background. In Figure 6, the crystal extraction results are obtained by using the image processing method. Figure 6(b) shows the segmentation result for the image of Figure 6(a) with removing non-uniform background. The intensity of crystals was 255 and the intensity of background was 0 in the segmented image. Thus the meaningful crystals were extracted in Figure 6(c) by particle sieving. As is seen from Figure 6(d), the minimum enclosing rectangles of the crystals are labeled and measured by using the best-fit rectangle method. The 2D sizes (i.e., length and width) were directly computed from these rectangles. In addition, for showing the superiority of the proposed method, multi-scale Canny method ^[11], and FCM method ^[14] as the popular crystal segmentation methods are also used for the preprocessed image shown in Figure 4(b). Figure 7 indicates that the proposed method obtains a better result than the other methods for complex conditions.

5.2.   Crystal size result

In the process of crystallization, it is very necessary that real-time size measurement of β-form LGA. For optimization controlling, information of crystal size growth should be provided timely and effectively. In this experiment, the minimum number of particles to be analyzed was 100 for the calculation of CSD. Generally, the growth rate is represented in crystal mean size. Therefore, for showing the measurement accuracy in this study, the comparison of the growth rate in the literatures ^{[17,23,25,26,27]} is given for the length direction in Table 1. These experiments were carried out in different conditions and using different measurement methods. It is seen from Table 1 that these growth rates are in the same order of magnitude, so the size measurement results should be mostly considered as having a good agreement.

In addition, accurate size distribution detection is important for production quality control. For statistical significance of crystal size distribution in the reactor, about 20 images were taken within a minute. The 2D crystal size distribution was established approximating to a 2D Gaussian probability distribution model. Figure 8(a)–(d) show the 2D size distribution evolution for β-form LGA over time for the comprehensive analysis of the 2D size growth. Note that the start time of crystal image capturing was denoted as '0 min'. To show the growth state of the 2D sizes for crystals more clearly, Figure 9 displays the evolution with the growth of the length and width for crystals. From Figure 9, it is indicated that the deviation of Gaussian distribution is growing, describing the size range of crystals become much wider along the crystallization process.

5.3.   Size prediction result

For the sizes of same crystals are so difficult to measure at different times in a stirring reactor, mean size may have statistical significance. In this study, the changes and trends in crystal length and width were investigated for β-form LGA. The measured curves for the mean length ${\mu _l}$ and width ${\mu _w}$ are shown in Figure 10. The mean length and width were fitted using the least squares method. The residual errors of the fitting were about 10%. As shown in Figure 10, the mean length of the crystal grows faster than the mean width. It is seen that the linear regression may be available to show a linear trend of crystal size growth both in length and width. However, it could not be relatively desired for size prediction due to the sharp fluctuation of mean sizes in a short period. Because fine engineering needs accurate real-time operation, more realistic prediction of crystal size is also necessary in a short term with unchanged operation conditions.

It is arduous and complicated to predict the sizes of crystals for a long time in crystallization, due to the changes of supersaturation, hydrodynamics, temperature etc. Therefore, the short-term prediction model for crystal size growth is feasible. For the size short-term prediction of crystals, the prediction results for the mean length and width are shown in Figure 11. In the following two minutes, the six values of crystal sizes were predicted by exponential smoothing model with the set of previous ten values. It is seen from Figure 11 that the prediction points are close to the measured ones which are realistic. The prediction results showed the size short-term prediction model for optimization control under unchanged conditions was effective. According to the trend of size change, the spatial search range of PSO was set as 0.2–0.8. Note that it was two minutes that could be sufficient for the subsequent operation in this experiment.

6.   Conclusions

In this work, a synthetic imaging measurement method, enabling us to estimate evolution of crystal size in a stirred reactor, has been developed for the crystallization case of β-form LGA. The raw crystal images were captured from the non-invasive imaging system. The image processing was used to extract the crystals effectively from the images influenced by uneven illumination and movement of particles and solution. Image segmentation process was speeded up with threshold searching based on PSO. The evolution analysis of 2D sizes (i.e., length and width) was conducted for probability distribution of crystal population. In addition, the short-term mean sizes of crystals were availably predicted with the exponential smoothing method for the subsequent control of crystallization. The experimental results demonstrated that the image analysis method was effective to facilitate successive processing. Future work will investigate a fast image analysis technique to monitor crystal growth rates of individual facets in crystallization processes.

Acknowledgments

This work is supported by the China Postdoctoral Science Foundation (grant number 2020M680979).

Conflict of interest

All authors declare no conflicts of interest in this paper.