Image-Based in Situ Identification of Face Specific Crystal Growth

Article pubs.acs.org/crystal

Image-Based in Situ Identification of Face Specific Crystal Growth Rates from Crystal Populations Christian Borchert,† Erik Temmel,† Holger Eisenschmidt,§ Heike Lorenz,† Andreas Seidel-Morgenstern,†,‡ and Kai Sundmacher†,§,* †

Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstraße 1, D-39106 Magdeburg, Germany Otto von Guericke University, Chemical Process Engineering, Universitätsplatz 2, D-39106 Magdeburg, Germany § Otto von Guericke University, Process Systems Engineering, Universitätsplatz 2, D-39106 Magdeburg, Germany ‡

ABSTRACT: Online microscopy has received much attention in the field of crystal shape control over recent years, since commonly used measurement techniques cannot provide enough information for this purpose. In this work, we present an estimation scheme that serves to reconstruct the 3D crystal shape from the measured 2D crystal projection. The boundary curves of the crystal projections are parametrized by Fourier descriptors, which are subsequently compared with a precomputed database. The procedure is evaluated in comparison with various effects that might impair the estimation. A good agreement between the true and estimated values is found in all cases. The presented methods are applied to batch cooling crystallizations of potassium dihydrogen phosphate (KDP) dissolved in water. As a result, the face specific growth rates are determined as a function of supersaturation. To validate the performance of this scheme, the calculated face specific growth rates are used in a model-based prediction of the supersaturation profiles, which agrees well with the experimental data.

1. INTRODUCTION Crystal shape is a product property that is desirable to be controlled with regard to various process steps and also with regard to the physico-chemical properties of the crystalline material. For instance, filterability, dissolution behavior, and conveying performance of crystallized solids are key factors for the design of downstream processes and are decisively influenced by the crystal shape. The evolution of crystal shapes is a process that is kinetically controlled. 1 Since the underlying growth kinetics and mechanisms are sensitive to the process conditions, it is desirable for control purposes to observe a process with regard to the crystal shape evolution rather than to rely on model predictions only. On the other hand, observations can serve as the basis for model calibration and identification. A method to directly measure the 3D geometry of a crystal is given by tomography. Here a wave is generated that penetrates the crystal and is partially absorbed. From every direction a different image is obtained, and the entire set of collected images is used to reconstruct the 3D crystal shape, which is represented by a 3D array of voxels.2−4 A numerical value (grayscale level) is assigned to every voxel to distinguish between the solid and the surrounding phase as well as to resolve the inner structure of the solid material. An example of the application of such a technique can be found in Merrifield et al.5 Here the authors combined highresolution sychrotron X-ray microtomography and X-ray diffraction to characterize powders in a size range of 20−50 μm to characterize the crystal shape and the molecular orientation of the crystal faces simultaneously. © 2014 American Chemical Society

Another method for the direct measurement of the 3D crystal shape is given by optical sectioning. A crystal is mechanically dissected into thin slices and fixed in a glassy solid. The reconstruction of the 3D crystal shape can then simply been done by stacking the 2D slices to a 3D object. This technique offers as well the opportunity to analyze the inner structure of the crystal and is frequently applied in volcanology and mineralogy.6−8 A nondestructive approach to sectioning is the use of a confocal microscope. The slicing is achieved by altering the focal plane of the microscope. The collected images can afterward be processed in order to reconstruct the 3D crystal shape. Successful applications of this technique in the field of crystallization can, for instance, be found in Singh et al.,9 Castro et al.,10 Conchello and Lichtman,11 and Wilson.12 Both approaches, tomography and sectioning, require a thorough sample preparation, which is timeconsuming and may alter the probe, as well as costly devices. Despite the direct measurement of the 3D crystal shape, several approaches have been presented in the literature that are concerned with the reconstruction of the crystal shape from microscopic images. Due to an easier sample collection, especially these approaches have received much attention over the last years,13 and in principle, three different sampling technologies can be distinguished: (a) Of f line Imaging of Dry Samples. A crystal sample is collected from the process, filtered, dried, and dispersed on a Received: July 19, 2013 Revised: December 13, 2013 Published: January 13, 2014 952

dx.doi.org/10.1021/cg401098x | Cryst. Growth Des. 2014, 14, 952−971

Crystal Growth & Design

Article

obtain 2D shape distributions of the population. They discuss particularly the transient behavior of crystal habits in an industrial crystallizer. Patchigolla and Wilkinson26 analyzed images of α and β L-glutamic acid and monosodium glutamate monohydrate to measure size and shape distributions. They compare the measured size distribution with results obtained from manual microscopic image analysis, ultrasonic attenuation spectroscopy, and laser diffraction spectroscopy and find all measurements in reasonable agreement. Alander et al.27,28 investigate with the help of offline microscopic images the agglomeration behavior of paracetamol (acetaminophen) in different solvents. Not only scalar shape measurements of the identified particles but also the intensity variation of the grayscale within the particle region are taken. Sucrose crystallization from water has been investigated using a custom flow-cell equipped with a microscope camera with regard to crystal morphology by Ferreira et al.18 The recorded data is used for a mass-based growth rate determination at different impurity concentrations (calcium chloride and D-fructose). Similarly to the aforementioned study of Alander et al.28 also the agglomeration behavior in chemically differently environments is studied. A couple of papers with qualitative image analysis studies have been published by the groups of Wang and Roberts: In De Anda et al.,29 the segmentation of images obtained from an in situ imaging probe and subsequent identification of particle regions is discussed. De Anda et al.30,31 used a pragmatic approach to distinguish α and β polymorphs of L-glutamic acid on the basis of Fourier descriptors together with a neural network. In Wang et al.,32 the needle shaped habit of β-L-glutamic acid crystals is approximated by equivalent cuboids varying in length and width. Neglecting the influence of crystal orientation, this shape model was used to identify the crystal growth laws for these individual directions, using the recorded video frames33 or external samples of the seed, and the final crystal size distribution,34 respectively. 1.2. Estimation of the State from 2D Projection. When the 3D geometrical state of the crystal has to be reconstructed from a single 2D image, only model based methods can be used in contrast to nonparametric shape quantifications of true 3D sensors. The model is needed because features of the projection can then be used as an indicator for the 3D shape, that is, the space of possible objects throwing the projection is confined to those that are as well obtainable by the shape model. Ma et al.35 determined the face specific growth rates of potash alum crystals in a hot-stage reactor. Assuming that all crystals were photographed from a known and constant perspective, he was able to derive characteristic lengths from a detailed shape model. By tracking of the evolution of these characteristic lengths over time, face specific growth rates for this model substance could be obtained. Shape recognition in grayscale in situ images of a crystal suspension has, for instance, been discussed by Larsen et al.36 and Larsen.37 The method presented therein mainly relies on the detection of linear features in the grayscale image and is designed for well faceted high aspect ratio crystals. A further development is capable of reconstructing wire frame models of the crystals with which crystal morphologies can be estimated.22 Mazzotti and coworkers developed a flow-through microscope in which the passing particles can be seen from two different directions.17 Based on former work,16 which estimated the height and width of needle-like crystals from a photograph from a single direction, which was also applied to ascorbic acid crystallization, they estimate the height, width, and depth distribution, that is, a

microscope slide. This sample can then be analyzed using a standard laboratory microscope. Here careful sample preparation is required, in order not to alter the sampled crystals.14 Due to this time-consuming sample preparation, this technique is not suited for process control, since it can be applied neither in real time nor in situ or in-line. (b) Inline Imaging of the Suspension. The suspension is pumped through a cuvette, which is installed on a microscope stage via an external sampling loop. With this technique, the crystallization process can in principle be monitored in real time, since no additional time for sample preparation is required. If the cuvette is carefully designed, the observed crystals will align in the focal plane of the microscope, which typically results in better image quality compared with in situ probes. The major drawback of these techniques is the external sampling loop, which is sensitive to clogging or encrustation if not designed and tempered properly. Several devices of this kind have been presented in the literature ranging from custommade systems as presented by Patience and Rawlings,15 Eggers et al.,16 Kempkes et al.,17 and Ferreira et al.18 to commercially available systems as presented by Borchert and Sundmacher.19 (c) In Situ Imaging of the Suspension. A video camera is installed directly in the crystallizer. Therefore, the problems related to the external sampling loop of the previous method are avoided. Compared with inline imaging systems, however, the quality of the acquired images is lower in the sense that only few of the photographed crystal particles will be located within the focal plane of the camera. If the suspension density is high, different overlapping particles can hardly be distinguished, and unlike the case for in-line imaging, it is not possible to dilute the photographed suspension. Encrustation can pose a serious problem for the in situ probe, see, for example, Li et al.,20 which can be avoided by installing the camera outside the crystallizer at an imaging window, as, for example, demonstrated by Li et al.21 or Larsen et al.22 The further processing of images can be distinguished in two directions. The first one deduces shape and size measurements from the 2D projection and uses this information for further applications, for instance, for process monitoring, model calibration, or quality control. This utilization of imaging techniques is a rather simple but powerful means to gain better insight into a crystallization process. The second direction that can be identified aims at reconstructing the true 3D shape from the 2D projection. 1.1. Simple 2D Observation. Imaging or better to say visual inspection of particles is the oldest and most trusted analysis practice in mechanical process engineering in general and specifically in crystallization. The automated observation of the evolving crystal morphology has been of interest to researchers for years over a widespread class of substances. Hence, the following survey cannot be comprehensive but claims to give a representative selection. Size distribution analysis is of course the most prevalent task in particle image processing. For example, nucleation detection based on histogram matching and hypothesis testing,23 as well as principal component based multivariate image analysis,24 was done by Simon and co-workers using bulk video imaging. Monnier et al.25 use automatic image analysis in order to determine the size distribution of samples that are taken online for adipic-acid crystallization. The observed size distribution along with a calorimetric supersaturation measurement has been used to identify a full population balance model of the crystallization process. Puel et al.14 measure hydroquinone crystals dispersed on a microscope slide in length and width to 953



Article

Burge.40 They are useful in quantifying the shape of a projection in a pragmatic and simple way. Particularly, shape distributions can be produced in a straightforward way, and the structure evolution of particles can be observed in such a manner, as shown, for example, by Borchert and Sundmacher.19 Multiple scalar shape descriptors can be assembled to a descriptor vector, which is a point in the descriptor space involving, for example, also a quantity for size. Though scalar shape descriptors can be used for classification, they are a rather coarse measure compared with the rich information that is actually contained in the contour. For the identification of the 3D state, a descriptor carrying more detailed shape information is requisite. Clearly, this information must be extracted from the boundary curve whose detailed deduction from the 3D object is discussed next. In sections 2.2−2.4, the signature and the abstraction of the same, making up an adequate information-rich descriptor vector, is introduced. 2.1. 2D Boundary Curve. The shape of the crystal is given by a point set defined by the geometrical state vector h that contains the perpendicular distances of the crystallographic faces to the crystal center, see Figure 1. The orientation of the

3D property distribution, for glutamic acid and ascorbic acid crystals using a 4D axis length distribution.38 Schorsch et al.39 extended the existing approaches by classifying the observed crystals as spherical, cuboid, or needle-like objects. Using generic properties of these shape classes, they were able to assign the corresponding characteristic lengths of these classes to each observed crystal. By this approach, Mazzotti and coworkers were able to reduce the necessary computation time significantly and to avoid the assumption of a uniformly distributed orientation distribution of the crystals. The aim of this work is to develop a scheme with which the 3D shape can be estimated from 2D photographs of the crystals. Based on a detailed crystal shape model, this scheme allows for the simultaneous estimation of the crystal shape together with the orientation from which the crystal was photographed. The methods are shown to be robust to a finite CCD-chip resolution, imperfect image segmentation, and crystal asymmetry. The proposed scheme can be applied to images from any video microscope. Here a flow-through microscope has been chosen because it yields high quality images with regard to particle sharpness.19 The only information that is extracted from the photo to estimate the 3D shape is the boundary curve (contour) of the 2D particle projection. If the 3D crystal is not attached to others, that is, it is a single, convex crystal, the projection from any orientation is a convex region as well. If the Miller indices of the faces are known, the set of projection shapes from different directions of a crystal with fixed geometrical state can be precomputed. This knowledge can in turn be used conversely to relate the 2D projection to the 3D object. Of course, this inversion is not always unique, that is, the estimation is subject to possibly large, but quantifiable, errors. Being able to observe the shape evolution allows then for the extraction of growth kinetics. This article is organized as follows. In section 2, the specification of the 2D projection in terms of parameters, socalled shape descriptors, is presented. From the shape descriptors of the 2D projection, the 3D shape shall be obtained, which is subject of section 3, which includes the especially important analysis on the error of this estimation scheme. In section 4, the proposed estimation scheme is applied to images recorded in experiments. This information can be used to identify simple models as demonstrated in section 5. Finally, section 6 concludes the work.

Figure 1. Quantification of the crystal shape using the distances h of the crystallographic faces.

faces is characterized by their unit normals whose rows make up the matrix N. The space (coordinate r) that is occupied by the crystal is thus S = {r: Nr ≤ h}

(1)

Alternatively, this set can be expressed as the convex domain bound by the crystal vertices, the so-called V-representation, see Ziegler:41

2. SHAPE DESCRIPTORS Crystals observed with microscopes are viewed from one perspective, and essentially only a projection of the particle can be seen. The projection is characterized by its boundary curve or contour. This curve is determined by the shape of the crystal and its orientation in space. The same crystal shape can produce clearly different contours depending on the orientation from which the projection is taken. However, different crystal shapes projected onto the 2D plane from the same direction usually produce distinct projections. Hence, the projection can be used as an indicator for the shape of the original 3D crystal shape. But the inversion from the projection to the 3D shape is usually not unique. For a convex crystal, the boundary curve is determined by the convex hull of the projection of the crystal vertices, which are connected by straight lines. That is, the crystal projection exhibits features reflecting the characteristics of the original 3D state. Essentially, scalar (solidity, aspect ratio, convexity, circularity, eccentricity) and vector-valued shape descriptors are distinguished. Definitions of scalar descriptors are, for instance, given by Gonzalez and Woods,2 Alander et al.,28 and Burger and

nS

S = {r: r =

nS

∑ αj vj, ∑ αj = 1} j=1

j=1

(2)

where the crystal vertices, vj ∈ 3, can be calculated from the plane orientations, N, and the geometrical state vector, h, see refs 1, 41, and 42. The number of vertices is denoted by nS and can be a function of the geometrical state h, but this is not further discussed here, see ref 42 for details. The projection of S onto the photoplane is given by nP

P = {x : x =

nP

∑ αj pj ,

∑ αj = 1}

j=1

j=1

(3)

where x = (x, y)T ∈ 2 is the coordinate vector on the projection plane and pj ∈ 2 are projected vertices of the crystal: pj = Pvj (4) 954



Article

Figure 2. Relationship between 3D crystal shape and Fourier descriptors. (left) Projection of a potassium dihydrogen phosphate (KDP) crystal onto a projection plane. The orientation of the crystal in terms of the Euler angles is ψ = [1.12,2,1], the geometrical state is hT = [0.177, 0.464]. The crystal is projected on an array simulating a CCD chip whose pixel size is 5 × 10−6 in both directions. The overall number of pixels of the boundary is 1485, which is a value comparable to well developed crystals as recorded in experiments. (right) The signature function ‘centroid to boundary versus angle’ of the projection evaluated for a sequence of 128 angles. Major features like corners can be clearly identified (top). The absolute values of the Fourier coefficients of the signature yield a description independent of rotations on the projection plane (bottom).

see Figure 2 (left). The projection matrix is composed of a true projection matrix and a product of rigid body rotation matrices, which specify the relative orientation between the crystal and plane in terms of the Euler angles ψ = (ψ,ϑ,ϕ):43,44

P = PPPψ

A=

(5a)

(5b)

xC =

⎡1 0 0 ⎤ ⎢ ⎥ ⎢ 0 cos ϑ −sin ϑ ⎥ ⎢⎣ 0 sin ϑ cos ϑ ⎥⎦

(5c)

(7)

j=1

1 6A

nP

∑ (pj + pj + 1)det[pj , pj + 1],

pn

j=1

r = |x − x C |

P+1

= p1 (8)

(9a)

x ⎧ +arctan ⎪ ⎪ θ=⎨ ⎪−arctan x ⎪ ⎩

B = {x : x = (1 − λ)pj + λ pj + 1, 0 ≤ λ ≤ 1, j = 1 ... nP}, p

nP

∑ det[pj , pj + 1]

With eq 6, the projection boundary B is deduced from the state vector h and the orientation. The next section introduces a more handy representation of this information. 2.2. Signature. The boundary surrounding the 2D projection is a 1D curve and can thus be represented by a 1D functional that can be defined in different ways.2 The simplest of these so-called signature functions is obtained by measuring the distance between the centroid of the object and the boundary curve as a function of the angle, see Figure 2. This signature function is referred to as centroid-to-boundary distance.45 It is transformed to polar coordinates whose origin is the centroid of the projection:43

see also Figure 2. In eq 3, not all nS vertices of the crystal are used but only the nP ones that form the convex hull of the projection, that is, those that lie on the boundary of the projection. Let the projected vertices, as Figure 2 depicts, be ordered in a counterclockwise direction, that is, the boundary of the projection is composed of linear sections between two successive projected vertices. With this, the boundary can be expressed as

pn + 1 = p1

1 2

and the centroid, xC:37,41

⎡1 0 0⎤ PP = ⎢ ⎣ 0 1 0 ⎥⎦ ⎡ cos ψ −sin ψ 0 ⎤ ⎢ ⎥ Pψ = ⎢ sin ψ cos ψ 0 ⎥ × ⎢⎣ 0 0 1 ⎥⎦ ⎡ cos ϕ −sin ϕ 0 ⎤ ⎢ ⎥ × ⎢ sin ϕ cos ϕ 0 ⎥ ⎢⎣ ⎥ 0 0 1⎦

Properties of the projection that can be computed from the projected vertices efficiently and shall be used below are the area, A,40

(6) 955

− xC for (y − yC ) ≥ 0 r − xC for (y − yC ) < 0 r

(9b)



Article

2.4. Abstraction of the Signature: The Fourier Transform. The centroid-to-boundary distance is in principle a good representation form of the projection shape. It characterizes the boundary in any detail, and in an identification scheme for the 3D state, it would be possible to take all these details into account. However, the boundaries in real crystal images are subject to noisy variations. In general, the crystal shape is not ideally formed but featured with irregular bumps, for instance, caused through smaller crystals attached to the surface or growth irregularities. Furthermore, the real projection boundary is not recorded on an infinitely well resolving plane but on a finitely resolving CCD-chip. That is, recording the boundary of a real crystal with a digital microscope allows only for a limited representation accuracy of subtle boundary features. For this reason, it is desirable to transfer the centroid-to-boundary function to a representation in which slight variations, compared with the major shape features, are captured with descriptors that assume values of a magnitude well below those descriptors that carry the principal shape data. Since the centroid-to-boundary function is a 1D signature of the 2D projection contour, the well developed 1D signal processing tools can be used in order to analyze and classify the projection. The Fourier transform is the most widely used concept in signal processing and has in fact been used for shape analysis some 50 years ago for the first time, see Zahn and Roskies46 and references therein. From the 1970s onward Fourier descriptors have been employed for pattern recognition. Of the numerous applications, the identification of handwritten characters has been the most vividly studied problem; see, for example, the papers by Zahn and Roskies,46 Granlund,47 and Persoon and Fu.48 A major difference, whether advantageous or disadvantageous, is that Fourier descriptors do not resolve the contour locally but all major and minor details are distributed globally onto the Fourier coefficients. Thus, a small number of Fourier coefficients suffices to reconstruct the gross shape.2 The Fourier transform of the sampled signature is given by43

For a convex polygon defined by the edges p1, ..., pnP, the function r(θ) is composed of linear connecting lines between the edges: r (θ ) = −

det[pj , pj + 1] (p Tj + p Tj + 1)q

for

0≤−

p Tj q (p Tj + p Tj + 1)q

≤1 (10a)

⎛−sin θ ⎞ q=⎜ ⎟, ⎝ cos θ ⎠

0 ≤ θ ≤ 2π

(10b)

The functional r(θ) is invariant under translation of the object on the projection plane but clearly depends on the orientation. This can be remedied by a transformation displacing r(θ) such that it starts at a characteristic feature, for instance, at the maximum. Such an operation corresponds to a shapepreserving rotation on the photoplane. A typical example of the centroid-to-boundary signature is shown in Figure 2. The signature of the projection on the left is depicted in the upper right corner. It exhibits two maxima corresponding to the projected apexes of the pyramids and two minima reflecting the drawn-out prismatic part of the crystal. Alternatives to the centroid-to-boundary approach use the tangent angle versus arc length function, directly a complex function in the x−y-plane, or the distance between the centroid and special feature points on the boundary, for instance, points of high curvature, which are associated with corners as can be seen in Figure 2.2,45 However, these variants require additional processing steps (e.g., an algorithm and decision criteria to detect high curvature points), which reduces the robustness compared with the centroid-to-boundary distance signature as used in the following. 2.3. Sampled Signature. So far the signature of the projection contour has been discussed in continuous terms, that is, the form given in eq 10 can be evaluated at an arbitrary accuracy. In practice, the boundary can be recorded only at a finite resolution on a CCD-chip. That is, the boundary is in practice given as a set of mean pixel coordinates, bj = (xj,yj)T ∈ 2 , which make up the discretized boundary, B = [b1, ..., bb]

N−1

Rk+1 =

(11)

n=0

k = 0, ..., N − 1 (14)

see also Figure 3. The sampling of the signature obtained from the discrete digital image is obtained from measuring the distance from the centroid to the boundary in the respective direction

where i denotes the imaginary unit. The Fourier coefficients Rj are in general complex-valued and specify amplitude and phase of periodic components in the original sequence r, see eqs 12 and 13. That is, a shift of the signature function, corresponding to a rotation on the projection plane, is reflected in the phase of the Fourier transformed signal. Therefore, a rotation-invariant

r(B) = [dist(x C, B , θ1), ..., dist(x C, B , θN )]T ∈ N (12)

where dist (·) is the Eucledian distance between xC and the intersection between (i) the ray with angle θ starting at xC and (ii) the straight line connecting the two boundary points bj and bj+1 that are closest to the ray, see Figure 3. N is the number of sampling points of the signature, which was in this work set to N = 128. The sampled signature function not taken from a discretized image but directly from the analytical form of the boundary, eq 9, evaluated at predefined angles θ1, ..., θn is simply obtained from solving eq 10 at the respective angles: r(B) = [r(θ1), ..., r(θN )]T ∈ N

⎛ 2πikn ⎞ ⎟, N ⎠

∑ r(θn+ 1) exp⎜⎝−

(13)

Of course, this sampling of the signature function is not used for the analysis of images but shall be utilized to build up a database with which images are analyzed.

Figure 3. Evaluation of the signature functional at defined angles for pixelated images. 956



Article

Figure 4. Computational flowcharts of the database generation (left) and the shape estimation routine (right).

Figure 5. “Real world” effects considered in the assessment of the estimation scheme: pixelated image (a), blurred crystal boundary (b), and asymmetric crystals (c).

and the orientation ψ, that is, the descriptor vector can be seen as a function of these inputs:

shape description is obtained by considering only the magnitude of the Fourier coefficients: |R k | =

Re(R k)2 + Im(R k)2 ,

k = 1, ..., N − 1.

d = fdescr(h, ψ )

(15)

In order to achieve scale invariance, the Fourier coefficients are resized by division through the largest one. The so obtained values make up the shape descriptor vector d = (d1 , ..., dN )T =

1 (R1 , ..., RN )T , s

s=

(17)

3. SHAPE ESTIMATION SCHEME Clearly, a loss of information is caused by the projection, and in general there is not a unique mapping from the descriptor vector d to the state vector h. That is, the same set of descriptors can be produced from different states and orientations so that the inverse operation to eq 17,

1 max(R j) j

(16)

(h, ψ ) = fdescr −1 (d)

which shall be used in the subsequent section for the identification of the 3D shape. As described in sections 2.1−2.4, the descriptor is obtained from the 3D state quantified for a particular crystal system by the geometrical state vector h

(18)

must not necessarily be unique. But if the geometrical state is given, it is an easy task to project the vertices on a plane and determine the descriptor vector d. In Borchert and 957



Article

Figure 6. Estimator performance tested on a KDP population being concentrated in three different areas of the state space. Each “island” consists of 200 individual crystals that have been photographed in silico from a random perspective. The images are taken with infinite resolution, that is, the contour is given in an analytical form, eq 10, and subsequently sampled according to eq 13. Left column, Given shape distribution with sketches of representative shapes for each region (top). Estimated distribution on the basis of synthetic images from which Fourier descriptors are extracted and used to match it with the closest entry in a precomputed database (bottom). Right column, Distributions of the relative error of the estimated values for h1 and h2 separately for all three regions in which the crystals are concentrated, that is, sample 1 (top), sample 2 (middle), and sample 3 (bottom). The size of the used database was nlut = 45000 (nlut,h = 300 and nlut,ψ = 150).

Sundmacher,49 we have described a shape estimation scheme that uses a precomputed database compiled from synthetic crystal images. Given a geometrical state h and orientation ψ relative to the photoplane, the descriptor vector d in terms of the Fourier transform of the projection boundary is determined as explained in section 2. This is repeated by randomly varying the crystal orientation and state, see Figure 4, left. All resulting descriptor vectors and the respective geometrical state and orientation are stored in a lookup table (or database). The so acquired lookup table is used to compare descriptor vectors of crystal projections that have been recorded from crystal images. The descriptor vector in the database that is closest to the measured one is taken as the hit, and the assigned geometrical state vector is taken as an estimate for the 3D state of the crystal that can be seen on the photo, see Figure 4, right. Since all entries of the database were generated by choosing random geometrical states and orientation, all

entries in this database will be unique. Therefore, we can ensure that the comparison of a descriptor vector from a recorded crystal projection to the descriptor database will lead to unique estimates for the geometrical state and the crystal orientation. The accuracy of this estimation scheme depends on the complexity of the crystal, the quality of the images and the degree of consistency with which the real crystals are reflected in the polyhedral crystal model. The following section assesses the quality of the estimates taking different effects into account that deteriorate the image quality. The following case study is based on synthetic images; that is, it serves to compare the measured with the actual state to test the estimation achievable in practice. A population consisting of three subpopulations of potassium dihydrogen phosphate (KDP) crystals, which are all normally 958



Article

Figure 7. Estimator performance tested on a KDP population being concentrated in three different areas of the state space. Each “island” consists of 200 individual crystals, which have been photographed in silico from a random perspective. The images of the asymmetric crystals are taken with finite resolution and noise added to the projection boundary and are subsequently sampled according to eq 13. Left column, Given shape distribution with sketches of representative shapes for each region (top). Estimated distribution on the basis of synthetic images from which Fourier descriptors are extracted and used to match it with the closest entry in a precomputed database (bottom). Right column, Distributions of the relative error of the estimated values for h1 and h2 separately for all three regions in which the crystals are concentrated, that is, sample 1 (top), sample 2 (middle), and sample 3 (bottom). The size of the used database was nlut = 45000 (each crystal resulting from the 300 state vectors has been imaged from 150 different perspectives).

distributed with means and variances and vanishing covariances, is examined.

The assessment of the estimation scheme is carried out by adding various “real world” effects (Figure 5) that make the estimation more difficult due to the resulting lower data quality. After examining the results for the highly idealized case of infinite chip resolution and perfectly imaged symmetric crystals, section 3.1, these effects are discussed in section 3.2. The effect of the database size on the estimation performance is analyzed in section 3.3. 3.1. Projection on a Plane with Infinite Resolution. Figure 6 (right) shows the error distribution for each of the three populations for both quantities h1 and h2. It can be seen that the majority of the estimates come with an error well below 5%. The standard deviation of the error distribution ranges from 5.12 × 10−2 to 1.5 × 10−1 as specified in the titles of the error distributions shown in Figure 6. Note also that the standard deviation of the error distribution is relatively high because erratic outliers for which the estimate is very poor exist in all populations. This is because it is not always possible to

μ1T = (2, 6) × 10−4 σ1,11 = 2 × 10−5 σ1,22 = 2 × 10−5 μ2T = (4, 3) × 10−4 σ2,11 = 5 × 10−6 σ2,22 = 1.25 × 10−5 μ3T = (2, 3) × 10−4 σ3,11 = 1 × 10−5 σ3,22 = 1 × 10−5 (19)

From each subpopulation, 200 samples are generated, see Figure 6 (top left), and randomly oriented, and subsequently the descriptor vector d = fdescr(h, ψ) is determined. The above given procedure is used to reestimate the geometrical state. The error of the estimate of the jth particle with respect to the quantity hk is measured by ejk ,est =

hk ,est − hk ,true hk ,true

(20) 959



Article

Figure 8. Mean error between the given and estimated crystal shape distribution as a function of database size.

3.2. Asymmetric Crystals Projected on a Plane with Finite Resolution and Blurred Boundaries. Clearly in practical imaging, the resolution of the crystal boundary is finite, and furthermore, it may be deteriorated, that is, blurred. In order to quantify these effects, at first, the images are taken on a plane with infinite resolution, which is then converted to an image with finite resolution (discretization) of the projection plane. The projection of the crystal can in practice become diffuse as depicted in Figure 5b due to two effects: (i) small crystals adhering to the “actual” large crystal or (ii) crystals are imaged out of focus. The crystal model that shall be identified is in a particular point highly idealized: In stirred systems, the crystals rarely grow perfectly symmetric but the distances of individual faces of one form usually fluctuate around a mean value, see Figure 5c. The joint impact of blurred projection boundaries and asymmetric crystals is shown in Figure 7. The error distributions are broader compared with the previously presented ideal case. It is this magnitude of error that we shall in the worst case be concerned with in experiments. As will be shown later, the blurring of the crystal boundary as used in these studies (Figure 5b) is far worse than that for the processed experimental images. The standard deviation of the error can reach values as high as 20%; however, the mean values of the populations (quantified by the almost zero median) are still surprisingly well reproducible. 3.3. Database Size. Of course, the estimation errors given in sections 3.1 and 3.2 are dependent not only on the quality of the observed crystals but on the size of the applied database as well. This size is defined by the number of simulated crystals, nlut,h, the number of orientations, nlut,ψ, from which these crystals were photographed, and the number of sampling points

reconstruct the 3D shape uniquely. For instance, if the three representative shapes, shown in Figure 6 (top left), are viewed from top, that is, only pyramidal faces can be seen, the state of these pyramidal faces cannot be determined. Such outliers can thus only be eliminated by analyzing the inner structure of the particle region or increasing the number of perspectives from which the crystal is imaged, see the works of Larsen37 and Kempkes.50 Both approaches are not further pursued in this work. The median as well as the mean of the error distribution is for all error distributions below 1 × 10−2, that is, below 1%; the mean value of a population, wherever located in the state space, can be measured with a high accuracy. Hence, the error distribution indicates that the estimation scheme is not working equally well in all areas of the state space. Particularly population 3 can be identified with a better quality than populations 1 and 2. The accuracy of the estimates of h2 in populations 1 and 2 is rather poor compared with estimates of h1. In population 1, this is due to the lower prominence of the h2 faces on the crystal surface, see examples of these shapes in Figure 6 (top left). Therefore, features of h2 faces are less expressed in the projections and thus the distinction between subtle differences in h2 requires highly resolved lookup tables to achieve better results. However, a standard deviation of the error, denoted by σek, of less than 10% shall be considered sufficient for our purposes especially with regard to the well measurable mean values of the population. This result has been achieved with a lookup table that was made up of nlut,h = 300 different state vectors whose resulting crystal shapes are each photographed from nlut,ψ = 150 different random directions. Next the case of pixelated, blurred images of asymmetric crystals is presented. 960



Article

apparatus for the image acquisition is presented in section 4.1. Since the so obtained images cannot be fed directly to the shape estimation scheme, several image processing steps have to be carried out, which are shortly sketched in section 4.2. After this, two case studies on evolving KDP crystal populations are presented in section 4.3. 4.1. Image Acquisition. Of the numerous methods and devices with which crystals could be imaged, the flow through microscope QICPIC (Sympatec GmbH, (Clausthal, Germany)) has been selected. This instrument records up to 25 grayscale images per second of the passing suspension in a flow cell of 2 mm width, see Figure 9. The field of view is 5000 μm by 5000 μm (also other modes are available but have not been used in our experiments) and is projected onto a CCD chip with a resolution of 1024 pixels by 1024 pixels. The objects are illuminated with laser light, which is transmitted through the suspension in the flow cell. Due to the uniformity of the distance to the camera, the particles are relatively well focused, see Figure 11 and also the discussion in the Introduction. The suspension that flows through the microscope’s cuvette is continuously withdrawn from a vessel in which the suspension is well stirred in order to minimize biasing of the sampling with respect to crystal size, see Figure 10. The peristaltic pump is operated at 200−250 mL/min. In principle, it can happen that the crystals break due to mechanical stressing in the pump. However, since the crystal system under consideration, KDP, has a relatively high Vickers hardness of 150 in contrast to comparable model substances such as potassium alum (56) or sucrose (64),51 breakage could not be observed even for prolonged continuous pumping. 4.2. Image Analysis. For the quantitative analysis of the images that are acquired with the flow through microscope, an image processing routine has been assembled that accounts for

Figure 9. Working principle of the flow-through microscope.

of the signature N. To assess the quality of the database, 1000 randomly sampled crystals from each subpopulation given by eq 19 were projected on artificially pixelated images with a resolution of 5 μm/pixel. The mean relative error of the estimated crystal sizes was used as a quality measure and computed for databases of various sizes. The results are depicted in Figure 8. The mean estimation error is about 10% with small databases, and stagnates at approximately 4% with larger databases. Especially the dependency of the estimation error on the number of sampling points N is worth noting. At low values of N, not enough details of the boundary curve are recognized, and no good estimation results can be reached, even for large databases. This effect is decreasing, as more sampling points are used to characterize the boundary curve, but at higher values of N (i.e., N = 128 and N = 256), no substantial effect on the estimation error is observable. In the following sections, the size of the database was set to nlut,h = 300, nlut,ψ = 150, and N = 128 to ensure a reasonably good estimation quality at moderate computational costs.

4. APPLICATION TO REAL IMAGES The shape estimation method that has been developed and examined in the previous section shall now be applied to real images rather than testing it only on synthetic ones. The

Figure 10. Sketch of the experimental setup. From the batch crystallizer, the suspension is continuously pumped through the cuvette in which the crystals are imaged. The suspension vessel can be tempered, which allows cooling crystallization experiments to be performed. 961



Article

Table 1. Parameters of Experiments 1 and 2 quantity

unit

symbol

mass of solvent

kg

mH2O

mass of dissolved KDP mass of seeds temperature gradient QICPIC frame rate pH value source of KDP, purity seed preparation

kg kg K/h Hz

mKDP,diss mKDP,seed ΔT f

value experiment 1 2

value experiment 2 2

0.682 0.682 5 × 10−4 5 × 10−4 −20 −5 10 10 3.99 3.99 Merck, CAS-No. 7778-77-0, ≥99% sieve fraction 150−200 μm taken from the as-delivered crystal population

performed. To extract quantitative information, and especially the contours of single crystals for the shape identification, from the images, the following main steps are taken out and discussed in subsequent sections: (1) image enhancement, (2) thresholding and region filling, and (3) particle identification and measurement. Details on the processing routine can be found in ref 19. The algorithm has been implemented in MATLAB 2010B, which was equipped with the IMAGE PROCESSING TOOLBOX V6.4. 4.3. Two Case Studies on KDP Crystal Populations. Up to now, the theoretical background of the image processing and shape identification has been discussed, which is in this section applied to real experimental data. A seeded batch cooling crystallization has been taken out for the system potassium dihydrogen phosphate (KDP) (crystalline phase) and water (solvent). In Figure 10, the experimental setup is sketched. The jacketed crystallization vessel has a volume of about 3 L, and a thermostat allows the application of temperature profiles to the vessel content. The crystal suspension is continuously pumped

Figure 11. Example of a processed image of experiment 1 at t ≈ 1500 s.

the special requirements due to the relatively dark background and partly transparent crystals, see Figure 11. The advantage over the software that comes with the device is full access to all image abstraction parameters necessary to identify particle regions and full access to all single-particle measures and their user-defined manipulation. Also the analysis of multivariate distribution data with respect to number or mass density can be

Figure 12. State of the continuous phase for experiment 1. Due to the relatively high cooling rate (left), the supersaturation (right) becomes large.

Figure 13. Evaluation of the amount of nuclei produced over the course of the experiment. (left) particles that are smaller than the average sphere equivalent diameter of the single crystals minus 200 μm are classified as nuclei. (right) Nuclei mass fraction. 962



Article

Figure 14. Total number of particles and number of particles classified as single crystals (left) and the number fraction of single crystals (right), measured within a sampling interval of 25 s for experiment 1.

Figure 15. Evolution of the shape distribution of experiment 1 using a kernel density estimator. 963



Article

An example image is shown in Figure 11. Colored lines are traced around identified particles. Red stands for particles that have been classified as aggregates (solidity < 0.95), and green refers to areas that have been identified as single crystals (solidity ≥ 0.95). Additionally, a list of numbers beneath or within each particle is added, which are in the depicted image from top to bottom: circle equivalent diameter, major and minor axes lengths, solidity, and eccentricity. All these descriptors can be calculated using the MATLAB IMAGE PROCESSING TOOLBOX’s command regionprops. The numbers drawn onto the processed image as well as the particle boundaries serve to check for the consistent operation of the algorithm. Though for storing the acquired information it is not necessary to endorse the data onto the image that has been processed, we have continuously filed every tenth so prepared image for visual inspection. The experiments we have undertaken have been inspired by the ones reported by Yang et al.52 They have shown that the shape of the KDP crystals grown from water−KDP solution strongly varies with the temperature profile and seed load, and thus supersaturation, that is applied. We have performed two experiments to verify their findings that at high supersaturations KDP assumes a more compact shape, see section 4.3.1 and that for lower supersaturations the form becomes more elongated, see section 4.3.2. Before the experiments are further discussed in a more technical detail, the unsurprising major results, confirming the results of Yang et al.,52 can be seen in the evolution of the mean states depicted in Figures 16 and 22. The harsher temperature profile of the first experiment, see Figure 12, yields higher supersaturation values, which produces relatively compact crystals (Figures 11 and 16). The relatively gentle temperature decrease of the second experiment, see Figure 16, leads to a lower supersaturation level producing more elongated crystals than in the first experiment (Figure 22). 4.3.1. KDP Growth at High Supersaturation (Experiment 1). Basic technical data of experiment 1 can be found in Table 1. The temperature profile was set to linear cooling of −20 K/h, which was not fully achieved as can be seen in Figure 12 (left). The seed crystals were added when the solution was slightly undercooled at t ≈ −60 s. Due to the increasing undercooling of the solution, the suspension becomes supersaturated, reaching a level slightly below σ = 0.12, see Figure 12 (right). The increase of the supersaturation is retarded by the transfer of KDP from the dissolved to the crystalline state by growth and nucleation. However, the amount of nucleated particles is rather low, see, for example, Figure 11. In terms of the mass fraction of nucleated crystals, this is quantitatively

Figure 16. Mean geometry of the crystals that are recognized as single crystals in experiment 1.

through the microscope’s cuvette in which the suspension is imaged, see section 4.1 and particularly Figure 9. The hardness of the crystals is relatively high51 so that disintegration of single crystals due to mechanical stressing in the peristaltic pump was not observable at the chosen flow rate of 250 mL/min. In order to track the composition of the fluid phase, clear solution is continuously withdrawn through a frit, heated to a reference temperature (here 40 °C) and subsequently fed to a densitometer (model Mettler Toldeo DE 40 density meter) in which the density of the clear mother liquor is measured. After the measurement pipe, the solution is fed back to the vessel. Since the fluid phase consists of water and KDP only, which is relatively well soluble in water, the current composition can be determined from the density. Before the start of the experiment, the densitometer was fed with clear solution of known composition in order obtain a calibration curve with which the time series of the density evolution, recorded during the course of crystallization, can be converted to a composition time series. With the help of the temperature that is measured in the vessel, the equilibrium composition is calculated by the empirical solubility curve: wKDP,sat(T ) = (5.5843 × 10−5)T 2 − 0.028159T + 3.6832 (21)

where wKDP,sat is the ratio of the maximally soluble mass of KDP per unit mass of water and thus has the unit [kgKDP/kgwater]. The parameters of this empirical approach were calculated based on the values given by Mullin51 and supported by our own measures. The supersaturation is defined by wKDP σ= −1 wKDP,sat (22) The crystal suspension was continuously photographed at a rate of 10 Hz in order to track the state of the crystalline phase.

Figure 17. State of the continuous phase for experiment 2. Due to the moderate cooling rate (left), the supersaturation (right) is lower than in experiment 1, see Figure 12. 964



Article

Figure 18. Total number of particles and number of particles classified as single crystals (left) and the number fraction of single crystals (right), measured within a sampling interval of 25 s for experiment 2.

Figure 19. Evaluation of the amount of nuclei produced over the course of experiment 2. (left) Particles that are smaller than the average sphere equivalent diameter of the single crystals minus 200 μm are classified as nuclei. (right) Mass ratio of the nuclei.

illustrated in Figure 13 (right). It can be observed that the mass of the nuclei remains well below 2%. However, this should be seen only as a rough estimate obtained by the following pragmatic approach. The number-based average sphere equivalent diameter of that part of the (seed) population that has been identified as single crystals is taken, see Figure 13 (left). Initially, the average size is about 210 μm; therefore particles whose size is lower than the average diameter minus 200 μm are classified as nuclei. The mass of a particle is assumed to be proportional to the power of 3/2 of the measured particle area for all particles regardless of particle shape, that is, the estimated mass-ratio between nuclei and all particles is obtained as n

∑i =nuc1 Ai 3/2 mnuc = n mparticles ∑i =particles Ai 3/2 1

(23)

where nnuc and nparticle are the number of nucleated and the number of all particles, respectively. Though this approach is rather simple, it can for instance be seen that the slope of the nuclei mass fraction reaches its maximum at t ≈ 1000 s, see Figure 13 (right), which corresponds to the maximum supersaturation level depicted in Figure 12 (right). The result of the application of the image processing routines discussed in detail in section 4 and the subsequent identification of

Figure 20. Example of a processed image of experiment 2 at t ≈ 4000 s. 965



Article

Figure 21. Evolution of the shape distribution of experiment 2 using a kernel density estimator.

the time direction. In Figure 15, besides the number density, the trajectory of the mean state of the seed crystals is drawn as a thick black line. This mean is calculated as follows: First, the maximum, hmax, of the estimated number density is determined. Then the first moment of the density function in the window of h1,max ± 100 μm, h2,max ± 100 μm is calculated and taken as the mean of the seed population. At t = 0, the population (and seed) mean is h1 ≈ 80 μm and h2 ≈ 110 μm. The seed crystals spread around this value where the standard deviation in h1- and h2-directions is about 20 μm in all directions. Keep in mind that not all crystals are used to estimate the h-distribution but only those that have been identified as single crystals and thus are eligible to be fed to the shape estimation scheme. For this reason, the unquestionable presence of aggregation producing considerably larger particles than can be seen here is not

the geometrical state vector using the routines developed in section 3 are depicted in Figure 15; The contours of the normalized number density function in h-space estimated using a kernel density estimator, see Borchert and Sundmacher19 and references therein, are shown at eight different instances starting at t = 0. The bandwidth of the Gaussian kernel was set to ⎡ 20 μm 2 0 0 ⎤ ⎢ ⎥ 2 ⎥ H=⎢ 0 20 μm 0 ⎢ ⎥ ⎢⎣ 0 0 60 s2 ⎥⎦

(24)

that is, the influence domain of a single measurement (one particle) in terms of the variance of the Gaussian kernel function is 20 μm in each direction of the h-space and 60 s in 966



Article

Gi = kg, iσ gi

observable. On the right side of the distribution snapshots, the geometry of the population mean of the single crystals is drawn. As time progresses and supersaturation is built up by undercooling (Figure 12), the population moves in state space to the upper right. Clearly, the height and width of the population distribution increases, which may be attributed to growth dispersion. Besides the measurement of the mean crystal shape evolution, this experiment was used to analyze the computational performance of the entire shape estimation procedure. The entire computation time for the processing of all 25000 frames of this experiment was about 8400 s on a standard PC (computation on a single core, 3.00 GHz, 4.00 GB RAM, Windows 7, MATLAB 2010B). Here, 2400 s were spent on the basic image processing and 550 s on the object quantification and classification. In total, about 103000 objects were identified as single crystals and subsequently passed to the shape estimation procedure, which took about 5450 s. These times correspond to an average performance of approximately 3.0 frames per second, which could be analyzed in real time using the proposed methods at their current status. 4.3.2. KDP Growth at Low Supersaturation (Experiment 2). In the preceding experiment 1, we have seen that at high supersaturation the crystals keep their compact shape. Yang et al.52 report that crystals grown at low supersaturation assume a more elongated shape. Therefore, we performed experiment 2 under practically equal conditions to experiment 1, see Table 1, except for the cooling rate, which at −5 K/h was set to a value of a quarter of that chosen in experiment 1. The resulting temperature and supersaturation profiles are shown in Figure 17. The supersaturation remains below σ = 0.08. Figure 18 (left) depicts the overall number of particles and the number of detected single crystals that are recorded within a sampling interval of 25 s. It can be seen that the number of particles that are observed remains over the course of the whole experiment around 7000 and the number of single crystals falls from around 2000 to about 200 at the end of the batch. In terms of the number fraction, this is a decrease from about 25% to below 5%. This is mainly due to crystal aggregation, which is more noticeable compared with experiment 1, see Figure 14, owing to the longer processing time. In the interval 3000−3400 s, the total number of particles drastically increases and the fraction of single crystals significantly drops. At this time, problems with the feeding of the microscope occurred, that is, smaller particles clogged the supply piping, which was dispensed during that interval. After 4000 s, similar problems occurred, which finally led to the abortion of the experiment. As in the previous case, also in this experiment, nucleation was relatively low as can be seen in Figure 19. The drastic increase of the nucleated crystals after about 3000 s is attributed to sampling problems due to clogging, which caused the preferential sucking of smaller particles. The lower supersaturation environment caused the seed crystals to grow to a more elongated shape, which means that the values for h2 (pyramidal faces) become significantly larger than those for the prismatic h1 faces as depicted in an image of the suspension in Figure 20. In Figure 21, this can be seen on the basis of the h-distribution evolution, but even clearer in the time series of the mean state shown in Figure 22 (compare with Figure 16 of experiment 1).

(25)

the measured evolution of the mean crystal shape h̅j can be compared with the mean crystal size following from the growth law together with supersaturation measurements: h̅ mod(t ) = h̅ 0 + g = (g1 , g2)T

∫0

t

k g σ g dt ,

k g = (k1 , k 2)T , (26)

The deviation between measured and modeled mean size evolution is quantified by the objectives

Figure 22. Mean geometry of the crystals that have been classified as single crystals in experiment 2.

Figure 23. Mean geometry of the crystals that have been classified as single crystals in experiments 1 and 2 with model fits.

5. ESTIMATION OF GROWTH RATES The quantitative observation of a crystallization process is of value in its own right. But even more desirable is the extraction of kinetic data from an observed process. Since only crystal growth has been taken into account in our analysis, we aim at estimating the underlying growth rates as a function of supersaturation. If a growth model of the power-law form is assumed,

Figure 24. Comparison of the KDP growth kinetics measured by Gunawan et al.54 and determined in this work. 967



Article

Figure 25. Temperature profile of experiment 1 (left) and the resulting supersaturation profile using the estimated growth parameters (right).

Figure 26. Temperature profile of experiment 2 (left) and the resulting supersaturation profile using the estimated growth parameters (right).

Figure 27. Temperature profile of experiment 1 (left) and the resulting supersaturation profile using the parameters from the literature.54 The supersaturation can obviously not be reproduced well for the conditions of experiment 1 (right). ns

e1 =

kg,2 = 2.07 × 10−6 m/s,

)2 ̅ ∑ (h1,̅ j − h1,mod (27a)

j=1

)2 ̅ ∑ (h2,̅ j − h2,mod (27b)

j=1

The parameters minimizing both objectives were found to be: kg,1 = 643.2 × 10−6 m/s,

g1 = 3.89

(28b)

The evolution of the measured mean crystal size together with the modeled mean size evolution is shown in Figure 23. It can be seen that obtained kinetic parameters reflect the measured mean size evolution well. The identified growth rates are comparable to previous studies on the single crystal level. See, for example, Alexandru53 and references therein who found for σ = 10% a value of approximately 0.1 μm/s for the pyramidal {101} faces. In Figure 24, the obtained growth laws obtained are depicted together with the growth laws reported

ns

e2 =

g2 = 1.23

(28a) 968



Article

by Gunawan et al.,54 which were found to predict significantly higher growth rates. It should be noted that for the determination of the growth rates only the measured supersaturation profile and the measured mean crystal size evolution is used and a full process model is not fitted. That is in particular important since the mass balance of the continuous phase was not used in the estimation procedure. Therefore, the simulation of a full population balance model (for details, see Borchert and Sundmacher49) coupled to the mass balance of the continuous phase constitutes a test whether the obtained growth kinetics can reproduce the measured supersaturation profile. Since the nucleation was found to be rather low (see section 4.3), only growth has been accounted for in the population balance model. The simulated supersaturation profiles are depicted in Figures 25 and 26 for experiments 1 and 2, respectively. It can be seen that the measured supersaturation profiles can be reproduced with high accuracy in both cases, whereas the simulation with the parameters taken from Gunawan et al.54 leads to a supersaturation profile that is significantly lower than the measured one, see Figure 27.

Because of the second observation, the quality of the orientation estimate will be significantly increased, which will result in a significantly better shape estimate.

■

AUTHOR INFORMATION

Corresponding Authors

*Kai Sundmacher, Tel: +49 391 61110 350, Fax: +40 391 6110 523. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Kai Sundmacher gratefully acknowledges the support of this work by the German Research Foundation (DFG) under the Grant SU 189/6-1 in the framework of the priority program SPP 1679.

■

NOTATION

Latin

A B B b d e f G g g H h kg kg m N N P P p R r

6. CONCLUSIONS In this article, an approach for the reconstruction of 3D crystal shapes from recorded microscopic images was presented. The methods account for an unknown orientation from which the crystal was photographed, which is typically the case when batch crystallization processes are monitored by video microscopy. It was demonstrated that the presented methods are robust even to blurred and pixelated images of asymmetric crystals. Although the estimation error distribution is significantly higher compared with the analysis of images with perfectly shaped crystals at infinite resolution, the mean crystal shape can be reproduced well. Due to that and a comparably low computational effort, the presented methods open the way for online process monitoring and thereby process control with regard to crystal shape. To demonstrate the applicability of the proposed methods, batch cooling crystallizations of KDP dissolved in water were monitored, and the supersaturationdependent and face-specific growth rates of the KDP crystals were determined, which could further be used for process modeling and optimization. Of course, KDP crystals have a rather simple geometry, which is composed only of prismatic {100} and pyramidal {101} faces. One can easily think of crystals exhibiting a more complex shape, for example, more different face types can appear on the crystal surface. It is expected that a larger database will be required in such cases. This is because a higher variety of qualitatively different geometries will be possible,42,55 and the features of individual faces are likely to be less prominent, if present at all, on the projection boundary. In such cases, the incorporation of features within the crystal projection in the estimation scheme, for example, as demonstrated by Ma et al.35 or Larsen et al.,36 may help as well in improving the estimation performance. Another challenging example is the shape estimation of needle or plate-like crystals as they are frequently seen in industrial crystallization processes. Due to the small crystal thickness, the estimate for the crystal orientation will be highly sensitive to the finite image resolution resulting in an erroneous shape estimate. An alternative in such cases is to observe the process from two different directions, for example, as demonstrated by Kempkes et al.17 and Schorsch et al.,39 and to equip the used database with a second set of descriptors accounting for the additional crystal projection.

r S s T t v w x x,

projection area matrix of boundary pixels boundary of the crystal projection coordinates of a boundary pixel vector of Fourier descriptors error measure frame rate growth rate vector of growth rate exponents growth rate exponent bandwidth matrix geometrical state vector vector of growth rate constants growth rate constant mass matrix of face normals number of sampled boundary pixels projection matrix crystal projection vector of the projected crystal vertices Fourier descriptor 3D space coordinate/vector of centroid to boundary distances centroid to boundary distance crystal polyhedron scaling factor temperature time coordinate of crystal vertex concentration 2D coordinate of projection plane y coordinates of the projection plane

Greek

α θ λ σ ψ ψ, ϑ, ϕ

■

weighting factor centroid to boundary angle weighting factor supersaturation vector of Euler angles Euler angles

REFERENCES

(1) Zhang, Y.; Sizemore, J. P.; Doherty, M. F. Shape evolution of 3dimensional faceted crystals. AIChE J. 2006, 52, 1906−1915.

969



Article

(2) Gonzalez, R. C.; Woods, R. E.: Digital Image Processing; Pearson Prentice Hall: Upper Saddle River, NJ, 2008. (3) Midgley, P. A.; Ward, E. P. W.; Hungria, A. B.; Thomas, J. M. Nanotomography in the chemical, biological and materials sciences. Chem. Soc. Rev. 2007, 36, 1477−1494. (4) Buzug, T. M. Computed Tomography − From Photon Statistics to Modern Cone-Beam CT; Springer: Berlin, 2008. (5) Merrifield, D. R.; Ramachandran, V.; Roberts, K. J.; Armour, W.; Axford, D.; Basham, M.; Connolly, T.; Evans, G.; McAuley, K. E.; Owen, R. L.; Sandy, J. A novel technique combining high-resolution synchrotron x-ray microtomography and x-ray diffraction for characterization of micro particulates. Meas. Sci. Technol. 2011, 22, No. 115703. (6) Jerram, D. A.; Mock, A.; Davis, G. R.; Field, M.; Brown, R. J. 3D crystal size distributions: A case study on quantifying olivine populations in kimberlites. Lithos 2009, 112, 233−235. (7) Higgins, M. D. Measurement of crystal size distributions. Am. Mineral. 2000, 85, 1105−1116. (8) Holden, E. J.; Moss, S.; Russel, J. K.; Dentith, M. C. An image analysis method to determine crystal size distributions of olivine in kimberlite. Comput. Geosci. 2009, 13, 255−268. (9) Singh, M. R.; Chakraborty, J.; Nere, N.; Tung, H.-H.; Bordawekar, S.; Ramkrishna, D. Image-analysis-based method for 3D crystal morphology measurement and polymorph identification using confocal microscopy. Cryst. Growth Des. 2012, 12, 3735−3748. (10) Castro, J. M.; Cashman, K. V.; Manga, M. A technique for measuring 3D crystal-size distributions of prismatic microtites in obsidian. Am. Mineral. 2003, 88, 1230−1240. (11) Conchello, J. A.; Lichtman, J. W. Optical sectioning microscopy. Nat. Methods 2005, 2, 920−931. (12) Wilson, T. Optical sectioning in fluorescence microscopy. J. Microsc. 2011, 242, 111−116. (13) Nagy, Z. K.; Fevotte, G.; Kramer, H.; Simon, L. L. Recent advances in the monitoring, modelling and control of crystallization systems. Chem. Eng. Res. Des. 2013, 91, 1903−1922. (14) Puel, F.; Marchal, P.; Klein, J. Habit transient analysis in industrial crystallization using two dimensional crystal sizing techniques. Chem. Eng. Res. Des. 1997, 75, 193−205. (15) Patience, D. B.; Rawlings, J. B. Particle-shape monitoring and control in crystallization processes. AIChE J. 2001, 47, 2125−2130. (16) Eggers, J.; Kempkes, M.; Mazzotti, M. Measurement of size and shape distributions of particles through image analysis. Chem. Eng. Sci. 2008, 63, 5513−5521. (17) Kempkes, M.; Vetter, T.; Mazzotti, M. Measurement of 3D particle size distributions by stereoscopic imaging. Chem. Eng. Sci. 2010, 65, 1362−1373. (18) Ferreira, A.; Faria, N.; Rocha, F.; Teixeira, J. A. Using an online image analysis technique to characterize sucrose crystal morphology during a crystallization run. Ind. Eng. Chem. Res. 2011, 50, 6990−7002. (19) Borchert, C.; Sundmacher, K. Crystal Aggregation in a Flow Tube: Image-Based Observation. Chem. Eng. Technol. 2011, 34, 545− 556. (20) Li, R. F.; Penchev, R.; Ramachandran, V.; Roberts, K. J.; Wang, X. Z.; Tweedie, R. J.; Prior, A.; Gerritsen, J. W.; Hugen, F. M. Particle shape characterization via image analysis: From laboratory studies to in-process measurements using an in situ particle viewer system. Org. Process Res. Dev. 2008, 12, 837−849. (21) Li, R. F.; Thomson, G. B.; White, G.; Wang, X. Z.; De Anda, J. C.; Roberts, K. J. Integration of crystal morphology modeling and online shape measurement. AIChE J. 2006, 52, 2297−2305. (22) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. An algorithm for analyzing noisy, in situ images of high-aspect-ration crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61, 5236−5248. (23) Simon, L. L.; Oucherif, K. A.; Nagy, Z. K.; Hungerbuhler, K. Histogram matching, hypothesis testing, and statistical control-chartassisted nucleation detection using bulk video imaging for optimal switching between nucleation and seed conditioning steps. Ind. Eng. Chem. Res. 2010, 49, 9932−9944. (24) Simon, L. L.; Merz, T.; Dubuis, S.; Lieb, A.; Hungerbuhler, K. In-situ monitoring of pharmaceutical and specialty chemicals

crystallization processes using endoscopy-stroboscopy and multivariate image analysis. Chem. Eng. Res. Des. 2012, 90, 1847−1855. (25) Monnier, G.; Fevotte, G.; Hoff, C.; Klein, J. P. An advanced calorimetric approach for population balance modelling in batch crystallization processes. Thermochim. Acta 1996, 289, 327−341. (26) Patchigolla, K.; Wilkinson, D. Crystal shape characterization of dry samples using microscopic and dynamic image analysis. Part. Part. Syst. Charact. 2009, 26, 171−178. (27) Alander, E. M.; Uusi-Penttilä, M. S.; Rasmuson, Ǻ . C. Characterization of paracetamol agglomerates by image analysis and strength measurement. Powder Technol. 2003, 130, 298−306. (28) Alander, E. M.; Uusi-Penttilä, M. S.; Rasmuson, Ǻ . C. Agglomeration of paracetamol during crystallization in pure and mixed solvents. Ind. Eng. Chem. Res. 2004, 43, 629−637. (29) De Anda, J. C.; Wang, X. Z.; Roberts, K. J. Multi-scale segmentation image analysis for the in-process monitoring of particle shape with batch crystallizers. Chem. Eng. Sci. 2005, 60, 1053−1065. (30) 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, 1406−1414. (31) De Anda, J. C.; Wang, X. Z.; Lai, X.; Roberts, K. J. Classifying organic crystals via in-process image analysis and the use of monitoring charts to follow polymorphic and morphological changes. J. Process Control 2005, 15, 785−797. (32) Wang, X. Z.; De Anda, J. C.; Roberts, K. J. Real-time measurement of the growth rates of individual crystal facets using imaging and image analysis − A feasibility study on needle-shaped crystals of L-glutamic acid. Chem. Eng. Res. Des. 2007, 85, 921−927. (33) 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, 1173−1184. (34) Ma, C. Y.; Wang, X. Z. Model identification of crystal facet growth kinetics in morphological population balance modeling of Lglutamic acid crystallization and experimental validation. Chem. Eng. Sci. 2012, 70, 22−30. (35) Ma, C. Y.; Wan, J.; Wang, X. Z. Faceted growth rate estimation of potash alum crystals grown from solution in a hot-stage reactor. Powder Technol. 2012, 227, 96−103. (36) 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, 1430−1441. (37) Larsen, P. A. Computer Vision and Statistical Estimation Tools for In Situ, Imaging-based Monitoring of Particulate Populations, PhD thesis, University of WisconsinMadison, 2007 (38) Eggers, J.; Kempkes, M.; Cornel, J.; Mazzotti, M.; Koschinski, I.; Verdurand, E. Monitoring size and shape during cooling crystallization of ascorbic acid. Chem. Eng. Sci. 2009, 64, 163−171. (39) Schorsch, S.; Vetter, T.; Mazzotti, M. Measuring multidimensional particle size distributions during crystallization. Chem. Eng. Sci. 2012, 77, 130−142. (40) Burger, W.; Burge, M. J. Digital Image Processing − An Algorithmic Introduction Using Java; Springer: New York, 2008. (41) Ziegler, G. M. Lectures on Polytopes; Springer: Berlin, 2006. (42) Borchert, C.; Sundmacher, K. Efficient formulation of crystal shape evolution equations. Chem. Eng. Sci. 2012, 84, 85−99. (43) Bronstein, I. N.; Semendjajew, K. A.; Musiol, G.; Mühlig, H. Taschenbuch der Mathematik; H. Deutsch: Frankfurt, 2001 (44) Eggers, J. Modeling and Monitoring of Shape Evolution of Particles in Batch Crystallization Processes, PhD thesis, ETH Zürich, 2008. (45) Loncaric, S. A survey of shape analysis techniques. Pattern Recognit. 1998, 31, 983−1001. (46) Zahn, C. T.; Roskies, R. Z. Fourier descriptors for plane closed curves. IEEE Trans. Comput. 1972, 21, 269−281. (47) Granlund, G. H. Fourier processing for hand print character recognition. IEEE Trans. Comput. 1972, 21, 195−201. (48) Persoon, E.; Fu, K. S. Shape discrimination using Fourier descriptors. IEEE Trans. Syst., Man, Cybern. 1977, 7, 170−179. 970



Article

(49) Borchert, C.; Sundmacher, K. Morphology evolution of crystal populations: Modeling and observation analysis. Chem. Eng. Sci. 2012, 70, 87−98. (50) Kempkes, M. Monitoring of Particle Size and Shape in Crystallization Processes; PhD thesis, ETH Zürich, 2009. (51) Mullin, J. A. Crystallization; Butterworth-Heinemann: Oxford, 2001. (52) Yang, G.; Kubota, N.; Sha, Z.; Louhi-Kultanen, M.; Wang, J. Crystal shape control by manipulating supersaturation in batch cooling crystallization. Cryst. Growth Des. 2006, 6, 2799−2803. (53) Alexandru, H. V. KDP kinetics and dislocation efficiency of growth. J. Cryst. Growth 1999, 205, 215−222. (54) Gunawan, R.; Ma, D. L.; Fujiwara, M.; Braatz, R. D. Identification of kinetic parameters in multidimensional crystallization processes. Int. J. Mod. Phys. B 2002, 16, 367−374. (55) Reinhold, A.; Briesen, H. Convex geometry for the morphological modeling and characterization of crystal shapes. Part. Part. Syst. Charact. 2011, 28, 37−56.

971


Image-Based in Situ Identification of Face Specific Crystal Growth

Recommend Documents