Langmuir 2006, 22, 2005-2015
2005
Quantitative Imaging of Aggregated Emulsions Robert Penfold,* Andrew D. Watson, Alan R. Mackie, and David J. Hibberd Institute of Food Research, Norwich Research Park, Colney, Norwich NR4 7UA, United Kingdom ReceiVed October 7, 2005. In Final Form: December 20, 2005 Noise reduction, restoration, and segmentation methods are developed for the quantitative structural analysis in three dimensions of aggregated oil-in-water emulsion systems imaged by fluorescence confocal laser scanning microscopy. Mindful of typical industrial formulations, the methods are demonstrated for concentrated (30% volume fraction) and polydisperse emulsions. Following a regularized deconvolution step using an analytic optical transfer function and appropriate binary thresholding, novel application of the Euclidean distance map provides effective discrimination of closely clustered emulsion droplets with size variation over at least 1 order of magnitude. The a priori assumption of spherical nonintersecting objects provides crucial information to combat the ill-posed inverse problem presented by locating individual particles. Position coordinates and size estimates are recovered with sufficient precision to permit quantitative study of static geometrical features. In particular, aggregate morphology is characterized by a novel void distribution measure based on the generalized Apollonius problem. This is also compared with conventional Voronoi/Delauney analysis.
1. Introduction Since the modern reintroduction of direct observation in 1973,1 optical videomicroscopy coupled with digital image processing methods have frequently been used to study the particle interactions, structure, and dynamics of colloidal suspensions.2,3 Interest has primarily centered on micron sized monodisperse systems as physical models for atomic and molecular manybody processes in condensed matter.4 In particular, the phase diagram of liquidlike and solidlike states with coexistence regions and the associated disorder-order structural transitions have received attention from several groups.5-7 Nonequilibrium behavior has also been the subject of recent research.8-10 In the dilute regime ( 0.55), (14) Marcus, A. H.; Rice, S. A. Phys. ReV. E 1997, 55, 637. (15) Marcus, A. H.; Schofield, J.; Rice, S. A. Phys. ReV. E 1999, 60, 5725. (16) Crocker, J. C.; Grier, D. G. J. Colloid Interface Sci. 1996, 179, 298. (17) Lu¨thi, Y.; Ricˇka, J. J. Colloid Interface Sci. 1998, 206, 302. (18) Bongers, J.; Manteufel, H.; Versmold, H.; Vondermassen, K. J. Chem. Phys. 1998, 108, 9937. (19) Bongers, J.; Versmold, H. J. Chem. Phys. 1996, 104, 1519. (20) Kasper, A.; Bartsch, E.; Sillescu, H. Langmuir 1998, 14, 5004. (21) Valentine, M. T.; Kaplan, P. D.; Thota, D.; Crocker, J. C.; Gisler, T.; Prud’homme, R. K.; Beck, M.; Weitz, D. A. Phys. ReV. E 2001, 64, 061506.
10.1021/la052719w CCC: $33.50 © 2006 American Chemical Society Published on Web 01/27/2006
2006 Langmuir, Vol. 22, No. 5, 2006
fluorescence CLSM and digital image analysis have been employed to quantify the local structure of colloidal glasses quenched from refractive index matched silica dispersions.23 Dynamical heterogeneities of dense and glassy states have also been studied by these methods24,25 as well as the kinetics of crystallization processes.26 An interesting attempt to apply confocal microscopy and image analysis methods for structurally characterizing a dense paste of cornflour starch with highly irregular and polydisperse particles has also appeared recently,27 while the quantitative study of pore and void morphology in acid-set milk protein gels as a function of sucrose content has also been reported.28 Moving beyond particulate systems with purely repulsive interactions, Dinsmore and co-workers29,30 have applied the combination of CLSM, the Crocker-Grier digital image processing algorithms and multiple particle tracking, all in 3D, to moderately concentrated (φ < 0.1) suspensions flocculated by short ranged attractive forces effectively induced through the depletion mechanism. Varadan and Solomon have similarly examined the high-density regime (0.26 < φ < 0.50), and carried out a Voronoi volume analysis to establish long-range structural heterogeneity,31 as well as the formation of voids and fractures under unconstrained uniaxial compression (or squeeze flow).9 The present work is primarily concerned with the quantitative estimation of center coordinates and sizes for large numbers of aggregated polydisperse spheres by optical tomography. In principle, this is a direct generalization of the previous methods,16 but polydispersity presents a significant additional challenge since a single well-defined length scale, which affords a simple mask for object location in the Crocker and Grier algorithm, is no longer available. A fixed mask based on a large scale will not differentiate a cluster of small features. Conversely, a small mask can amplify noise by identifying random brightness variations with a cluster of artifacts in place of a single large object. Brightness smoothing to reduce variations across large features in turn degrades small details elsewhere in a polydisperse image. Even where a local brightness maximum is successfully associated with a single object, the size relationship of the feature and the fixed mask is unclear. Furthermore, aggregated systems exacerbate the segmentation problem where poorly resolved features in close proximity need to be discriminated. It is especially difficult to resolve an aggregate comprising several small features neighboring a single large object. Brujic´ et al.32 have tackled these problems in a CLSM study of interparticle forces in a concentrated polydisperse emulsion compressed to a “jammed” state where all droplets are in contact with their neighbors. By considering the emulsion as a collection of spheres convolved with a set of delta functions, the Fourier filtering method was applied to extract droplet positions and sizes with sub-voxel accuracy. An alternative algorithm is presented here based on the Euclidean distance map (EDM)33 to locate closely situated polydisperse features in 3D. Given a binary image that segregates (22) Chestnut, M. H. Curr. Opin. Colloid Interface Sci. 1997, 2, 158. (23) van Blaaderen, A.; Wiltzius, P. Science 1995, 270, 1177. (24) Kegel, W. K.; van Blaaderen, A. Science 2000, 287, 290. (25) Weeks, E. R.; Crocker, J. C.; Levitt, A. C.; Schofield, A.; Weitz, D. A. Science 2000, 287, 627. (26) Gasser, U.; Weeks, E. R.; Schofield, A.; Pusey, P. N.; Weitz, D. A. Science 2001, 292, 258. (27) Bromley, E. H. C.; Hopkinson, I. J. Colloid Interface Sci. 2002, 245, 75. (28) Pugnaloni, L. A.; Matia-Merino, L.; Dickinson, E. Colloid Surf. B 2005, 42, 211. (29) Dinsmore, A. D.; Weeks, E. R.; Prasad, V.; Levitt, A. C.; Weitz, D. A. Appl. Opt. 2001, 40, 4152. (30) Dinsmore, A. D.; Weitz, D. A. J. Phys.: Condens. Matter 2002, 14, 7581. (31) Varadan, P.; Solomon, M. J. Langmuir 2003, 19, 509. (32) Brujic´, J.; Edwards, S. F.; Grinev, D. V.; Hopkinson, I.; Brujic´, D.; Makse, H. A. Faraday Discuss. 2003, 123, 207.
Penfold et al.
objects of interest from background, the corresponding EDM associates pixel intensity with the distance to the nearest background pixel. As the name suggests, the metric is just the usual norm on a Euclidean space. The distance transformation required can be efficiently computed for 2D or 3D images. In 2D, for example, applied to a bright disk on a dark background, the center pixel of the corresponding EDM has a value proportional to the disk radius (within the limits of discretisation). The next neighbor pixels in all directions are assigned values approximately one unit less than the central maximum and so on, until at the feature edge pixels are assigned values of one. Local EDM maxima and their positions provide, respectively, a measure of spatial extent and centroid coordinate estimates for image features. Since construction of the EDM requires a binary image, an intensity threshold must first be applied to the raw micrographs. Reliable results are crucially dependent on a correct choice of threshold parameter. If the threshold value is too high, many features will be missed, the sizes of those that are located will be underestimated, and single large features will be interpreted as an agglomerate of smaller objects. If the threshold is too low, discrimination between features will be poor: a cluster of closepacked small features will appear as a single large feature, ultimately as a consequence of limitations in image resolution. The threshold value must be tuned to maximize the utility of the EDM. Efficiency is gained by exploiting the a priori knowledge of emulsion systems that image features are spherical and nonoverlapping, though point contacts are allowed. A single dark pixel within an otherwise bright region is enough to differentiate a cluster of several small objects from a single large feature. Moreover, the EDM is better able to resolve features in 3D rather than in two-dimensional slices of three-dimensional specimens. 2. Experimental 2.1. Emulsion System. The emulsion formulation is described in detail elsewhere.34 Briefly, an oil-in-water emulsion has been developed with a disperse phase mixture of n-hexane and tetradecafluorohexane. Since the perfluorinated oil has a refractive index less than that of water, an emulsion with a good match in both buoyancy and refractive index can be achieved by appropriate adjustment of the oil composition. A diblock copolymer F6H6 was also added to lower the consolute temperature and suppress liquidliquid-phase separation of the mixture. The emulsion was stabilized using the nonionic surfactant Brij-35, contained sodium azide as a preservative, and was flocculated using poly(ethylene oxide) (PEO) of molecular weight Mw ≈ 900 kDa at concentrations up to the overlap limit C/p ≈ 0.9% w/w in the continuous phase. The disperse phase was stained with the fluorescent dye Nile Red. Homogenization by the low shear cross-flow membrane (XME) process35 was crucial for obtaining a relatively tight size distribution (mean particle diameter of 6.0 µm with a standard deviation of 2.5 µm) and for preventing vaporization of the volatile oil phase. 2.2. Image Acquisition. Three-dimensional images of the emulsion system were obtained using CLSM and optical tomography. A Bio-Rad 1024 confocal system mounted on a purpose built horizontal z-axis microscope using Nikon optics gave improved resolution in the direction of gravity (x-axis in a right-handed coordinate system). Micrographs were collected using an × 60 oil immersion objective with a numerical aperture of 1.4. An image stack comprised 512 × 512 × 64 voxels, each a rectangular prism with dimensions 0.21 × 0.21 × 0.48 µm, giving a total viewing field of 107.3 × 107.3 µm in the focal x-y plane and 30.7 µm along the (33) Russ, J. C. The Image Processing Handbook, 4th ed.; CRC Press LLC: Boca Raton, FL, 2004. (34) Hibberd, D.; Watson, A.; Mackie, A. to be published. (35) Joscelyne, S. M.; Tra¨gårdh, G. J. Membrane Sci. 2000, 169, 107.
QuantitatiVe Imaging of Aggregated Emulsions
Langmuir, Vol. 22, No. 5, 2006 2007
optical z-axis. Acquisition time for a complete scan was approximately 250 s. The sample was contained in a modified glass cuvette with dimensions 10 × 4 × 37 mm where one x-y face was ground off, replaced by a microscope cover slip and sealed with epoxy cement. The cuvette was filled to a height of 35 mm. Close refractive index matching permitted visualization of the sample interior rather than the surface layer. To limit edge effects, samples were probed at least 10 µm from the container walls and well below the emulsion/air interface.
3. Methods 3.1. Image Analysis. Following Crocker and Grier,16 we identify four stages of processing: noise discrimination, image restoration, segmentation (object identification) and refinement of feature locations. With fluorescence microscopy, correction for nonuniform background illumination is largely unnecessary, but random pixel noise remains, associated with detector response and digitization. Boxcar averaging is not appropriate for polydisperse and clustered objects where a single masking length scale is not well-defined. Indeed, without cascading, the box filter optical transfer function (OTF) exhibits negative values leading to contrast inversion, whereas the intrinsic anisotropy is exacerbated in 3D. Although specific for CLSM, we have reduced the effects of noise by employing a theoretical OTF to account for inherent physical limitations of the image formation system. Regularized deconvolution procedures enable 3D filtering of the optical slices to improve resolution along the specular axis. After rudimentary restoration operations, a novel application of the EDM combined with suitable thresholding steps provides effective segmentation of clustered objects. Feature location estimates are refined by appealing to the local intensity distributions. 3.1.1. Noise Discrimination and Reduction. To model the optical system, the microscope objective lens is assumed to be an apochromatic, linear and shift-invariant (isoplanatic) element. Therefore, the measured image radiance distribution g(r) can be expressed as a convolution of the actual fluorochrome distribution f(0)(r) and a total point spread function (PSF) K(r) thus
g(r) )
∫dr′K(r - r′)f(0)(r′) + w(r)
(1)
but where g is also polluted by random noise w(r) associated with the detection and/or recording processes, as well as the inherently stochastic (Poisson) nature of fluorescence emission. The total PSF describes the system response to an impulse (Dirac) distribution wherever it is located and accounts for image blurring or smearing as a result of diffraction limited optical resolution. For epi-illumination fluorescence imaging, where the short but variable fluorescence decay time produces spatially incoherent radiation, K(r) can be decomposed into the uncorrelated product of partial PSFs both characterized by the same lens numerical aperture (N.A. ) n sin R where n is the medium refractive index) but differing in the monochromatic light wavelength λ of excitation and emission
K(r) ) KR(r,λex)KR(r,λem)
(2)
An explicit expression for the partial OTF K ˆ R(k,λ) of a circular aperture A(R) in 3D can be obtained by analytically evaluating the autocorrelation integral of the Debye image amplitude36
K ˆ R(k,λ) )
∫ dk′ ∫ dr ∫A(R) dq ∫ dr′ ∫A(R) dq′
( (
exp -i r‚k + (r + r′)‚k′ + (36) Philip, J. J. Mod. Opt. 1999, 46, 1031.
2π (r‚q + r′‚q′) λ
)) (3)
where ψ ˆ (k) denotes the Fourier transform of ψ(r). By avoiding small angle approximations, this result is crucial for calculating the theoretical PSF of oil-immersion objectives with very high numerical aperture that are typically required for high-resolution fluorescence imaging. The straightforward Fourier inversion of (1) fails, however, since the bandlimited OTF K ˆ (k) has only finite bounded support that precludes unique solutions (i.e., arbitrary “invisible” objects exist which cannot be imaged) and also denies the existence of an inverse Fourier transform for ˆf(0)(k) in the presence of randomly fluctuating noise. To combat this ill-posed problem (in the Hadamard sense), the notion of an exact solution must be rejected in favor of an approximate f(r) from a set which is constrained by appealing to additional a priori knowledge. We have implemented such a Tikhonov regularisation scheme by introducing a linear filter based on the minimization of a cost functional which constrains the discrepancy, between approximate solution and data, with a smoothness condition on the norm of the solution and its first derivative, to obtain37
K ˆ *(k) gˆ (k) ˆf µ(k) ) 2 |K ˆ (k)| + µ(K ˆ max2 + ν k‚k)
(4)
where K ˆ max2 is the maximum square modulus of the total OTF. By interactively choosing the regularisation parameters µ > 0 and ν g 0, very substantial improvements have been achieved in the quality and resolution of the CLSM images as demonstrated on a test case protein gel shown in Figure 1. Panel (a) displays an image of a 4%w/w caseinate/water gel acidified by hydrolysis of 1 of g/g of glucon-δ-lactone (GDL) on a protein basis and dyed with rhodamine-B. Fluorescence was excited by the 514 nm line of an argon ion laser and emission measured at 598 nm. The raw data consisted of 1024 × 1024 × 16 voxels corresponding to a field of 48.8 × 48.8 × 4.8 µm collected through an oilimmersion lens with N.A. ) 1.4. The result of unconstrained deconvolution (µ ) ν ) 0) is indicated by Figure 1b where amplification of the high spatial frequency noise contribution has entirely corrupted the signal. This is confirmed by the corresponding 2D Fourier transforms shown in panels d and e, respectively, where “invisible” objects beyond the resolution of the OTF have been unphysically exaggerated. The dramatic improvement conferred by even the most simply constrained deconvolution (µ ) 10 and ν ) 0) is clearly demonstrated in panel (c) and verified in reciprocal space by the recovery of important short wavelength information as shown in panel f. 3.1.2. Image Restoration. Small numbers of extremely bright pixels biased the contrast of typical deconvolved “raw” images (RI), with potentially damaging consequences for the thresholding step. For each optical slice, therefore, contrast was enhanced by truncating the high-intensity tail and rescaling the histogram to the full brightness range. For successful application of the Euclidean distance transform, true spherical features must be represented as spheres in the corresponding image. That is, the diameter along any of the three coordinate directions must have the same measure in image units. To minimize the acquisition time of raw images (important particularly for moving objects), voxels are generally not regular (37) Bertero, M.; Boccacci, P. Introduction to InVerse Problems in Imaging; Institute of Physics Publishing: Bristol, U.K., 1998.
2008 Langmuir, Vol. 22, No. 5, 2006
Penfold et al.
Figure 1. CLSM image stack of an acidified dairy protein gel. Panel a shows the central slice raw image, whereas panels b and c illustrate the results of 3D deconvolution by a theoretical optical transfer function without and with regularisation, respectively. For visual clarity, panels d-f show the corresponding Fourier transform in 2D.
cubes, but are somewhat larger along the optical z axis. Moreover, the imaging process is not strictly confocal with a consequent lack of resolution in the specular direction so that the image of a sphere actually appears oblate. Isotropy is achieved by applying a linear scaling to interpolate voxels in the z direction. The scale factor, which depends on the microscope configuration, was determined manually by selecting a value yielding regular disks for scaled images of spheres viewed either in the x-z or y-z planes. This strategy assumes that all image features have the same aspect ratio of major to minor axes. Images were smoothed by convolution with an everywherepositive Gaussian ball in 3D, since multiple one-dimensional smoothing can introduce errors for a close-packed system. Finally, a cut was applied to remove background noise and arrive at a “processed image” (PI). A flowchart detailing the restoration algorithm is presented in Figure 2. 3.1.3. Segmentation. The segmentation algorithm proceeds as illustrated in Figure 2. First, a minimum feature size Rmin is specified and the number of size bins, indexed by n, is fixed. The maximum feature size can be set arbitrarily, but computational efficiency is optimized by a single pass through the EDM loop to establish Rmax. By applying a given brightness threshold tm to the PI, a binary image (BI) is formed, and the EDM subsequently obtained under the distance transform d. The EDM itself was used to first identify and then remove all features in the radius range Ri < R e Rf in decreasing order. This procedure used the EDM property that feature radius estimates appear as local maxima R0(r) at position r. Each feature at r is eliminated from the PI by blanking with a mask of size not less than Rf, before decrementing the current size bounds (Ri and Rf) and returning for the next iteration. In turn, this procedure is nested inside a brightness threshold loop, that was effective for balancing feature discrimination and overall capture with reasonable size estimates, to finally produce a quantified image (QI). For the present applications, three
successively smaller threshold values were chosen by inspecting line brightness profiles of the PI. The impact of a candidate threshold value was initially verified against the RI to ensure that at least minimal information of interstitial structure is retained in the BI. Since features gave a brightness distribution that was Gaussian rather than top-hat in cross section, bright spheres resulting from applying a threshold cut were automatically of smaller diameter than the droplets themselves. The final radii estimates likewise underestimate the true feature size. It was necessary that no “residue” remain after object removal, since this might be discovered on later iterations as a set of smaller features. To ensure complete removal, the size of the blanked region was artificially increased by an adjustable additive parameter. The output radius values were those resulting from the EDM analysis without the additive term. 3.1.4. Coordinate Refinement. Relative to conventional masking, use of the EDM algorithm gave a good estimate of the location of features sufficiently close to the image boundary that part of the feature was lost. In fact, if the center point of the feature was actually within the image volume the location estimate was as robust as that for any feature which appeared fully within the image. This was implemented by working with a brightnessinverted binary image, and computing the EDM of background instead of foreground features so that the exterior to the image volume was equivalent to a feature of infinite extent. Location estimates for features away from the boundaries of the image volume were refined by applying brightness masks of the appropriate size to the PI and calculating the brightness “center of mass”, the first moment of the brightness distribution. The second moment was not used to estimate sphere radii since this would have resulted in an inconsistent set of values between those features wholly visible in the image, for which brightness moment calculations were robust, and those partially outside the image volume where they were not. Instead, we have used feature radius estimates from the EDM itself. For each feature, the
QuantitatiVe Imaging of Aggregated Emulsions
(
Langmuir, Vol. 22, No. 5, 2006 2009
) ( ) ( )
(x1 - x2) (y1 - y2) (z1 - z2) M ) (x1 - x3) (y1 - y3) (z1 - z3) , (x1 - x4) (y1 - y4) (z1 - z4) ∆1 - ∆2 s ) ∆1 - ∆3 , ∆ 1 - ∆4
R1 - R 2 d ) R 1 - R3 R 1 - R4
(8)
with ∆i ) (xi2 + yi2 + zi2 - Ri2)/2. The matrix M is singular if and only if the Bi are coplanar whence B either does not exist or is not unique. Suppose det M * 0 and form M-1 to solve (7)
r ) r0 - Ru
(9)
where r0 ) M-1s is the center position of B in the special case that the Bi are monodisperse (d ) 0). More generally, u ) M-1d * 0 is a dimensionless direction that gives the required center as a displacement from r0 in units of R. It now remains only to solve the single quadratic equation in R obtained by substituting (9) into the first of (6)
(1 - u‚u)R2 + 2(R1 - u‚v)R + (R12 - v‚v) ) 0 (10)
Figure 2. Starting from a deconvolved raw image (RI), this flow diagram details algorithms for restoration (to arrive at a processed image (PI)), and subsequent segmentation by looping over feature lists (sorted by size R) generated from iterative application of thresholds tm (to give binary images (BI)) and the metric transformation d (to give Euclidean distance maps (EDM)). A quantified image (QI) is finally output.
intensity contribution was calculated as the zeroth moment of the brightness distribution using a suitable mask. 3.2. Morphological Characterization. 3.2.1. Void Distribution by Tangent Sphere Construction. We consider a 3D generalization for the problem of Apollonius38 that originally concerned the tangency of objects in the plane. Given four spherical balls Bi ∈R3 with radius Ri and centered at ri ) (xi, yi, zi) for i ) 1, 2, 3, 4, we require the “osculating spherical ball” B that is tangent to each Bi and remains external to all the Bi satisfying 4 B∩(∪i)1 Bi) ) Ø
(5)
This corresponds to exactly one of the 24 (possibly degenerate) solutions of the generalized Apollonius problem. Let B be located at r ) (x, y, z) with radius R so the system of constraints becomes
(x - xi)2 + (y - yi)2 + (z - zi)2 - (R + Ri)2 ) 0, (i ) 1, 2, 3, 4) (6) Eliminating quadratic terms in the unknowns by subtracting the first equation (i ) 1) obtains the underdetermined 3 × 3 linear system
Mr ) s - Rd where (38) Gisch, D.; Ribando, J. M. Am. J. Undergrad. Res. 2004, 3, 15.
where v ) r1 - r0. A straightforward analysis shows that the degenerate case of a mutual tangent plane yields the condition u‚u ) 1. After discounting this exception, an explicit solution is obtained as the positive root of (10). The position r is finally obtained from (9). From the configurational data liberated by the image analysis, the set {B} of osculating balls is constructed for every possible selection of four particles. Periodic boundary conditions are applied on the image volume and members of {B} that intersect any particle or that have a diameter exceeding the field of view are pruned away to obtain the subset {O} ⊂ {B} of overlapping holes (O-holes). The union of O-holes characterizes the pore space morphology of the configuration. A clearer representation of this void structure emerges by forming a further subset {I} ⊂ {O} of independent holes (I-holes) where O-holes are first ranked according to size, then those that intersect larger O-holes are systematically removed. Thus, I-holes provide a coarse indication of local pore size and position. Standard cluster analysis methods39 can also be applied to {O} for examining the larger scale distribution of voids. 3.2.2. Voronoi Analysis. Given a configuration of sites, the Voronoi tessellation of space associates a unique convex polyhedron to each such that any interior point is closer (in the Euclidean metric) to the site than to any other. In 3D, the Voronoi construction has been effectively applied to quantify the heterogeneous structure of a broad range of disordered materials40 including colloidal gels.9,31 Here, the Voronoi diagram is obtained from the convex hull generator (qhull) of Barber et al.,41 specifically designed for computational geometry with floating point arithmetic. Metric properties are subsequently calculated using the Gram determinant formula for the measure of a simplex with known vertexes.
4. Results and Discussion 4.1. EDM and Skeletons. Figure 3a shows a single complete 2D slice from a PI image stack, focused at a depth of 12 µm into the image volume, corresponding to 22 µm into the physical
(7) (39) Rapaport, C. C. The Art of Molecular Dynamics Simulation; Cambridge University Press: Cambridge, 1995. (40) Montoro, J. C. G.; Abascal, J. L. F. J. Phys. Chem. 1993, 97, 4211. (41) Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. T. ACM Trans. Math. Software 1996, 22, 469.
2010 Langmuir, Vol. 22, No. 5, 2006
Penfold et al.
Figure 3. Two-dimensional slice of a 3D CLSM image of a 30% oil in water emulsion flocculated with 0.9% w/w PEO is displayed in panel a. The lower left-hand quadrant indicated is expanded in panel b. Panel c shows the 2D Euclidean distance map corresponding to the image segment in panel b, whereas panel d shows the result of a morphological thinning operation (the white lines), applied in 2D to the EDM and superimposed onto the original image.
sample. The network of aggregated, polydisperse droplets is clearly visible. The lower left-hand quadrant is reproduced at twice the scale in Figure 3b. To illustrate the present segmentation method, a 2D EDM of this slice is shown in Figure 3c, where the region depicted corresponds once again to the lower lefthand quadrant. Bright centers correspond to large features of the PI. We emphasize that this 2D example was not used for 3D feature location, but serves only to visually demonstrate the method principles more clearly. Intensity maxima in the EDM (Figure 3c) are apparently linked by suggestive “brightness bridges”. These are ultimately a consequence of finite pixel size and microscope resolution, so that two droplets in direct contact appear to be joined by a fluidlike neck. Applying a conventional 2D morphological thinning operation to the EDM generates the skeleton depicted in Figure 3d, where it is superimposed on the original image. Similar skeletons can also be generated in 3D. 4.2. Coordinate Extraction. Although the present algorithm was computationally intensive for large numbers of small features, the segmentation stage (where an EDM was used to iteratively remove all features in a given size range) was found to be relatively efficient for small numbers of large features. Moreover, this
procedure addressed the inevitable problem of multiple “hits” where the EDM contained equal brightness points in close proximity. An alternative strategy, more efficient for small features, involved the well-known technique of comparing an image with a grayscale morphological dilation of itself. Conversely, this technique was computationally expensive for large features. A combination of the two processes introduced consistency problems. The EDM feature-finding algorithm was successfully used to quantify several 3D image stacks. For the example presented here, the technique located 3667 droplets ranging in diameter from σ ) 1.7 to 7.7 µm, where these “raw” EDM values systematically underestimate the true size. This is consistent with a mean diameter σ j ) 4.0 µm on the number distribution obtained by electric zone sensing (Coulter counter), and a nominal volume fraction of 30%. No attempt has been made to eliminate droplets on account of low brightness, though a reasonable minimum size cutoff is arbitrarily applied. Small droplet identification may not be robust, partly because they are difficult to distinguish from the residue left after larger features were removed. In principle, the EDM approach works for a single
QuantitatiVe Imaging of Aggregated Emulsions
Langmuir, Vol. 22, No. 5, 2006 2011
Figure 4. Panel a shows the same emulsion image quadrant as in Figure 3. Panel b shows the corresponding quantified image (QI), based on features located in the processed image (PI) by the feature-finder. The regions marked A and B indicate features that appear in the PI but which are apparently absent in the QI. The gray arrow indicates the location of the slice in Figure 5. Panel c shows the PI with all located features removed. The brightness has been rescaled to full range.
Figure 5. Image volume viewed perpendicular to the optical axis in the x-z (constant y) plane. Plane a shows the quantified image (QI), whereas panel b is the corresponding real processed image (PI). Note that the space axes are isotropic. The gray arrow indicates the large feature discussed in the text. The horizontal gray line marks the z value corresponding to the x-y slices of Figure 4. The absence of the large feature in region A of the QI which is nevertheless present in the PI indicates that the gray line intersects the true image of the large feature but not the measured counterpart. Note the overall smearing along the optical z axis in the PI.
isolated bright pixel, which in the present case would correspond to a feature of diameter σ ) 0.21 µm, though clearly unreachable in practice. Figure 4 shows the result of droplet location. Panel a is the same image quadrant depicted in Figure 3b. Figure 4b shows the corresponding QI, in which balls of appropriate radius were drawn at the locations output by the EDM feature-finding algorithm. The radii were scaled up by a factor of 1.4 relative to the “raw” EDM values to compensate for the inherent underestimation. A casual comparison of panels (a) and (b) shows that many features were successfully identified and located. The paucity of overlapping features in Figure 4b is consistent with the absence of multiple hits. Noting that the bottom and left-hand edges of both panels correspond to the actual image volume boundary, we observe that boundary features can be located correctly, though objects with true centers lying outside the image volume are interpreted as small features. It is also clear from comparing panels a and b that some features are apparently missing from the QI. For example, in the region marked A, there is no hint of a large feature which is clearly visible in the corresponding location in PI. Similarly, the two significant features in region B of the real image are absent in panel b. This behavior stems from the loss of resolution in the specular direction, intrinsic to confocal imaging, so that spherical objects appear significantly elongated along the optical z axis. The effect was accounted for by the feature-finder, as demonstrated in Figure 4c, which shows the image remainder after all features were located and removed. Bright regions are those portions of the image that were not identified as belonging to located features and were therefore not removed. In fact, the removal is rather better than this diagnostic image suggests, since the brightness has been scaled to the full range in order to reveal the residue, which for the most part is actually low brightness
portions of the original features. The features missing from A and B in panel b are blanked out in panel c, confirming that the feature-finder has indeed located them. The absence of features in A and B is best understood by considering a projection perpendicular to the optical axis as depicted in Figure 5. In both panels the physical distance scale is isotropic along the two axes, so spherical objects appear as disks in the QI (panel a) but the PI (panel b) clearly exhibits the distorted character of droplets along the optical axis (vertical). The gray arrow indicates the large feature that is missing in region A of the QI in Figure 4b. The gray horizontal line corresponds to the location of the x-y slice shown in Figure 4a, and clearly exhibits an intersection with the smeared “tail” of the large feature (Figure 5b). In the corresponding QI (Figure 5a), the x-y slice does not intersect the reconstructed sphere, so this feature does not appear in this particular plane of the QI. The feature has been located, however, and is present in x-y planes of smaller z values. For completeness, the gray arrow in Figure 4b indicates the plane of the x-z slices of Figure 5. In general, real images appear much “fuller” than the quantified versions because features are strongly smeared in the z direction. Correspondingly, single x-y planes of the QI may not show structure appearing in the PI that contains smeared tails at z above and below the feature sphere. The a priori knowledge that features are spheres is thus essential to this feature-finding procedure. Image slices became dimmer further into the sample, a consequence of scatter which itself is linked to imperfect refractive index matching and the presence of the fluorescent dye. This degradation of image quality with depth into the sample contributed yet another constraint, over and above lens working length and container wall thickness, on how far a sample could be visually penetrated. The deterioration of image contrast must
2012 Langmuir, Vol. 22, No. 5, 2006
Penfold et al.
Figure 7. Droplet size distribution plotted by volume on linear axes, from experiment (solid line) and from the EDM (dotted line). The large-diameter tail has been arbitrarily truncated for the experimental plot.
Figure 6. Same emulsion image quadrant as Figure 3 indicating features used for error evaluation. Those features marked with a star were located to a better than average accuracy by the feature-finder, and the EDM estimates of their radii were closer (between 1 and 2 times) to the “true” radii values as determined using a brightness profiling tool. Identification of those features marked with a white spot was less robust.
be taken into account by the thresholding protocol. In timeseries experiments, photobleaching of the dye is also a consideration. Figure 5b shows a distortion of imaged features to the right with increasing z. This is due to an overall movement in the sample. A slight buoyancy mismatch caused sedimentation at the rate of approximately 10 cm over 129 days. Since an entire stack took about 250 s to acquire, there was sufficient time for regular drift between slices and this is manifest as the inclined blur in Figure 5(b). We have not attempted to correct for this intra-stack drift. Figure 5 also suggests a good identification “strike rate”, where all of the features in the PI are apparently represented in the corresponding QI. In addition, edge features on the boundary of the image volume are correctly located and not biased by a center of mass algorithm. 4.3. Location Accuracy. Errors in the feature location estimates are partitioned into two classes. Class I collects errors arising from imperfect feature location, for example due to the proximity of other features leading to threshold artifacts. Class II errors are stochastic in nature and relate to random effects such as noise. To assess the magnitude of class I errors, the locations of 20 characteristic features determined by the automatic feature-finder have been tested against direct visual inspection of the images. The identified features are shown in Figure 6. To extract “true” values from the RI, a single slice is interrogated using an intensity profiling tool in one dimension that provides pixel coordinates of feature centers and edges by visual inspection. Although useful for the analysis of single features, that are not necessarily isolated, this method is subjective, laborious and resistant to automation in contrast to the EDM scheme. By expressing the error in terms of the RMS quantity dxy ) [∑((∆x)2+(∆y)2)/n]1/2 where ∆x and ∆y are the differences between true and automatically determined coordinates, respectively, and n is the number of features, we find dxy ) 0.7 µm. An estimate of the class I error in z was not attempted since the limited resolution parallel to the optical axis
precluded straightforward estimates of true locations without recourse to additional image processing. “True” values of feature radius, obtained from the onedimensional intensity profile, are also compared with those returned by the feature-finder. Since the EDM automatically returns a radius value REDM not greater than the true radius Rk, we define a multiplicative parameter R > 1 such that Rk ) RREDM. For the same sample of 20 features, the parameter lies in the range 1 e R e 3, with a mean value R j ) 1.56. Good radius estimates are returned for relatively large and clearly distinguished features, while there is also evidence that high values of R stem from poor location predictions. This can be the outcome of a thresholding artifact caused by poorly resolved features, for example. An estimate of class II (stochastic) error can be made by comparing two sets of feature locations and radii extracted from the same sample. With a single microscope and sample, this means that the two sets must be taken at different times, and since the system can be expected to undergo at least a minimal temporal evolution such an estimate is an upper bound on stochastic error. To estimate the level of class II error we have compared the same 20 features from two different image stacks with a 300 s time interval between the start of each stack acquisition. In this time, residual buoyancy mismatching results in a sample-wide sedimentation that has been factored out by estimating the average sample-wide drift over the time interval. The RMS x, y variation between the two stacks is 0.53 µm, and in the z coordinate is 0.44 µm. These figures are an upper bound on stochastic location error. The RMS variation in radius estimate is 0.1 µm. We note that the accuracy of edge feature determination is in principle distinct from that of nonedge features. This is because edge features do not have a location refinement stage based on a brightness center of mass calculation. We reiterate the point that the EDM technique is better than brightness masking at returning realistic estimates of edge feature locations. This is significant since, according to the coordinates of the located features in our sample, one-third of features (1196 out of 3667) qualify as edge features. This is mainly due to the fact that the sample volume measures only four or five “typical” droplet diameters in the z-direction. 4.4. Droplet Size Distribution. Figure 7 compares the droplet size distribution (by volume bins) physically measured by electric zone sensing (Coulter counter) with that obtained by our image quantification scheme, each normalized to give a maximum value of unity. The lower cutoff at ≈1 µm diameter in the distribution from direct visualization corresponds to a minimum feature radius of 3 pixels. Evidently, the sizes extracted by image analysis follow a realistic distribution but are, as expected, displaced toward smaller values due to the inherent tendency of the EDM
QuantitatiVe Imaging of Aggregated Emulsions
Figure 8. Size distribution of located droplets as a function of changing segmentation threshold cut. Open circles (joined by smooth line) correspond to the threshold values tm (m ) 1, 2, 3) used elsewhere in this paper. Diamonds linked by a dotted line indicate the effect of a 10% reduction in t1, the upper of the three thresholds used. Squares linked by a dashed line illustrate the result of increasing t1 by 10%.
to underestimate radii. A scale factor in the range 2-3, as suggested by the direct analysis of section 4.3, would serve to approximately superimpose the two curves. Uniform scaling is a crude approximation, however, while a factor of this magnitude is too large, as an analysis of the volume fraction reveals. The experimentally determined volume fraction is φ ) 30%. After excluding background noise, a simple pixel count in the PI gives φ ) 61%, an upper bound determined by direct visualization. This apparent inflation arises because each truly spherical feature is smeared in the z direction, and thus contributes a larger bright region than expected. By assuming that the overestimate is proportional to the z distortion of the droplet radius, then the aspect ratio should be ξ ) 61/30 ) 2.03. Smeared features will overlap, however, for droplets closely situated in z, so the true ξ should be larger than the ratio of the pixel counts. Inspection of x-z image slices suggests a value 2.6 < ξ < 2.9, consistent with the voxel interpolation used in the restoration algorithm. By applying a uniform radius scaling factor R j ) 1.56, estimated manually from a sample of one-dimensional intensity profiles (section 4.3), then the bright pixel count of the reconstructed image corresponds to φ ) 25.5%. The droplet size distribution presents one way of indicating the impact of different choices of threshold value. For the data set discussed in section 4.2, Figure 8 shows the change in the droplet number distribution when the first and largest of the three selected threshold values t1 is increased or decreased by 10% relative to the preferred value tref 1 . All remaining parameters are held fixed. With t1 ) 1.1tref , the total number of particles 1 increases by 4.1%, while setting t1 ) 0.9tref 1 decreases the count by 11.4%. Most of the change is to be found in the lower size classes; as expected since decreasing t1 shifts the size distribution to smaller droplet diameters. Scaling all three threshold values by 10% produces similar changes to those seen by changing just t1 alone. A threshold change of 10% typically corresponds to 20 units for bright features in an 8-bit image. An entirely manual process, combining subjective thresholding with one-dimensional brightness profiling, is unlikely to select values outside this window. At the lower limit, feature differentiation deteriorates as objects of the BI appear to merge, whereas increasing numbers of objects are “top-sliced” or missed altogether at the upper limit. An inappropriate choice of threshold only weakly perturbs feature position coordinates but impacts significantly on radius estimates and, more seriously, on feature identification. 4.5. Void and Pore Morphology. The generalized Apollonius construction was used to analyze the pore structure of the moderately concentrated emulsion depicted in Figure 3a. Since
Langmuir, Vol. 22, No. 5, 2006 2013
Figure 9. Size distribution of O-holes (s) and I-holes (---) are compared for an aggregated emulsion at 30% volume fraction.
Figure 10. Void space construction using the generalized Appollonius algorithm is shown for a subsample from an image stack of an emulsion gel at 30% volume fraction where 2798 particles are identified (panel a). The complementary space of overlapping holes is included in panel b and shown alone in panel c. Panel d shows the subset of non-overlapping independent holes (gray), whereas the isolated holes are also indicated (black).
the computational effort scales such as N4 for a system of N particles, it was necessary to form a subsample by extracting the central prism that occupies 80% of the imaged volume and encloses 2798 particles. Size distributions for the overlapping and independent holes are compared in Figure 9 where the respective average diameters are 13.3 and 10.7 µm. With its mode significantly displaced below the expected value, the O-hole distribution is noticeably more asymmetric and there is some evidence for a broad shoulder straddling the mean. Although the I-hole distribution appears more strongly centered, some left skewness remains and the variance is somewhat larger, as expected. The slowly decaying tail, extending to a few instances of holes more than four times larger than average, is inherited from the underlying O-hole space. The corresponding structure is shown explicitly in Figure 10, where the particle network is indicated in panel a. Panel b exemplifies the utility of the generalized Apollonius method by demonstrating the high complementarity of the particle configuration and the O-hole morphology. To clarify, the O-hole structure alone is shown in panel c while the subset of nonoverlapping I-holes is seen in panel d. A cluster analysis on the same subsample is detailed in Table 1 and clearly confirms
2014 Langmuir, Vol. 22, No. 5, 2006
Penfold et al.
Table 1. Cluster Analysis Results for a Subsample from an Image Stack of an Emulsion Gel at 30% Volume Fraction Where 2798 Particles Are Identifieda no. holes/cluster no. clusters no. holes
1 79 79
2 38 76
3 8 24
4 5 20
9 1 9
10 1 10
12 1 12
1282 1 1282
a The generalized Apollonius algorithm found 1512 overlapping holes, of which 464 are independent, with 79 isolated holes.
Figure 12. Three dimensional metric distributions of an aggregated 30% o/w emulsion from Voronoi analysis of droplet coordinates obtained by quantitative imaging. The quantity M denotes polyhedral perimeter (s), surface area (---), and volume (‚‚‚), respectively and the abscissa is normalized by the corresponding average value 〈M〉.
provide useful characteristic relations. For packings built numerically by random sequential adsorption, Oger et al.43 determined that K1 and K2 were independent of volume fraction and approximately equal to 0.75 and 14.3, respectively. On the other hand, Yang et al.42 reported that DEM generated packings exhibit a strong density dependence in these parameters, with K1 increasing and K2 decreasing as volume fraction grew. For the 30% volume fraction emulsion gel studied here, somewhat smaller values of these parameters are observed with K1 ) 0.558 and K2 ) 13.7.
5. Conclusions Figure 11. Topological measures of an aggregated 30% o/w emulsion from Voronoi analysis of droplet coordinates obtained by quantitative imaging.
percolation of the void space at this volume fraction with a single box spanning cluster of 1282 O-holes. At the opposite extreme of isolated holes (unit cluster size) and small clusters, the distribution appears to decay geometrically. The subset of isolated holes is also indicated in Figure 10(d). Following Yang et al.,42 we distinguish between topological characteristics of the Voronoi tessellation, in particular: 1. the number of edges for each polyhedron face; and 2. the number of faces for each polyhedron; as well as metric quantities: 3. the perimeter and area of a polyhedron face; and 4. the perimeter, area and volume of a polyhedron. Figure 11 exhibits the topological measures of the 30% volume fraction emulsion corresponding to Figure 3a following analysis of droplet coordinates obtained by quantitative imaging. The results are entirely consistent with computer simulation of sedimenting random packings using the discrete element method (DEM) and incorporating long-range van der Waals forces.42 This is also confirmed by the analysis of 3D metric quantities shown in Figure 12 where the polyhedral perimeter distribution appears highly symmetric, whereas the surface area and volume distributions are progressively more right-skewed with modes significantly smaller than the corresponding mean values. Given the average perimeter 〈P〉, surface area 〈S〉, and volume 〈V〉 of the Voronoi polyhedra, both the sphericity coefficient
36π〈V〉2
K1 )
〈S〉3
(11)
and the parameter
K2 )
1 4π 2/3 〈P〉 3 3 〈V〉1/3
( )
(12)
Quantitative 3D imaging of aggregated emulsions is demonstrated using a novel formulation designed specifically for good matching of both refractive index and buoyancy. The refractive index matching is necessary to reduce scatter and allow optical penetration. The buoyancy matching is essential to slow gravity-induced overall drift to such an extent that microscopic imaging becomes viable. The conflicting constraints of image slice acquisition time, optical resolution, lens working distance and residual light scattering have been set against the need for high resolution over the largest possible image volume in the shortest possible time. For our microscope system, an image volume large enough to embrace droplet aggregates in a time frame short enough that intra-stack motion is modest results in a compromised z resolution. Samples were imaged along a horizontal optical z axis, perpendicular to the gravitational field, rather than in the more conventional vertical or inverted modes. To maximize the image quality, we have successfully implemented a 3D deconvolution algorithm, based on an analytic optical transfer function, without recourse to small angle approximations that are inappropriate for a high numerical aperture lens. This is distinct from simple noise filtering. Since the microscope resolution is finite, however, information is lost by the image formation system. The impact of this information loss and of the high-frequency noise on the Fourier inversions has been mitigated by using a regularizing scheme based on a linear filter. A feature-finding algorithm has also been developed, based on the Euclidean distance map, for both locating polydisperse, aggregated features in 3D and for estimating their sizes. The technique relies on the a priori knowledge that the features are spheres, and exploits the fact that local maxima of the EDM provide a measure of feature size and position. Due attention must be paid to the impact of incomplete confocality and relatively large z-step values, that both compromise resolution along the (42) Yang, R. Y.; Zou, R. P.; Yu, A. B. Phys. ReV. E. 2002, 65, 041302. (43) Oger, L.; Gervois, A.; Troadec, J. P.; Rivier, N. Philos. Mag. B 1996, 74, 177.
QuantitatiVe Imaging of Aggregated Emulsions
optical axis. The prescription works best for separated, clearly imaged features but gives acceptable results, with quantifiable errors, for aggregated systems. The technique offers a significant improvement over simple brightness masking for features near the boundary of the image volume. In passing, we have noted that the EDM followed by thinning may offer a route to sophisticated structural analysis in 3D. The resulting skeletons may provide a means of mapping inter-droplet contacts and thereby prove useful in the analysis of stress chains. This would give quantitative insight into a number of properties such as droplet connectivity and coordination number, as well as aggregate size and orientation. No attempt is made to develop this theme in the present work. Having demonstrated that the EDM approach can be used to extract quantitative information from CLSM images, we offer a number of suggestions for improvements. First, the deconvolution stage could be more fully integrated into the voxel interpolation step and use of the a priori knowledge to objectively fix parameters of the regularisation operator. In this way, the need for subjective manual intervention could be reduced and the complete image processing cycle rendered automatic with
Langmuir, Vol. 22, No. 5, 2006 2015
clear benefits for handling large data sets, especially relevant for tracking dynamical behavior. Second, where we have used Gaussian smoothing an alternative that better preserves feature boundary sharpness might be appropriate. Third, and most importantly, we believe there is scope for an improved thresholding protocol. Explicit droplet coordinates in 3D allows structural analysis at a level that has previously been impossible. To illustrate the potential, we have used the coordinates of 3667 droplets extracted from a 3D image to study the void morphology using a novel generalization of the Apollonious problem and conventional Voronoi tessellation. This work also promises the direct measurement of many-body correlation functions, particularly beyond the pairwise level accessible to scattering experiments. Acknowledgment. This work is funded by the UK Biotechnology and Biological Sciences Research Council, under Grant 218/D17326. We are grateful to Grant Calder and the John Innes Centre for making the confocal microscope available to us. Insightful discussion with Gary Barker is also acknowledged. LA052719W
