Shape Identification of Primary Particles in Potash Alum Aggregates

Article pubs.acs.org/crystal

Shape Identification of Primary Particles in Potash Alum Aggregates Using Three-Dimensional Tomography Data Tijana Kovačević,† Jonathan Schock,‡ Franz Pfeiffer,‡ and Heiko Briesen*,† †

Lehrstuhl für Systemverfahrenstechnik, Wissenschaftszentrum Weihenstephan, Technische Universität München, 85354 Freising, Germany ‡ Lehrstuhl für Biomedizinische Physik, Physik-Department & Institut für Medizintechnik, Technische Universität München, 85748 Garching, Germany

ABSTRACT: The degree of agglomeration and the aggregate shape influence the quality of crystalline products and the ease of downstream processing. Studying the shape of primary particles in an aggregate can lead to a better understanding of the underlying aggregation mechanism. We present an automatic image processing procedure for identifying the shape, size, and position of each primary particle in microcomputed tomography (μCT) images of potash alum aggregates. Splitting an aggregate into primary particles is based on recombining watershed-transform regions, where concavity points are considered as indicators of correct segmentation. The shape identification algorithm uses the Hough transform to identify visible face normals and matches them to the set of face normals defined by a crystal model. In principle, the algorithm is applicable to other crystalline compounds provided that sufficient symmetry is present to determine the shape of a primary particle from its visible part.

1. INTRODUCTION Crystal shape influences the quality of the final product, such as the solubility or bioavailability of pharmaceuticals, as well as the ease of downstream processing required in manufacturing.1 Significant effort is being invested in designing processes that yield a product with the desired properties,2 where agglomeration is usually undesirable. In recent years, sophisticated methods for both growth modeling3−5 and detailed characterization6−13 of the size and shape of single crystals have been developed. Algorithms for processing images obtained by in situ probes even enable segmentation and size measurement of elongated particles in the presence of some degree of particle overlap.6,14,15 However, in situ imaging probes cannot be used to study aggregation as it is not possible to determine whether the overlapping crystals represent aggregates or merely overlap when projected to the image plane. Other imaging techniques enable studying aggregation by measuring particle size and classifying crystals based on shape descriptors.16−20 However, apart from diffraction pattern studies aimed at determining primary particle orientation in calcite,21,22 no attempts to fully characterize the size and shape of primary particles that constitute the aggregate are known to the authors. Aggregation modeling involves considering the aggregation rate, which is defined by the rate of collisions and their efficiency. The efficiency of collisions depends on the strength of the neck formed between the particles and the hydrodynamic forces.23−26 Such models typically assume spherical particles and, apart from one known case,27 they do not consider the influence of particle shape on aggregation efficiency, even when modeling aggregate morphology.25,28 The ability to exper© 2016 American Chemical Society

imentally determine the orientation between primary particles in aggregates can lead to a better understanding of aggregation efficiency. Furthermore, it can enable the study of other potentially occurring mechanisms that may lead to joint particles, such as twinning. The aim of this work is to present an imaging and image processing tool for such investigations. Herein, we extend our method for extracting full threedimensional (3D) crystal shape from ex situ μCT images10 to identifying the shape and size of each primary particle in aggregates. In order to perform the shape identification, it is first necessary to separate an aggregate into primary particles. The developed algorithm must assume that primary particles have different shapes and cannot be based on searching for spherical structures, as has been successfully implemented for maltodextrin aggregates measured by μCT.29 This problem statement is similar to that of splitting biological cell clusters into individual cells or nuclei, which, among others, has an important role in cancer research.30 Many employed methods are based on watershed transform,31 where oversegmentation is resolved by finding region markers or by merging and splitting the obtained regions.32−34 Another group of methods uses concavity information and splits clusters in two-dimensional (2D) images by connecting concavity points.35−39 The two approaches were coupled by using concavity information to extract watershed markers.40,41 Received: December 22, 2015 Revised: March 10, 2016 Published: March 15, 2016 2685

DOI: 10.1021/acs.cgd.5b01806 Cryst. Growth Des. 2016, 16, 2685−2699

Crystal Growth & Design

Article

crystal geometry and image processing functionality previously developed by the Briesen group,10,42,44 and the cddlib library.45 2.2. Distance Transform-Based Watershed Segmentation. Distance transform-based watershed segmentation is a classical image processing tool that represents the foundation for the developed algorithm for primary particle identification. An overview can be found in the work of Vincent and Soille,31 whereas the book by Gonzalez et al.46 gives Matlab interface details. The name originates from the analogy with a topographic surface where watersheds separate catchment basins containing water. In image processing, the topographic height is defined by the grayscale value of a pixel (voxel). Distance transform can be used to assign to each pixel (voxel) of a binary image a grayscale value corresponding to its distance from the background. The resulting values are inverted to represent depth, and background voxels are assigned a value − ∞ to ensure separation at the boundary. The obtained image is used as input for the watershed transform. This procedure results in a labeled image of segmented objects that are often, but not always, separated near concavity points.

In this work, we combine some of the techniques proposed for cell cluster splitting and adapt them for finding the primary particles in 3D images of crystal aggregates. Furthermore, we modify the single particle shape identification algorithm10 to handle incomplete primary particles. The method is tested on a μCT image of potash alum aggregates as well as on several simulated data sets. This article is structured as follows. Section 2 introduces basic concepts such as the used crystal representation and watershed transform. The methods used for obtaining aggregates and their images, including crystallization, μCT setup, and simulation, are presented in Section 3. Section 4 describes the procedure for segmenting an aggregate into primary particles. The algorithm for identifying the shape of each primary particle is described in Section 5. The results are demonstrated in Section 6. Section 7 concludes the work.

2. BASIC CONCEPTS 2.1. Crystal Representations. A crystal can be described by the matrix A holding nH crystal face normals ai, ∥ai∥ = 1: A = [a1, ..., a nH]T

3. OBTAINING AGGREGATE IMAGES 3.1. Simulated Data. 3.1.1. Model Aggregates. In order to assess the algorithm itself, we used simulated 3D images. This allows full control of aggregate size, degree of primary particle overlap, and image quality. All three are important in order to understand the properties and limitations of the proposed method. We simulated 100 potash alum aggregates consisting of purely octahedron-shaped primary particles, using the procedure described below. We obtained 58 aggregates consisting of two primary particles, 13 containing three primary particles, 21 containing four, and 8 aggregates containing five primary particles. About half of the aggregates contained two primary particles as this was the most common configuration in the observed real data. Primary particles were generated from a normal distribution with the dimensionless mean face distance μh = 250 and standard deviation σh = 30. In the first step, two particles are generated and placed at the origin. One of the particles is moved in a randomly chosen direction, where a bisection algorithm is used to find a position in which the particles lightly touch. In the case where more than two primary particles should be generated, new primary particles are added successively using the same procedure so that they lightly touch the existing aggregate. In order to investigate the influence of primary particle overlap on the algorithm properties, each primary particle is grown according to

(1)

and the vector h of their distances hi to some chosen origin. The crystal polyhedron P is then described by points x such that P = {x | A·x ≤ h}

(2)

This is known as H-representation.42 Alternatively, the polyhedron P can be described by its nV vertices vi, which is called V-representation. As crystals grow symmetrically, a dimensional reduction can be achieved if the origin is centered in the crystal. Then, only one face distance hi for each of the nC face groups is considered to obtain the vector hC. The fulldimensional vector h is related to the lower dimensional hC by the group mapping matrix MhC→h:42 h = M hC → h · h C

(3)

This is referred to as HC-representation. A reverse mapping is obtained using a pseudoinverse matrix denoted by +:42

h C = M +hC → h ·h

(4)

In this work, we use potash alum crystals which are typically modeled considering up to three face groups, as shown in Figure 1. Unless stated otherwise, computations were performed in Matlab 2014a and Matlab 2015b,43 using the

hi = hi + g ·t

(5)

where hi is the distance of the ith face from the crystal origin, g = 1 is the growth rate, and t simulates the dimensionless time that passed since the moment the aggregate was formed. We obtained three data sets by setting t to 50, 150, and 300 in order to simulate low, middle, and high primary particle overlap, as shown in Figure 2. 3.1.2. Generating Artificial 3D Micrographs. The resolution of artificial images was chosen based on the mean size as ⎡ μ + g ·t ⎤ r=⎢ h ⎣ 20 ⎥⎦

Figure 1. Potash alum crystal model with Miller indices for each face group. 2686

(6) DOI: 10.1021/acs.cgd.5b01806 Cryst. Growth Des. 2016, 16, 2685−2699


Article

Figure 2. Illustration of the influence of growth time on particle overlap. The aggregate was generated using t = 50 (a, low overlap), t = 150 (b, middle overlap), and t = 300 (c, high overlap) with proughness = 0.6. The identified shape of each primary particle is shown as a black outline.

Figure 3. Illustration of different degrees of surface roughness. The aggregate was generated using t = 150. (a) proughness = 0; (b) proughness = 0.2; (c) proughness = 0.4; (d) proughness = 0.6. The identified shape of each primary particle is shown as a black outline.

distance no larger than √2 from the center, resulting in an 18connected neighborhood of the central voxel. The considered values of parameter proughness were 0.2, 0.4, and 0.6, as well as 0, which represents the case where no surface roughness was simulated. Therefore, together with different values of t, 12 data sets, each containing 100 aggregates, were obtained, resulting in 1200 simulated aggregates. The influence of proughness on the image quality can be seen in Figure 3. Note that the uneven appearance of the particle surface in all images, including the proughness = 0 case seen in Figure 3a, does not represent boundary roguhness but is caused by the discretization into voxels and the Marching Cubes rendering algorithm.47 3.2. Experimental Data. 3.2.1. Crystal Preparation and Measurement. Potash alum crystals were grown by seeded cooling crystallization in a batch reactor. A sample of aggregates, chosen by visual inspection, was scattered over a microscope slide so that there is no direct contact between particles. The scattered crystals were covered by an adhesive tape which was then rolled and placed in a 0.65 mL micro test tube. The middle part of the test tube was imaged by μCT. The obtained 3D image, rendered by AVIZO Fire,48 as well as one vertical slice, are shown in Figure 4. The observed crystals exhibit a pronounced octahedron face group with Miller indices

so that a similar image quality is obtained for all data sets. For each aggregate, a 3D image matrix of length lj in dimension j was created, where ⎡ ⌈vj ,max ⌉ − ⌊vj ,min ⌋ ⎤ ⎥ lj = ⎢ r ⎥ ⎢

(7)

Here, vj,min and vj,max represent the aggregate vertex coordinates with the minimal and maximal value in the jth dimension, respectively. Translation and scaling with r are performed in order to map the voxel with coordinates (x, y, z) to a point p in the aggregate coordinate system: p = ((x , y , z) − 0.5)·r + ⌊(vx ,min , vy ,min , vz ,min)⌋

(8)

A voxel (x, y, z) is assigned the value 1 (foreground) if the point p is inside one of the primary particles. 3.1.3. Simulating Surface Roughness. The obtained images simulate perfect particles measured by a perfect μCT machine, using a given resolution. However, experimental data show rough surfaces due to both imperfect particles and measurement inaccuracy. We simulate this roughness by inverting the values of some of the voxels close to the aggregate surface. The degree of roughness is controlled by the parameter proughness. Voxels that may be inverted are found in the space between a larger and a smaller version of the aggregate, obtained by reducing and increasing the length of each primary particle by 10%. This space is divided into the space above the surface, which, in the beginning, is populated by the background voxels, and the space below the surface, which is populated by the foreground voxels. Each voxel in the space above the surface, if connected to at least nine foreground voxels, has a probability proughness of becoming a foreground voxel. Similarly, each voxel in the space below the surface, if connected to fewer than nine foreground voxels, has a probability proughness of becoming a background voxel. This ensures a rough boundary where the newly added foreground voxels are not isolated and the removed voxels do not create holes inside the crystal volume. Finally, morphological opening is performed to smooth the obtained surface. The structuring element is a 3D matrix of length 3 with values 1 in the center and for each voxel with a

Figure 4. Experimental data. (a) 3D μCT image; (b) One vertical slice. 2687



Article

{111}, whereas the remaining two face groups with respect to the model shown in Figure 1 correspond to small faces. The acquisition of the sample was carried out with a ZEISS X-Radia Versa XRM-500 μCT machine with a tungsten target. The measurement was performed at a peak voltage in the X-ray tube of 60 kV and a power of 5 W. A total of 1601 projections were taken with an exposure time of 6 s, each with an intrinsic detector pixel binning of 2 by 2. The source-to-sample distance was 50 mm, and the detector-to-sample distance was 200 mm, resulting in an effective voxel edge length of 13.60 μm using the 0.39× magnification objective of the system. The CT reconstruction was performed with the provided proprietary software XRM Reconstructor, v 10.7.3245. 3.2.2. Isolating Aggregates. Aggregates were isolated from the binarized 3D μCT image of the entire sample by labeling connected components using the software MAVI.49 Binarization threshold was chosen manually by inspecting the results on different slices, which is enabled in the graphical user interface of the MAVI software. Due to transport and the simple preparation procedure, it is still possible that two or more aggregates lightly touch and are therefore considered one connected component. Lightly touching aggregates may pose problems for the chosen concavity-finding algorithm and must be split into aggregates that contain overlapping particles. The aggregate image is eroded using a cube structuring element of length 15, which ensures separation but also shrinks the separated aggregates. The obtained image is labeled to check whether more than one connected component is obtained. If so, foreground labels are used as a mask for imposing minima on the inverted distance-transform image of the original crystal. Performing watershed separtaion on this image results in the desired separation, similar to that of the eroded image. The difference is that the previously eroded voxels will be assigned to one of the obtained labels. Watershed voxels, as well as foreground voxels that were assigned to background by watershed transform, are reassiged to the nearest connected component, as will be explained in Section 4.3.1. The resulting objects represent aggregates with overlapping particles and are stored as separate images. Particles consisting of fewer than 10 000 voxels are discarded from consideration due to insufficient resolution. 3.3. Preprocessing. The procedure above results in one image per aggregate for both real experimental data and simulated particles. Before proceeding with segmentation into primary particles, each aggregate image is filtered. First, the image is eroded with a 3D structuring element of length 3 with a value 1 in the matrix center and at each voxel with distance 1 from the center (6-neighborhood). All connected components obtained thereafter with the exception of the largest are then deleted. Dilation with the same structuring element and a holefilling operation follow. The results of the filtering operation for one crystal are presented in Figure 5, where it can be seen that surface irregularities present in Figure 5a were successfully removed in Figure 5b. After filtering, the surface and volume point clouds, which will represent the basis for the following algorithm, are extracted. The volume point cloud consists of points whose coordinates represent the indices of the foreground voxels in the image matrix. The surface point cloud contains the voxels on the aggregate surface. These voxels are obtained using the Matlab function bwperim which considers a foreground voxel to be on the object surface if it is connected to at least one background voxel. Furthermore, a connected component search

Figure 5. Surface irregularities present in (a, unfiltered) were removed by filtering to obtain (b, filtered).

is performed on the obtained surface boundary, and all voxels that are not connected to the largest connected component are deleted. This ensures there are no voxels inside the crystal volume that are considered to belong to the surface due to surface artifacts.

4. SEGMENTATION INTO PRIMARY PARTICLES Before shape identification, it is necessary to split each aggregate into primary particles. An overview of the developed algorithm is shown in Figure 6, and the application of each step

Figure 6. Algorithm for splitting an aggregate into primary particles.

to one crystal is illustrated in Figure 7. The algorithm is based on seeded watershed segmentation and region recombination using concavity points as indicators of primary particle intersection. 4.1. Detecting Concavity Points. Concavity points are detected using a method similar to the 2D approach seen in the works of Fernandez et al.40 and Indhumathi et al.38 A cubic mask C is centered at each aggregate surface point p, and the concavity value cp is calculated as the ratio between the number of foreground points inside the mask, Nforeground mask, and the total number of voxels in the mask, Nmask:

cp =

Nforegroundmask Nmask

(9)

This is illustrated in Figure 8a. In a system with infinite resolution, a cp value larger than 0.5 would indicate a point that belongs to a concave region. However, due to discretization, 2688



Article

Figure 7. Separation of an experimentally measured aggregate into primary particles. Parts (a) and (c) show the watershed separation result before and after the concavity expansion (b) was performed, respectively. Concavity points are marked in red and the associated concavity value computation masks in black. The following images show the watershed regions obtained after reassigning watershed voxels (d), after merging small regions (e) and the first (f) and second (g) iterations of recombining the remaining large watershed regions to obtain the final separation into primary particles. Part (h) shows the identified primary particle shape and position.

Figure 8. Concavity expansion illustrated on one concavity point on a zoomed-in crystal surface. (a) A concavity computation mask (black), concavity point pc (red), and the middle point pb of the background part of the mask (blue). (b) Concavity expansion results in a thin hole through the crystal in the direction defined between the blue and the red point in (a).

2.2. However, if the concavity points are not sufficiently pronounced, this will not be the case and undersegmentation will occur. To ensure that the watershed segmentation splits near the concavity points, we expand the concavities in the concavity direction, using a simplified 3D version of the approach proposed by Zhang et al.41 for 2D images. This is illustrated in Figure 8. The concavity direction d is defined between a concavity point pc (red) and the center pb (blue) of the background part of the mask C (black) used for concavity value computation: p − pb d= c pc − pb (12)

this value depends on the mask size and orientation. We use 1.2 ·ct as a concavity threshold, where ct =

(2a + 1)2 ·(a + 1) (2a + 1)3

(10)

represents the concavity value of a point p on a flat surface parallel to the surface of the mask of length 2a + 1. The factor 1.2 is chosen to compensate for boundary roughness and different orientations. Mask length 2a + 1 is chosen according to the aggregate size using

a = [ 3 0.0013Nvoxels ]

(11)

The concavity of the point pc is expanded by changing the voxels p in the direction d from foreground to background. Foreground voxels p are considered to lie on the concavity direction d if

where Nvoxels is the total number of foreground voxels, and therefore represents the aggregate volume. In order to obtain one concavity point per concavity region, a nonmaximum suppression peak search50,51 over concavity values cp is performed. A cube-shaped window with length 2b + 1, where b = ⌊a/2⌋, is centered at each point p of the aggregate surface. The algorithm considers a point p to be a concavity point if its concavity value cp is the maximal concavity value among the surface points in each window that contains the point p. Identified concavity points are deleted if their concavity value cp is smaller than the threshold 1.2 ·ct. Figure 7a illustrates the obtained concavity points (red) for one aggregate as well as the masks used to compute the concavity values (black cubes). 4.2. Concavity Expansion. Segmentation into primary particles can in certain cases be achieved with a distance transform-based watershed algorithm, as explained in Section

p − pb tan(α) ≤ 1

(13)

Here, α is the angle between d and p − pb. In order not to split the object, the switching of a foreground voxel p to background is only performed if the distance from p to pb is smaller than one-quarter of longest such distance found in the direction d: ∥p − pb∥ < 0.25·

max

∥ p − p b tan(α) ≤ 1

∥p − pb∥

(14)

This procedure results in thin holes in the crystal volume at concavity points, as illustrated in Figures 8b and 7b. A distance transform-based watershed transform is performed on the 2689



Article

resulting image to obtain regions that are separated at concavity points. However, the foreground voxels that were deleted in the expansion procedure are not present in this image. Therefore, the obtained watershed regions are used as a mask to impose minima for the distance-transform based watershed separation performed on the original image, similarly to the procedure presented in Section 3.2.2. This ensures that the voxels deleted by the concavity expansion procedure will be present in the final separation result. Figure 7c shows that the obtained regions are indeed separated at concavity points, in contrast to the watershed segmentation obtained without performing concavity expansion, illustrated in Figure 7a. However, it can be seen that the procedure resulted in oversegmentation, which is resolved by region merging under the consideration of concavity points. 4.3. Watershed Region Merging. Before proceeding with region merging, a matrix that defines neighborhood relationships between regions is created. Two regions are considered neighbors if they contain voxels that are within a distance of √3 of each other (26-neighborhood). 4.3.1. Assigning Watershed Voxels to Watershed Regions. The segmentation obtained by the procedure in Section 4.2 includes watershed voxels that separate watershed regions and appear as holes in Figure 7c. Each of these voxels is assigned to the closest neighboring region to close the holes, as seen in Figure 7d. 4.3.2. Merging Small and Large Watershed Regions. Small watershed regions, which are defined as regions containing less than 3% of the total number of voxels in the aggregate, are considered as artifacts caused by surface roughness and irregularities. To prevent such regions from causing issues in concavity-based region recombination described below, they are merged with neighboring regions, similarly to the approach of Adiga and Chaudhuri.32 A small region i is merged with a neighboring region j which will lead to minimizing the concavity of the resulting region, defined as cr = 1 −

Concavity points are considered indicators of correct segmentation and will inhibit regions from merging by deleting neighbor relationships. The influence of each concavity point is determined by considering to which regions the voxels in its neighborhood, defined by the cubic mask C of length 2a + 1, belong. If there are exactly two regions in the vicinity of a concavity point, their merging will be forbidden, as in the case of the light green and yellow regions in Figure 7e. However, if more than two regions are present in the concavity point neighborhood, some region pairs may be suitable for merging, while merging others should be inhibited. To determine their suitability, the regions are considered pairwise, and a concavity value for each pair is defined similarly to eq 9 c p, i , j =

Nmask

(16)

Here, Nvoxels of i in C and Nvoxels of j in C represent the number of voxels of regions i and j within the mask C, respectively, and Nmask is the total number of voxels in mask C. Merging the regions i and j is forbidden if the value cp,i,j is larger than ct. Note that the inhibition threshold is relaxed in comparison to the case where only two regions are found in the concavity point vicinity. In this case, the region merging is always forbidden, which means that the threshold of 1.2·ct, as defined in Section 4.1 for finding concavity points, is implicitly employed. The same region pair will be considered multiple times, as there may be multiple concavity points between the regions, as seen for the blue and green region in Figure 7c. Their merging will be forbiden if there is at least one concavity point where the concavity value is too high. The two iterations of the large-region merging procedure can be seen in Figure 7f and g, resulting in the correct identification of primary particles. The merging is finished if there no longer exists a pair of regions whose merging is allowed. The obtained regions represent incomplete primary particles, and the algorithm proceeds with the identification of their shapes.

Ni + Nj V

Nvoxels of i in C + Nvoxels of j in C

(15)

5. PRIMARY PARTICLE SHAPE IDENTIFICATION In an ideal case, the procedure described in the previous section and summarized in Figure 6 results in the segmentation of the aggregate into primary particles. Below, we present the basic algorithm for shape identification, which is followed by postprocessing steps to treat non-ideal cases. In the section dedicated to the basic algorithm, we first summarize the algorithm for the identification of single crystal shapes,10 which is used if the segmentation procedure resulted in only one primary particle. We then describe the modifications necessary to accommodate incomplete primary particles. 5.1. Basic Algorithm. 5.1.1. Single Crystals. Hough transform52,53 is performed to find face normals, which are denoted by af,i and gathered in the matrix AF. To account for rotation and missing faces, the relationship between AF and face normals described by the crystal model A is then established as

Here, Ni and Nj are the numbers of voxels in regions i and j, respectively, and V is the volume of the convex hull of the merged regions. The absolute value ensures a positive value in cases where V is smaller than the sum of Ni and Nj due to discretization error, which may occur for regions that consist of only a few voxels. The merging is performed iteratively until no more small regions remain. It is possible for multiple small regions to be merged into a large region. The result of this procedure is shown in Figure 7e, in which the violet region present in Figure 7d is merged with the yellow region. 4.3.3. Recombining Large Watershed Regions. Before proceeding with merging watershed regions, neighborhood relationships are redefined. Two regions are now considered neighbors if they contain surface voxels that are at a distance 1 from each other (6-neighborhood). This ensures that the surface of each region obtained by merging is connected. Here the surface of a watershed region is the part of the aggregate surface that belongs to that region. The algorithm iteratively merges the remaining watershed regions to form primary particles. In each iteration, a pair of neighboring regions whose merging would minimize the average region concavity, where region concavity is defined in eq 15, is identified, similarly as in the work of Long et al.34

SD ·A F ≈ (SM ·A) ·RT

(17)

Here, R is a rotation matrix, SD removes extra faces obtained, for example, by breakage, and SM describes which model faces were matched to the identified faces. The concept of matched faces is illustrated in Figure 9, where the crystal model is shown on the left and the considered real particle is illustrated on the right. Here, the octahedron {111} 2690



Article

In this work, we use 40 Hough bins per considered spatial direction in order to achieve faster processing. Other parameters remain as in our previous work.10 5.1.2. Primary Particles in Aggregates. The surface of a primary particle within an aggregate is incomplete as it is defined as the corresponding part of the aggregate surface. This ̃ may lead to additional small peaks bpks,i when searching for crystal faces using the Hough transform. Note that b̃ is the reduced Hough space obtained by considering only the maximal bin count for each spatial search direction, as explained ̃ in our previous work.10 Therefore, a Hough peak of height bpks,i is considered an artifact and is deleted if

Figure 9. Face normals of a crystal model (a) and identified face normals of a real particle (b) are marked by arrows. Blue face normals from the model (a) are matched to the identified face normals of the real crystal (b), whereas the green and orange face normals are not matched.

̃ i < max bpks, j

1 N ∑k =dir1 bk̃ Ndir ̃ j· bpks, 1 N ̃ k ∑k =pks1 bpks, Npks

(19)

Here, Ndir is the number of considered spatial directions, bk̃ is the maximal bin count for the direction k, and Npks is the number of identified peaks, where the identification procedure is presented in our previous work.10 After performing this filtering operation, the rotation matrix R as well as correspondences SM and SD between the crystal model and the remaining face normals are identified similarly to those in the single particle algorithm.10 The slight modification is motivated by the fact that primary particles may exist where only three face normals can be identified by Hough transform to form AF. This would result in stopping the search upon the first rotation matrix candidate Rcand that results in L(Rcand) < 0.02 because all identified face normals would have been matched. As this is not necessarily the best solution, such stopping is no longer possible. Furthermore, if AF consists of three faces, the best rotation matrix candidate is chosen not according to the number of matched faces but as the minimum of the mean error

faces from the model, drawn in blue, are matched to the faces identified by Hough transform, marked by black arrows in the right part of Figure 9. The green and orange model faces are not matched because no corresponding faces were identified by Hough transform in the crystal on the right. These faces are either grown-out or are too small. The procedure for finding matrices R, SD, and SM is described in detail in our previous work10 and summarized here so that the changes required for primary particles can be explained. First, three largest compatible identified faces, whose normals are summarized in AF, are identified. Three faces are compatible if all pairs of their normals form angles between 0.9·αmin and 180° − 0.9·αmin, where αmin is the smallest angle between two face normals in the crystal model A. The chosen faces are then paired with some triplet of model faces from A, filling the corresponding rows of SM. A rotation matrix Rcand is determined using the Markley54 solution to Wahba’s problem55 for the considered three rows of the matrices (note: the equation stated in our previous work10 is missing the squaring symbol) 1 L(R) = ∑ wi∥[SDA F]i − [(SM A)RT]i ∥2 2 i (18)

Lm(R) =

1 2nM

nM

∑ ∥[SDAF]i i=1

− [(SM A)RT]i ∥2

(20)

where nM is the number of matched faces. In the next step of the algorithm, the crystal middle point and the face distances hi, should be identified for each face in the model A. Here, the concept of existing and nonexisting faces is needed. Certain rotated model faces do not exist on the incomplete primary particle surface because they are inside another primary particle. Whether faces exist on the incomplete surface is determined as follows: (1) All matched faces are considered existing. (2) All unmatched {111} faces are considered nonexisting. (3) An unmatched face from a face group other than {111} is declared existing if all neighboring faces from the {111} face group were matched. Otherwise, it is considered nonexisting. These considerations are specific to the observed potash alum crystals, assuming {111} to be the main face group, containing the largest number of matched faces. Existing faces are described using the nE × nH existence mapping matrix SE. This matrix contains a row with the value 1 at position i and zeros otherwise if the model face with rotated normal Rai exists on the primary particle surface. Here, nE is the number of existing faces.

Here, weights wi are chosen to be equal so that ∑i wi = 1. This procedure is repeated for all possible triplets of model face normals. If the rotation error L(Rcand) is smaller than 0.02, the algorithm proceeds by matching the remaining identified face normals from AF to the remaining rotated model face normals in A ·RTcand by considering the angles between them. This defines the remaining rows of SM. The procedure stops if all identified face normals from AF have a match in A. Otherwise, Rcand, which results in the largest number of matched faces from AF, is chosen as the best candidate. The matrix SD is filled by considering which faces in AF have a match in A, and the final matrix R is obtained by solving the problem in eq 18 using full matrices SD ·AF and SM ·A. Hough transform also delivers the distance hi of the matched faces to the particle center. The distance hi of an unmatched face is found as the Hough bin with the highest point density when projected to the plane perpendicular to the face normal. Grown-out faces are identified, and their distances hi are modified such that these faces touch the polytope. Finally, symmetry conditions are applied to the vector of face distances h to obtain a symmetrical crystal denoted by hC, compliant with the given model A. 2691



Article

For existing faces, the distance hi is determined as in the single-particle algorithm. It is not possible to determine the hi distance for nonexisting faces. However, not all hi values are needed to determine the symmetric crystal shape hC and the origin x0. They are obtained similarly to that in the singlecrystal algorithm, with the difference that only existing faces are used: M C,E = InE − (SEM hC → h) ·(SEM hC → h)+

(21)

x 0 = (M C,ESEART )+ M C,ESEh

(22)

SEh′ = SEh + SEART x 0

(23)

and h C = (SEM hC → h)+ SEh′

(24)

If no faces are identified for a face group, a value 0 will result in the corresponding place in hC. Instead, this face group is assigned some large value so that the faces of this group are grown out. A nonredundant set of inequalities that define the polyhedron is then identified,45 so that the remaining, grownout faces can be detected. Their distance is modified to a valid value, as in the single-particle algorithm: hi =

max

1 ≤ i ≤ nH ,1 ≤ j ≤ nV

⟨a i , vj⟩

Figure 10. Shape identification: the basic algorithm followed by postprocessing steps.

(25)

where the mapping according to eq 3 was performed beforehand. The inverse mapping according to eq 4 finalizes the shape identification. 5.2. Postprocessing of Non-Ideal Cases. Due to surface roughness and irregularities, the segmentation into primary particles may be imperfect and result in an extra or a missing region. Furthermore, the primary particles may overlap to a very high degree, leaving little information for shape identification. Therefore, the basic shape-identification algorithm is followed by a set of postprocessing steps, as illustrated in Figure 10 and described below. 5.2.1. Deleting Faulty Regions. If the basic shapeidentification algorithm fails to provide a reasonable shape, the underlying region is considered faulty and is deleted. One reason for this is that not enough face normals are present to find a rotation matrix. Furthermore, the shape is not considered reasonable if some face group other than {111} is considered as the main face group, or the faces of this group are found to be grown-out. An example of a faulty region is illustrated in Figure 11. 5.2.2. Merging Overlapping Particles. Due to an inconvenient concavity point, which is usually caused by a surface irregularity, a primary particle may be split into two regions. In this case, the shape identification would result in two highly overlapping particles as the same shape is approximated from its two parts. Therefore, such regions are merged and the basic shape-identification algorithm is performed on the newly obtained region. If this leads to a region where no shape identification is possible, this faulty region is deleted as in the procedure described above. These cases are detected by considering all identified shapes pairwise and computing the volume of their overlap, Voverlap,i,j. If the degree of overlap with respect to particles i or j is larger than 0.9, regions i and j will be merged. The degree of overlap with respect to the particle i is computed as

Figure 11. Crystal regions before (11(a)) and after (11(b)) deleting the faulty orange region.

oi =

Voverlap, i , j Vi

(26)

where Vi is the volume of the particle i. This is illustrated in Figure 12. 5.2.3. Moving and Resizing the Identified Shape. A primary particle can overlap with another particle to such a degree that it is not possible to identify its size, even though it is possible to determine its orientation, as seen in Figure 13. This is because three face normals are needed to compute the rotation matrix but may not be enough to completely define the size. We consider the particle size to be identifiable if there exists a model {111} face that was matched and whose neighboring {111} faces were also matched. Otherwise, the identified shape is moved and resized while maintaining the orientation in order to identify the best fit for the point cloud. The direction for the origin movement is chosen as n

dM =

−∑i M,oct Ra i n

∑i M,oct Ra i

(27)

where Rai corresponds to the nM,oct matched crystal normals of the {111} face group. The sign of this vector may be inverted if 2692



Article

Figure 12. Large primary particle was split into the green and blue region (a), leading to highly overlapping large primary particles (c). The problematic regions are detected and concatenated (b) to provide only one primary particle (d).

lapping particles are merged, and this is followed by the same postprocessing steps described above. These postprocessing steps involve shape identification and deletion in the case of identification failure. As the merged particles may also require resizing, the procedure is iterated until no pairs of particles to be merged are found, as illustrated in Figure 10. This completes the shape identification of each primary particle in an aggregate. 5.3. Measuring Fit Quality. 5.3.1. Simulated Data. For the simulated data, the goodness of fit can be measured by comparing the simulated polytopes with the identified polytopes after adjusting the identified polytopes for the simulated μCT resolution. This enables an evaluation of the accuracy of shape identification. For each simulated primary particle, the algorithm finds a match among the identified primary particles by comparing the distance between particle centers. Fit quality is evaluated using two measures. The volume deviation is defined as

Figure 13. Shape identification before (a) and after (b) resizing the unidentifiable smaller particle.

determined that the origin should be moved toward the background. In each iteration, the crystal origin x0,i is moved according to x 0, i + 1 = x 0, i + dM

(28)

In the first iteration, x0,1 = x0. The search is stopped if the origin movement does not result in an improvement in comparison to the previous iteration according to the following measure. qp =

1 Nsurf

qv,s =

i=1

(30)

where Vi and Vs are the volumes of the identified and simulated polytope, respectively. Overlap deviation is measured by

Nsurf

∑ ∥pi − pi ,proj∥

|Vi − Vs| Vs

(29)

qo,s = 1 −

Here, pi are the Nsurf surface points of the primary particle and pi,proj are their projections to the polytope with the chosen best size. Resizing the shape with the currently considered origin x0,i is performed using a bisection algorithm. In the first step, a large crystal denoted by the distance vector hL and a small crystal denoted by hS are chosen as bisection boundaries. All surface points must be inside hL, and none may be inside hS. In each bisection iteration j, a new crystal defined by hj = 0. 5 ·(hL + hS) is evaluated. If more boundary points are outside than inside the crystal denoted by hj, this shape is set as the new small-sized boundary, hS. Otherwise, it is set as the new large-sized boundary hL. The procedure is stopped if the size difference between the large and small crystals is smaller than 0.1. In general, this procedure results in one of the shapes that correctly describe the primary particle orientation, but not its size, as discussed above. However, small parts of additional faces that were not identified by the Hough transform, but represent the sought boundary of the identified {111} faces, are often visible on the primary particle surface. In these cases, the quality measure defined in eq 29 will ensure that the best fit also covers these points, resulting in the correct description of the particle size as well. Moving and resizing is illustrated in Figure 13. The particle overlap issue described in the previous section may only become evident after resizing the particles. Therefore, over-

Vo,s,i Vs

(31)

where Vo,s,i is the volume of the overlap between the identified and simulated polytope. Combining the measures enables capturing the displacement as well as the deviation in size and shape. Shape identification is considered successful if qv,s < 0.15 and qo,s < 0.15 for each primary particle in an aggregate and if the number of primary particles was identified correctly. Furthermore, if at least one primary particle did not have a match among the identified particles, while the shape identification was successful for others, the aggregate is declared to contain particles for which shape fit was impossible. If a simulated primary particle overlaps to a degree of at least 80% with another primary particle, it is not considered since shape identification in this case is deemed impossible. 5.3.2. Real Data. As it is no longer possible to compare the obtained polytopes to the simulated polytopes, fit quality measures for the real data are based on comparison with the original point cloud. The mean quadratic distance deviation is defined as

qr,all =

1 Nsurf,agg

N

∑ j surf,agg min ∥pj − pj ,proj, i∥2 i

max hi ,oct i

2693

(32) DOI: 10.1021/acs.cgd.5b01806 Cryst. Growth Des. 2016, 16, 2685−2699


Article

Table 1. Shape-Identification Success Evaluated for the 12 Simulated Datasets Using the Quality Measures for Simulated Data, Defined in Section 5.3.1, and the Quality Measures for Experimental Data, Defined in Section 5.3.2a simulation parameters

quality measures - simulated data

quality measures - real data

growth time t

surface roughness proughness

[%] agg. all well identified

[%] agg. with badly identified

[%] agg. with unidentified

[%] agg. all well identified

[%] agg. with badly identified

[%] agg. with unidentified

50 50 50 50 150 150 150 150 300 300 300 300

0 0.2 0.4 0.6 0 0.2 0.4 0.6 0 0.2 0.4 0.6

95 95 97 96 94 95 92 94 82 85 85 83

4 4 3 4 6 4 7 5 14 11 11 12

1 1 0 0 0 1 1 1 4 4 4 5

95 95 97 96 96 97 97 97 95 94 93 91

5 4 3 4 4 2 3 2 5 4 7 7

0 1 0 0 0 1 0 1 0 2 0 2

a The first two columns show parameters used in sample generation. The third and sixth columns represent the percentage of aggregates where the shape of all primary particles was correctly identified. The fourth and seventh columns represent the percentage of aggregates containing at least one badly identified primary particle. The fifth and eight columns show the percentage of aggregates where the shape of some primary particles could not be identified.

Here, each of the Nsurf,agg aggregate surface points pj is projected to that identified polytope i which results in the smallest distance between pj and its projection pj,proj,i. The point can be projected to the polytope face, edge, or vertex, as explained in our previous work.10 The obtained error value is scaled by the {111} face distance of the largest identified primary particle, max hi,oct .

for different simulation parameters are shown in Figures 2 and 3. The percentage of cases in which the shape of at least one primary particle was identified badly or could not be identified is shown in columns four and five of Table 1. The main reasons for the observed issues are discussed below. Unsuccessful shape identification is often caused by faulty segmentation, which results from allowing the wrong pair of regions to be concatenated when considering a concavity point at the boundary of more than two regions. This is illustrated in Figure 14, where the concatenation of the blue regions leads to

i

A similar measure qr,fit is defined using only the Nsurf,fit surface points that belong to watershed regions for which shape identification was possible:

qr,fit =

1 Nsurf,fit

N

∑ j surf,fit min ∥pj − pj ,proj, i∥2 i

max hi ,oct i

(33)

Finally, a displacement measure for the real data is defined as qr,v =

(Vagg − Ninside) Noutside + Nvoxels Nvoxels

(34)

Here, Nvoxels is the total number of voxels in the aggregate, representing its volume; Ninside and Noutside are the number of voxels inside and outside all identified polytopes, respectively; and Vagg is the volume filled by the identified polytopes. This measure identifies displacement. Shape identification is considered successful if qr,all < 0.15 and qr,v < 0.2. Aggregates with particles where no shape fit was possible are identified as those where successful identification conditions are not fulfilled, but where qr,fit < 0.15 and qr,all > qr,fit.

Figure 14. Example of a badly identified primary particle (a) due to the faulty segmentation (b). The concavity points and their associated cubes are shown in red and black, respectively. The concavity point between the two blue and one green region wrongly allowed the concatenation of the blue regions. The simulation parameters were t = 50 and proughness = 0, the obtained error measures are qr,all = 0.186, qr,v = 0.310.

6. RESULTS AND DISCUSSION 6.1. Simulated Data. Goodness of fit was evaluated for the 12 obtained data sets, each of which consisted of 100 particles, using the comparison to the simulated polytope, as explained in Section 5.3.1, along with the measures developed for real data, as defined in Section 5.3.2. The results are summarized in Table 1. In between 82% and 97% of the cases, the shape of each primary particle was successfully identified, as seen in the third column of Table 1. Examples of successful shape identification

the badly identified primary particle. A similar issue may occur if a pair of primary particles lightly touch, as this can lead to an insufficient number of foreground voxels in the concavity computation mask which was introduced in Section 4.1. In the case of multiple primary particles, this issue cannot always be resolved by the steps described in Section 3, as both of the problematic particles may overlap to a high degree with another particle in the aggregate. 2694



Article

contribute little to the overall aggregate volume, as shown in Figure 15. This situation becomes more probable if particles overlap strongly. Furthermore, the measures are less sensitive than the ones defined in Section 5.3.1 because they are based on the aggregate as a whole and not on evaluating the fit quality for each primary particle. Therefore, inconsistencies created by a faulty shape identification of only one particle may have a small influence on the overall aggregate containing multiple particles. The number of aggregates containing particles for which shape identification was impossible was also underestimated by the quality measures developed for the real data. This can be seen by comparing the fifth and eight columns of Table 1. This can happen if primary particles were either not separated or the regions were wrongly concatenated, but the shape of one of the primary particles is identified correctly, as illustrated in Figure 16. If no region for which shape identification was impossible exists, the aggregate will be classified as having either a successful or unsuccessful shape identification, depending on how much the particle with missing fit contributes to the total particle volume. The percentage of successful shape identifications for aggregate classes consisting of two to five primary particles is illustrated in Table 2. The values are computed by considering

The middle part of Table 1 illustrates that simple boundary roughness does not have a large effect on the number of wellidentified aggregates. This is because the boundary roughness may influence the concavity value of unfavorably placed concavity points both positively and negatively, although the number of such cases is small. However, small objects on the surface, which may create problems for real crystals, were not simulated here. The degree of overlap between primary particles, resulting from growth, exhibited a much stronger influence on the number of successfully identified aggregates. This can be expected because both segmentation and shape identification become more difficult as the area of the particle surface which remains for shape identification becomes smaller, as seen in the example shown in Figure 15. Furthermore, the concavities between highly overlapping particles can in certain cases become less pronounced so that no concavity points are identified, as shown in Figure 16.

Table 2. Percentage of Successful Shape Identifications Per Aggregate Class Consisting of Two to Five Primary Particlesa number of primary particles

Figure 15. Shape identification (a) and segmentation into primary particles (b) for one aggregate. A part of the small yellow particle was concatenated with the orange region, leaving little information for its shape identification. The simulation parameters were t = 300 and proughness = 0.2. Shape identification was considered unsuccessful by comparison to simulated polytopes. However, it was misclassified as successful when using real-data-based quality measures (qr,all = 0.034 and qr,v = 0.077).

growth parameter

2

3

4

5

t = 50 t = 150 t = 300

100 100 98.71

96.15 94.23 88.46

82.14 82.14 60.71

100 78.13 28.13

a Values are computed for each growth parameter t by considering all four datasets obtained by varying the surface roughness parameter.

all four data sets obtained by varying proughness for some growth parameter t. It is seen that for all degrees of overlap, defined by the growth parameter t, successful shape identification becomes less probable for an aggregate with a large number of primary particles. This occurs because not all primary particles remain clearly visible or concavity points and their influence on region merging cannot be determined correctly, as illustrated in examples above. Surprisingly, for t = 50, all 8 aggregates consisting of 5 primary particles for all 4 different proughness levels were identified correctly. Note that the presence of {100} and {110} faces in the identified shapes in images of simulated octahedral particles is an artifact of searching for the plane with the highest point density in order to find the distance of unmatched existing faces, as explained in Sections 5.1.1 and 5.1.2. 6.2. Real Data. After separating lightly touching aggregates as described in Section 3.2.2, the real data contained 85 aggregates. The algorithm classified 65 shape identifications as successful, 9 as unsuccessful, and 11 as containing primary particles for which no shape identification was possible. Examples of successfully identified primary particle shapes can be found in Figure 17. While most primary particles in an aggregate are of similar size, aggregates containing particles with very different sizes, such as in Figure 17e, are also possible. The last row contains examples of particles whose shapes were judged to be successfully identified using the chosen thresholds,

Figure 16. Example of an aggregate where no separation between a pair of primary particles occurred (orange) due to a missing concavity point. However, the identified shape fits one of the two particles well. The simulation parameters were t = 300 and proughness = 0.6. Comparison to the simulated polytopes revealed a particle for which no shape identification was possible, whereas the shape identification was misclassified as successful when using real-data-based quality measures (qr,all = 0.090 and qr,v = 0.103). (a) Identification. (b) Segmentation.

A comparison between the middle and right parts of Table 1 shows that the successful classification of the identified shapes using the measures developed for real data in Section 5.3.2 also depends on the degree of particle overlap. The main reason for this behavior is that these measures are not capable of detecting problems if both the actual particle and its identified shape 2695



Article

Figure 17. Examples of successfully identified particle shapes from the real data. The values of qr,fit are not shown as they are equal to the qr,all values for all displayed particles. The identified shape of each primary particle is shown as a black outline. (a) qr,all = 0.063, qr,v = 0.120; (b) qr,all = 0.040, qr,v = 0.099; (c) qr,all = 0.034, qr,v = 0.091; (d) qr,all = 0.043, qr,v = 0.110; (e) qr,all = 0.067, qr,v = 0.146; (f) qr,all = 0.058, qr,v = 0.106; (g) qr,all = 0.033, qr,v = 0.074; (h) qr,all = 0.046, qr,v = 0.112; (i) qr,all = 0.087, qr,v = 0.198; (j) qr,all = 0.033, qr,v = 0.090; (k) qr,all = 0.053, qr,v = 0.114; (l) qr,all = 0.045, qr,v = 0.147.

aggregate for which it is difficult to determine whether the larger of the particles is in fact an aggregate or not. In such a case, a fit to a symmetric shape model will create a large error qr,v. Similarly, it is difficult to determine how many primary particles are present in Figure 19b, where the surface defects also cause concavity points that split the largest particle. The primary particles in Figure 19c are asymmetric, causing an error in qr,v. The particle in 19d was split into multiple parts due to surface roughness and attached small particles. Some of the obtained regions were not fittable, which resulted in them being excluded before measuring qr,fit. The attached small particles also caused qr,all and qr,fit to exceed 0.15. Finally, the last row of Figure 19 shows two particles for which the shape identification was bad due to faulty segmentation. 6.3. Time Requirements. The times required for sample preparation and imaging are in the order of a few hours. Both are strongly dependent on the chosen image resolution, sample size, and whether aggregates need to be prechosen for the measurement, or both single crystals and aggregates are imaged. Binarization, isolating aggregates, and splitting lightly touching aggregates are procedures that last up to a few minutes. Time requirements for the procedure consisting of preprocessing, segmentation into primary particles, and shape identification were measured on a desktop computer with 16

even though some issues were observed. In Figure 17j, a primary particle was divided into two parts; however, the overlap of the two shapes was not sufficient to concatenate the two parts back into one as described in Section 5.2.2. In Figure 17k, a small part of a third primary particle is visible but is too small to be identified. Figure 17l shows an aggregate of three primary particles in which the size of the smallest primary particle was overestimated. These observed issues are undetectable by the classification procedure because the problematic primary particle contributes little to the aggregate volume, as discussed in the previous section. Figure 18 illustrates 6 out of 11 aggregates that are considered to contain primary particles for which shape identification was impossible. Unlike the undetected small particle in Figure 17k, the particle shown in Figure 18e was segmented correctly, thus contributing to different values between qr,all and qr,fit, allowing the classification of a missing primary particle fit. The algorithm classified the aggregate in Figure 18f as containing a missing primary particle fit, even though this is not the case. This is due to surface imperfections, which caused an additional region to be considered as a primary particle, along with crystal asymmetry, which caused the value qv,r to be larger than 0.2. Figure 19 illustrates 6 out of 9 aggregates for which the shape identification was considered unsuccessful. Figure 19a shows an 2696



Article

Figure 18. Examples of aggregates for which the shapes of some primary particles were not identified. The identified shape of each primary particle is shown as a black outline. (a) qr,all = 0.166, qr,fit = 0.060, qr,v = 0.233; (b) qr,all = 0.163, qr,fit = 0.077, qr,v = 0.274; (c) qr,all = 0.205, qr,fit = 0.078, qr,v = 0.281; (d) qr,all = 0.474, qr,fit = 0.079, qr,v = 0.477; (e) qr,all = 0.076, qr,fit = 0.073, qr,v = 0.202; (f) qr,all = 0.108, qr,fit = 0.108, qr,v = 0.251.

Figure 19. Examples of badly identified primary particle shapes. Values qr,fit differ from those of qr,all in certain cases, meaning that there were regions in which no shape identification was possible. The identified shape of each primary particle is shown as a black outline. (a) qr,all = 0.188, qr,fit = 0.188, qr,v = 0.537; (b) qr,all = 0.087, qr,fit = 0.087, qr,v = 0.229; (c) qr,all = 0.105, qr,fit = 0.105, qr,v = 0.266; (d) qr,all = 0.181, qr,fit = 0.159, qr,v = 0.319; (e) qr,all = 0.053, qr,fit = 0.053, qr,v = 0.749; (f) qr,all = 0.043, qr,fit = 0.044, qr,v = 0.389.

GB of RAM and the Intel Core(TM)i5−3470 CPU @ 3.20 GHz, running Matlab 2015b on the 64-bit Windows 8.1 Enterprise. Computation was parallelized using the Matlab parallelization toolbox with four active workers, corresponding to four cores. Aggregates are processed by different workers as

there is no interdependence between different aggregate images. The average execution time for simulated data sets, each containing 100 aggregates, was about 35 min, whereas the execution time for the data set containing 85 imaged real particles was about 33 min. This leads to the sequential 2697



Article

(2) Simon, L. L.; et al. Assessment of Recent Process Analytical Technology (PAT) Trends: A Multiauthor Review. Org. Process Res. Dev. 2015, 19, 3−62. (3) Ma, C. Y.; Wang, X. Z.; Roberts, K. J. Morphological population balance for modeling crystal growth in face directions. AIChE J. 2008, 54, 209−222. (4) Singh, M. R.; Ramkrishna, D. A Comprehensive Approach to Predicting Crystal Morphology Distributions with Population Balances. Cryst. Growth Des. 2013, 13, 1397−1411. (5) Reinhold, A.; Briesen, H. High dimensional population balances for the growth of faceted crystals: Combining Monte Carlo integral estimates and the method of characteristics. Chem. Eng. Sci. 2015, 127, 220−229. (6) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. Model-based object recognition to measure crystal size and shape distributions from in situ video images. Chem. Eng. Sci. 2007, 62, 1430−1441. (7) Singh, M. R.; Chakraborty, J.; Nere, N. J.; 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. (8) Schorsch, S.; Ochsenbein, D. R.; Vetter, T.; Morari, M.; Mazzotti, M. High accuracy online measurement of multidimensional particle size distributions during crystallization. Chem. Eng. Sci. 2014, 105, 155−168. (9) Borchert, C.; Temmel, E.; Eisenschmidt, H.; Lorenz, H.; SeidelMorgenstern, A.; Sundmacher, K. Image-Based in Situ Identification of Face Specific Crystal Growth Rates from Crystal Populations. Cryst. Growth Des. 2014, 14, 952−971. (10) Kovačević, T.; Reinhold, A.; Briesen, H. Identifying Faceted Crystal Shape from Three-Dimensional Tomography Data. Cryst. Growth Des. 2014, 14, 1666−1675. (11) Schorsch, S.; Hours, J.-H.; Vetter, T.; Mazzotti, M.; Jones, C. N. An optimization-based approach to extract faceted crystal shapes from stereoscopic images. Comput. Chem. Eng. 2015, 75, 171−183. (12) Zhang, R.; Ma, C. Y.; Liu, J. J.; Wang, X. Z. On-line measurement of the real size and shape of crystals in stirred tank crystalliser using non-invasive stereo vision imaging. Chem. Eng. Sci. 2015, 137, 9−21. (13) Le Borne, S.; Eisenschmidt, H.; Sundmacher, K. Image-based analytical crystal shape computation exemplified for potassium dihydrogen phosphate (KDP). Chem. Eng. Sci. 2016, 139, 61−74. (14) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. An algorithm for analyzing noisy, in situ images of high-aspect-ratio crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61, 5236−5248. (15) Ahmad, O. S.; Debayle, J.; Gherras, N.; Presles, B.; Févotte, G.; Pinoli, J.-C. Quantification of overlapping polygonal-shaped particles based on a new segmentation method of in situ images during crystallization. J. Electron. Imaging 2012, 21, 021115−1−021115−11. (16) Borchert, C.; Sundmacher, K. Crystal Aggregation in a Flow Tube: Image-Based Observation. Chem. Eng. Technol. 2011, 34, 545− 556. (17) Faria, N.; Pons, M. N.; Feyo de Azevedo, S.; Rocha, F. A.; Vivier, H. Quantification of the morphology of sucrose crystals by image analysis. Powder Technol. 2003, 133, 54−67. (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) Ochsenbein, D. R.; Vetter, T.; Schorsch, S.; Morari, M.; Mazzotti, M. Agglomeration of Needle-like Crystals in Suspension: I. Measurements. Cryst. Growth Des. 2015, 15, 1923−1933. (20) Terdenge, L.-M.; Heisel, S.; Schembecker, G.; Wohlgemuth, K. Agglomeration degree distribution as quality criterion to evaluate crystalline products. Chem. Eng. Sci. 2015, 133, 157−169. (21) Collier, A. P.; Hetherington, C. J. D.; Hounslow, M. J. High resolution imaging of crystalline agglomerates. Chem. Eng. Res. Des. 1996, 74, 759−764.

execution time of between 1 and 2 min per aggregate. It depends strongly on the image resolution and aggregate complexity.

7. CONCLUSIONS AND OUTLOOK We presented an approach for identifying the size and shape of primary particles from μCT images of crystal aggregates. The method is based on merging watershed transform regions by considering concavity points, which serve as indicators of particle intersection. To evaluate its efficiency, the presented algorithm was tested on a sample containing 85 potash alum aggregates as well as on 12 data sets, each containing 100 simulated aggregates. The aggregates contained up to 5 primary particles. The results obtained using the simulated data show that the shapes of all primary particles were identified successfully for between 82% and 97% of the aggregates in the data set. The success rate decreased as the degree of overlap and the number of primary particles in the aggregate increased. An analysis of the real data shows that while shape identification was successful for most aggregates, the analyzed crystals do not always exhibit perfect symmetry, which may lead to a deviation from the identified symmetric shape. The examples presented in Figure 17 indicate that some aggregates contain primary particles that have the same orientation, such as in Figure 17g, whereas in others the primary particles seem to be oriented randomly toward each other, such as in Figure 17h. Therefore, the presented shape identification enables the analysis of the primary particle orientation and can lead to conclusions about the underlying aggregation mechanism. However, the method is limited to aggregates containing a small number of primary particles. Our future work will focus on investigating primary particle orientations in potash alum using improved sampling and preparation methods that ensure statistical relevance and no contact between different aggregates. In principle, the algorithm is applicable to other crystalline compounds, provided that sufficient symmetry is present to determine the shape of a primary particle from its visible part.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) through the Priority Programme SPP 1679, “Dynamic Simulation of Coupled Solids Processes”. The authors would like to thank Daniela Eixenberger, Sophie Bindl, and Florian Ernstberger for their assistance during the experimental and image processing part of this work, as well as Alexander Reinhold for fruitful discussions. F.P. and J.S. acknowledge financial support through the DFG Gottfried Wilhelm Leibniz program and the support of the TUM Institute for Advanced Study, funded by the German Excellence Initiative.

■

REFERENCES

(1) Dandekar, P.; Kuvadia, Z. B.; Doherty, M. F. Engineering Crystal Morphology. Annu. Rev. Mater. Res. 2013, 43, 359−386. 2698



Article

(43) MathWorks, MATLAB, and v. R2015b. (44) Reinhold, A. Crystal Shape Modeling and Convex Geometry Analysis of Geometric State Spaces. Ph.D. thesis, Technische Universität München, 2015. (45) Fukuda, K. cddlib, v 094g; 2012; https://www.inf.ethz.ch/ personal/fukudak/cdd_home/. (46) Gonzalez, R. C.; Woods, R. E.; Eddins, S. L. Digital Image Processing Using Matlab, 2nd ed.; McGraw Hill Education, Gatesmark, LLC: India, 2010; pp 542−550. (47) Hammer, P. Marching Cubes, 2011; http://www.mathworks. com/matlabcentral/fileexchange/32506-marching-cubes/content/ MarchingCubes.m. (48) FEI VSG, AVIZO Fire, v 8.1.0. (49) Fraunhofer ITWM, MAVI - Modular Algorithms for Volume Images, v 1.4.1. (50) Neubeck, A.; Van Gool, L. Efficient non-maximum suppression. 18th International Conference on Pattern Recognition, Vol 3; Hong Kong, China; August 20−24, 2006; pp 850−855. (51) Burger, W.; Burge, M. Principles of Digital Image Processing - Core Algorithms; Springer, Undergraduate Topics in Computer Science: London, 2009; p 59. (52) Hough, P. V. C. Methods and Means for Recognizing Complex Patterns; US Patent No 3,069,654, 1962. (53) Borrmann, D.; Elseberg, J.; Lingemann, K.; Nüchter, A. The 3D Hough Transform for plane detection in point clouds: A review and a new accumulator design. 3D Res. 2011, 2, 32:1−32:13. (54) Markley, F. L. Attitude Determination Using Vector Observations and the Singular Value Decomposition. J. Astronaut. Sci. 1988, 36, 245−258. (55) Wahba, G. A. Least Squares Estimate of Satellite Attitude. SIAM Rev. 1965, 7, 409−409.

(22) Collier, A. P.; Hetherington, C. J. D.; Hounslow, M. J. Alignment mechanisms between particles in crystalline aggregates. J. Cryst. Growth 2000, 208, 513−519. (23) Mumtaz, H. S.; Hounslow, M. J.; Seaton, N. A.; Paterson, W. R. Orthokinetic Aggregation During Precipitation: A Computational Model for Calcium Oxalate Monohydrate. Chem. Eng. Res. Des. 1997, 75, 152−159. (24) Hounslow, M. J.; Mumtaz, H. S.; Collier, A. P.; Barrick, J. P.; Bramley, A. S. A micro-mechanical model for the rate of aggregation during precipitation from solution. Chem. Eng. Sci. 2001, 56, 2543− 2552. (25) Briesen, H. Aggregate Structure Evolution for Size-Dependent Aggregation by Means of Monte Carlo Simulations. KONA 2007, 25, 180−189. (26) Hounslow, M. J.; Wynn, E. J. W.; Kubo, M.; Pitt, K. Aggregation of growing crystals in suspension: I. Mumtaz revisited. Chem. Eng. Sci. 2013, 101, 731−743. (27) Ochsenbein, D. R.; Vetter, T.; Morari, M.; Mazzotti, M. Agglomeration of Needle-like Crystals in Suspension. II. Modeling. Cryst. Growth Des. 2015, 15, 4296−4310. (28) Briesen, H. Hierarchical characterization of aggregates for Monte Carlo simulations. AIChE J. 2006, 52, 2436−2446. (29) Spettl, A.; Bachstein, S.; Dosta, M.; Goslinska, M.; Heinrich, S.; Schmidt, V. Bonded-particle extraction and stochastic modeling of internal agglomerate structures, Adv. Powder Technol. 2016 Submitted for publication. (30) Irshad, H.; Veillard, A.; Roux, L.; Racoceanu, D. Methods for Nuclei Detection, Segmentation, and Classification in Digital Histopathology: A Review - Current Status and Future. IEEE Rev. Biomed. Eng. 2014, 7, 97−114. (31) Vincent, L.; Soille, P. Watersheds in digital spaces - an efficient algorithm based on immersion simulations. IEEE Trans. Pattern Anal. Mach. Intell. 1991, 13, 583−598. (32) Adiga, P. S. U.; Chaudhuri, B. B. An efficient method based on watershed and rule-based merging for segmentation of 3-D histopathological images. Pattern Recognit. 2001, 34, 1449−1458. (33) Wählby, C.; Sintorn, I.-M.; Erlandsson, F.; Borgefors, G.; Bengtsson, E. Combining intensity, edge and shape information for 2D and 3D segmentation of cell nuclei in tissue sections. J. Microsc. 2004, 215, 67−76. (34) Long, F.; Peng, H.; Myers, E. Automatic segmentation of nuclei in 3D microscopy images of C. elegans; 2007 4th IEEE International Symposium on Biomedical Imaging: Macro to Nano, Vols 1−3; Arlington, VA; April 12−15, 2007; pp 536−539. (35) Kumar, S.; Ong, S. H.; Ranganath, S.; Ong, T. C.; Chew, F. T. A rule-based approach for robust clump splitting. Pattern Recognit. 2006, 39, 1088−1098. (36) Wienert, S.; Heim, D.; Saeger, K.; Stenzinger, A.; Beil, M.; Hufnagl, P.; Dietel, M.; Denkert, C.; Klauschen, F. Detection and Segmentation of Cell Nuclei in Virtual Microscopy Images: A Minimum-Model Approach. Sci. Rep. 2012, 2, 1−7. (37) Farhan, M.; Yli-Harja, O.; Niemistö, A. A novel method for splitting clumps of convex objects incorporating image intensity and using rectangular window-based concavity point-pair search. Pattern Recognit. 2013, 46, 741−751. (38) Indhumathi, C.; Cai, Y. Y.; Guan, Y. Q.; Opas, M. An automatic segmentation algorithm for 3D cell cluster splitting using volumetric confocal images. J. Microsc. 2011, 243, 60−76. (39) Qi, J. Dense nuclei segmentation based on graph cut and convexity-concavity analysis. J. Microsc. 2014, 253, 42−53. (40) Fernàndez, G.; Kunt, M.; Zrÿd, J.-P. A new plant cell image segmentation algorithm. Image Analysis and Processing. 8th International Conference, ICIAP ’95; September 13−15, 1995; San Remo, Italy; pp 229−34. (41) Zhang, C.; Sun, C.; Pham, T. D. Segmentation of clustered nuclei based on concave curve expansion. J. Microsc. 2013, 251, 57−67. (42) Reinhold, A.; Briesen, H. Convex Geometry for the Morphological Modeling and Characterization of Crystal Shapes. Part. Part. Syst. Charact. 2011, 28, 37−56. 2699


Shape Identification of Primary Particles in Potash Alum Aggregates

Recommend Documents