## Abstract

Computerized image-analysis routines deployed widely to locate and track the positions of particles in microscope images include several steps where images are convolved with kernels to remove noise. In many common implementations, some kernels are rotationally asymmetric. Here we show that the use of these asymmetric kernels creates significant artifacts, distorting apparent particle positions in a way that gives the artificial appearance of orientational crystalline order, even in such fully-disordered isotropic systems as simple fluids of hard-sphere-like colloids. We rectify this problem by replacing all asymmetric kernels with rotationally-symmetric kernels, which does not impact code performance. We show that these corrected codes locate particle positions properly, restoring measured isotropy to colloidal fluids. We also investigate rapidly-formed colloidal sediments, and with the corrected codes show that these sediments, often thought to be amorphous, may exhibit strong orientational correlations among bonds between neighboring colloidal particles.

© 2013 Optical Society of America

Investigations of dense particle systems have yielded significant insight into the physics of collective phenomena, including crystallization [1], glass transitions [2–6], gelation [7], granular packing [8–10], and nonequilibrium self-organization [11], in both 2D and 3D [12]. Most of these studies focus on colloids, micron-sized particles suspended in a fluid, which are small enough that thermodynamics governs their behavior, yet large enough to be resolved with optical techniques. Optical microscopy, in particular laser-scanning confocal microscopy, can image many thousands of such particles simultaneously and at high speeds, providing information on structure and dynamics of individual particles in large systems. [1,2,4–7,10–12].

Typical microscopy-based colloid studies collect large numbers of images, then process them with algorithms to determine the location of the particles. The most widely used, foundational codes were first created by John C. Crocker and David G. Grier (CG) more than 15 years ago [13], with some subsequent improvements in accuracy and performance [14–17]. These codes have had an enormous positive impact within the fields of soft-matter physics and biology; however, they were developed and optimized using languages and computer architectures that have since been superseded; in particular, some of the algorithmic choices made in the original code, based on then-current technologies, may no longer be ideal.

In this paper, we show that the rotationally-asymmetric convolution kernels for image processing within the standard particle-location codes can create significant distortions in the positions of located particles, in ways that create artifactual appearence of crystalline order even in physically isotropic systems, such as a neutrally-buoyant, dense fluid suspension of diffusing hard-sphere colloidal particles. The standard codes measurably distort not only the angles of bonds between particles, but also the radial distribution function *g*(*r*), in samples with high particle volume fraction *ϕ*. These artifacts cannot be removed by varying the input parameters, nor by derivatives of the algorithm which employ similar asymmetric kernels. However, by using rotationally-symmetric kernels in all stages of the image processing pipeline, we restore quantitatively-correct values known from theory for hard-sphere systems; our changes neither modify the algorithm globally, nor affect execution performance on modern architectures. We use the corrected algorithm [17] to explore the positions of particles in randomly-packed sediments of colloidal spheres; intriguingly, we find that when these sediments are prepared by centrifugation in thin capillaries, the bonds between neighboring particles adopt a strong hexagonal order, which is not expected for these ostensibly amorphous packings.

We suspend colloidal hard spheres of polymethylmethacrylate, sterically stabilized by poly-12-hydroxystearic acid, in a solvent mixture of 1:2.1 (by mass) cis/trans decahydronaphthalene and tetrachloroethylene [10]; this solvent mixture closely matches the density of the particles, which do not sediment appreciably over several weeks, and their refractive index, so that confocal microscopy can resolve individual particles more than 80 *μ*m into the sample. The particles have a diameter *σ* = 2.4 ± .05 *μ*m and are labeled with Nile Red; their diameters have poly-dispersity less than 0.05, measured with a combination of static and dynamic light scattering, confocal microscopy, and scanning electron microscopy of dry particles under vacuum. For the colloidal fluids, we load suspensions into capillaries (Vitrocom) of internal size 0.1×2×5 mm and seal with 5-minute epoxy; for colloidal sediments, we resuspend the particles in decalin alone before loading, and centrifuge the sealed capillary at 3000 rcf, aligning the capillary’s long axis along the direction of effective gravity.

We image these samples with a laser-scanning confocal fluorescence microscope (Nikon A1R) and 100× 1.4NA oil immersion objective (Nikon) using fast resonant scanning (30 fps) for colloidal fluids, and slow scanning (0.1–1 fps) for the sediments. We excite with *λ* = 514 nm and collect fluorescence images whose pixels correspond to either 0.166 or 0.12 *μ*m in the sample. Spatial resolutions, determined by the confocal pinhole, are 1.6 and 0.15 *μ*m along and perpendicular to the optical axis, respectively. We orient the capillary’s long axis along the fast resonantly-scanned axis, corresponding to the horizontal image axis, unless noted otherwise.

We show a typical image of a disordered fluid in Fig. 1(a). Our frame acquisition time of 0.033 sec is far faster than the time for the colloidal spheres to diffuse their own size: ∼10 sec in the dilute limit; far longer in our dense samples with *ϕ* = 0.39 ± 0.04. Therefore, we collect many slices through each particle at different depths before it moves significantly, permitting accurate estimation of each particle’s 3D position with suitable software, particularly the pioneering CG method [13], and subsequent derivatives engineered to improve usability, performance [14,16,17], and, in principle, accuracy [15]. We use primarily a C++ implementation (PLuTARC) [14,17], though our results do not change appreciably with other codebases.

The first stage of the CG image processing procedure [13] involves a “bandpass” filter, convolving in real-space the original image *I* with two sets of kernels: a small-wavelength gaussian kernel *G*, to remove high-frequency noise; and a much larger square kernel *T*_{sq} [18] to remove long-wavelength fluctuations in background brightness, illustrated in the inset to Fig. 1(b); the square kernel shape was originally chosen for performance reasons when the CG method was developed [19]. The results of these two operations are combined to yield the filtered image *I*_{f}*= I* * *G* − *I* * *T*_{sq}, where the * operator denotes 2D realspace convolution. This procedure preserves intensity variations on the order of the colloid size; the centers are then located in the second stage of CG algorithm [13]. We apply this standard CG analysis using *T*_{sq} to the image in Fig. 1(a), marking the center of each particle’s located position with a white dot in Fig. 1(b).

From these 2D positions, we calculate the standard radial distribution function *g*(*r*) [20], shown with green squares in Fig. 2(a); as expected, *g*(*r*) has a peak at *r* = *σ*, oscillates with the shell structure of the fluid [20] at intermediate values of *r*, and asymptotes to 1 for *r* ≫ *σ*. To evaluate the accuracy of particle location algorithm, we calculate *g*(*r*) for a hard-sphere fluid at *ϕ* = 0.39 by solving the Ornstein-Zernike equation [20] with the Percus-Yevick (PY) approximation [21], and plot this theoretical *g*(*r*) with a solid grey curve in Fig. 2(a). Strikingly, the overall agreement between the theoretical prediction and the experimental data, processed with the standard CG algorithm, is extremely poor. The deviation for *r* < *σ* is a known, small artifact: occasional overlaps in the projected images of closely-approaching spheres in a single plane [22] increase the apparent density of particles separated by less than the contact distance.

By contrast, the discrepancy at intermediate *r* cannot be explained away. The two-dimensional density distribution, which when averaged azimuthally yields *g*(*r*), is highly anisotropic and exhibits a square symmetry not physically present in the sample, as shown in the inset to Fig. 2(a). This anisotropy suggests that the discrepancy in *g*(*r*) could be caused by a step that breaks the circular symmetry and introduces a square ordering, e.g. the *T*_{sq} filtering. Indeed, varying the size of *T*_{sq} changes *g*(*r*) significantly, but in no case leads to even reasonable agreement with the theoretical prediction, as shown in Fig. 2(a). When *T*_{sq} is large, the first peak is low, yet centered around *r* = *σ*; therefore, effects that would change the average interparticle separation, such as electric charge repulsion, are insignificant. Shrinking *T*_{sq} raises the first peak’s height, but then the intermediate-*r* oscillations do not match the theoretical *g*(*r*).

To explore whether the cause of the *g*(*r*) deviations also affects local orientational order, we identify the nearest-neighbors (NN) of each particle with the Delaunay triangulation [23], and determine the probability *P*(*θ*) for a particle to have a NN bond at an angle *θ* relative to the long capillary axis (i.e. the horizontal image axis). *P*(*θ*) is normalized by its average; we ignore particles near image edges to remove boundary effects. For our disordered dense colloidal fluid sample, *P*(*θ*) should be circular: all of the NN bonds should be distributed isotropically through all *θ*. Strikingly, however, *P*(*θ*), calculated from the particle positions located in *T*_{sq}-filtered images, is highly anisotropic; it is not circular, but instead has a dramatic four-lobed appearance with peaks at *θ* = 0°, 90°, 180° and 270°, shown with the red dashed curve in Fig. 3(a). We quantify this anisotropy with the difference between the maximum and minimum of *P*(*θ*), which should be 0 for an isotropic system; here we find a value of 0.36.

These data suggest that what should be a fully disordered sample might actually have strong orientational order. To test this possibility, we physically rotate the fluid sample by 45° on the microscope, so that the direction corresponding to long axis of the capillary is at an angle of *χ* = 45° relative to the horizontal image axis. Since any forces outside of the capillary are unlikely to perturb the fluid, rotation of the capillary with respect to the lab frame should rotate the lab-frame orientation of the NN bonds; that is, if the orientational order is a physical attribute of the sample, *P*(*θ*) should concertedly rotate by 45°, with lobes at *θ* = 45°, 135°, 225°, and 315°. Conversely, if the orientational anisotropy is not physically present, but instead introduced by some aspect of the particle-location algorithm, physical sample rotation would leave *P*(*θ*) unchanged. Indeed, *P*(*θ*) for a sample rotated to *χ* = 45°, shown with a black curve in Fig. 3(a), perfectly overlaps with the original *P*(*θ*) obtained at *χ* = 0°. These data demonstrate that the appearance of orientational anisotropy in our physically-isotropic disordered fluid is an artifact of how the particles are located with software, not an actual physical characteristic.

To determine precisely when these artifacts are introduced, we rotate the experimental data by 45° at each stage of image collection and processing, and determine the resulting *P*(*θ*). We find that rotation at any stage before the *T*_{sq} convolution leaves *P*(*θ*) unchanged, as shown by the blue dotted curve in Fig. 3(a); by contrast, rotating the *T*_{sq}-filtered image by 45° causes *P*(*θ*) to rotate by the same amount, as shown with the green curve in Fig. 3(a). Therefore, it must be the *T*_{sq} convolution that introduces the orientational anisotropy. Removing artifacts of finite pixel size using the recent CG variant “fracshift” [15] does not change the cross-shaped appearance of *P*(*θ*) appreciably, shown by the purple curve in Fig. 3(b).

To explore the effects of shape of the background-subtraction kernel, we replace the square *T*_{sq} by a rotationally-symmetric circular pillbox kernel *T*_{cir} [18], shown in the inset to Fig. 1(c). The filtered image and located particle positions appear similar to the eye, as shown in Figs. 1(b) and 1(c); however, the differences in the structural metrics derived from the particle positions are strikingly large: filtering with the rotationally-symmetric *T*_{cir} kernel in place of *T*_{sq} restores rotational symmetry to both *P*(*θ*), shown with filled circles in Fig. 3(c), and the two-dimensional density distribution, shown in the inset to Fig. 2(b). Moreover, the resulting *g*(*r*) matches closely the peak, oscillations and asymptote of the theoretical PY prediction at the correct *ϕ*, as shown with red symbols and solid grey curve in Fig. 2(b); using the experimental *σ*, there are no free fitting parameters. This close agreement demonstrates that Coulombic charge has no significant effect in our particle suspensions [1, 10, 24], contrasting its role in fluids of colloidal ellipsoids [25, 26]. Together, these data demonstrate that the broken orientational symmetry is an artifact introduced by *T*_{sq}: physically, our sample is orientationally isotropic.

Although changing the size of *T*_{sq} does not fix the discrepancy in *g*(*r*), as previously shown, a sufficiently large *T*_{sq}, averaging over a sufficiently large number of particles, might restore isotropy. To test this, we first reduce the relative particle size (in pixels) by half with a 2 × 2 binning of the image; the resulting *P*(*θ*) is essentially unchanged. Furthermore, we repeat the analysis of the original images with far larger *T*_{sq} kernels, corresponding to sizes of tens of microns. We find that some anisotropy in *P*(*θ*) always persists; moreover, these much larger kernels create other significant problems: when the *T*_{sq} size exceeds the average interparticle spacing of a few microns, the software fails to locate all particles; thus, metrics derived from these positions may not be physically correct. Consequently, the *P*(*θ*) artifact appears to be a fundamentally inherent consequence of geometry: the diagonal of a square is longer than its side, so that the *I* * *T*_{sq} convolution biases those particles with NNs in a diagonal direction. Because this term is subtracted off to yield *I*_{f}, particles appear more likely to have NNs in the horizontal and vertical directions, evidenced by the *P*(*θ*) lobes in Figs. 3(a) and 3(b).

We choose the circular pillbox kernel to demonstrate the importance of rotationally-symmetric kernels, because *T*_{cir} is the closest circularly-symmetric kernel to *T*_{sq}, but in principle other circular kernels might also work. To test this, we repeat the analysis on the same image data shown in Fig. 1(a), using a 2D circular gaussian kernel *T*_{gau} [18], shown in the inset to Fig. 1(d), of similar size to *T*_{cir}. The filtered image and particle positions are quite similar to those derived using *T*_{cir}, as shown in Figs. 1(c) and 1(d); moreover, the *g*(*r*) and *P*(*θ*) curves are nearly identical, as shown in Figs. 2(b) and 3(c), respectively. These data demonstrate nearly identical, physically-correct results by filtering with rotationally-symmetric kernels, irrespective of specific functional form; our results are easily generalized to any well-behaved, rotationally-symmetric kernel, which can be expressed as a suitable linear combination of pillbox and gaussian basis functions. Furthermore, the gaussian has significant performance advantages over the pillbox, as it is separable: exp(−*r*^{2}) = exp(−*x*^{2} − *y*^{2}) = exp(−*x*^{2})exp(−*y*^{2}); 2D gaussian convolution can be carried out as successive 1D convolutions, which are substantially faster [18]. Thus, in addition to being more accurate, removing backgrounds via the two 1D gaussians should incur no performance penalty relative to the existing algorithms.

While rotationally-symmetric kernels are absolutely necessary for quantitative structural analysis of dense, disordered fluid samples (such as in Fig. 2), the symmetry of the background-subtraction kernel may play a different role in less-dense or strongly-ordered systems. To explore this, we repeat the experiment and analysis for a fluid at lower *ϕ* = 0.05, and a hexagonal close-packed crystalline layer. For the low-density fluid, *T*_{sq}-filtering yields a less anistropic *P*(*θ*), as shown with symbols in Fig. 3(b); for the crystalline layer, *P*(*θ*) is a very strong, six-pointed star pattern that reflects its strong orientational order, as shown with the red symbols in Fig. 3(d). *T*_{cir}-filtering of the crystalline sample yields a qualitatively-similar figure, as shown with blue symbols in Fig. 3(d); however, significant quantitative differences remain: the *T*_{sq}-filtered data are anisotropic, with slightly stronger peaks at *θ* = 0° and 180°, and an overall standard deviation of 9% for all six peaks; by contrast, data filtered with rotationally-symmetric kernels show no such bias, and the standard deviation for *T*_{gau} is a far smaller 3%.

*Our results therefore demonstrate that, to locate particles accurately, only rotationally-symmetric kernels no larger than the average interparticle spacing should be used.*

Finally, using these rotationally-symmetric, properly-sized kernels, we investigate dense, disordered systems where orientational anisotropy of NN bonds may have important physical consequences, such as the solid non-crystalline colloidal packings used to model the behavior of glasses [2]. We prepare these solids using a centrifuge to sediment rapidly the colloidal particles in a capillary [10], and image with confocal microscopy, as shown in Fig. 4(a). In this system, the effective gravity field breaks spatial symmetry in a way fundamentally different from atomic and molecular glasses, which also lack the hydrodynamic interactions between sedimenting colloids and the wall. Interestingly, we observe that the orientations of NN bonds exhibit significant rotational anisotropy: the bonds are more likely to point along the axis of sedimentation or at an integer multiple of 60° with respect to this direction, as shown with a solid blue curve in Fig. 4(b). To confirm that the obtained six-fold symmetric *P*(*θ*) is due to a true breaking of symmetry within the sample, not an imaging artifact, we physically rotate our sample to *χ* = 30°; this drives a similar rotation of *P*(*θ*), demonstrating that the ordering is real, as shown in Fig. 4(b). The observed orientational anisotropy of NN bonds, manifest as six-fold symmetry of *P*(*θ*), with one of the lobes oriented along the axis of (effective) gravity, differentiates our sediments from other amorphous solids, such as the molecular glasses, where no external fields usually exist to break the global symmetry of the structure, so that all spatial orientations are equivalent; a full physical mechanism of the observed breaking of rotational isotropy in sediments is currently being investigated.

## Acknowledgments

We are very grateful for helpful discussions with David Grier (NYU), John Crocker (UPenn), Eric Weeks (Emory), David Weitz (Harvard), Moshe Schwartz (Tel Aviv) and Yuwen Wang (Swarthmore). We thank Andrew Schofield (Edinburgh), for colloid synthesis; and Erez Janai (Bar Ilan), for code. M. S., E. S. and A. V. B. thank the Kahn Foundation for equipment, and the Israel Science Foundation ( #85/10, #1668/10) for financial support; P. J. L. thanks NASA ( NNX08AE09G, NNX08AE09G S11, NNC08BA08B, NNX13AQ48G), the NSF ( DMR-1310266), and the Harvard MRSEC ( DMR-0820484) for financial support.

## References and links

**1. **U. Gasser, E. R. Weeks, A. B. Schofield, P. N. Pusey, and D. A. Weitz, “Real-space imaging of nucleation and growth in colloidal crystallization,” Science **292**, 258–262 (2001). [CrossRef] [PubMed]

**2. **A. van Blaaderen and P. Wiltzius, “Real-space structure of colloidal hard-sphere glasses,” Science **270**, 1177–1179 (1995). [CrossRef]

**3. **E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, “Three-dimensional direct imaging of structural relaxation near the colloidal glass transition,” Science **287**, 627–631 (2000). [CrossRef] [PubMed]

**4. **Z. Zhang, N. Xu, D. T. N. Chen, P. Yunker, A. M. Alsayed, K. B. Aptowicz, P. Habdas, A. J. Liu, S. R. Nagel, and A. G. Yodh, “Thermal vestige of the zero-temperature jamming transition,” Nature **459**, 230–233 (2009). [CrossRef] [PubMed]

**5. **Z. Zheng, F. Wang, and Y. Han, “Glass transitions in quasi-two-dimensional suspensions of colloidal ellipsoids,” Phys. Rev. Lett. **107**, 065702 (2011). [CrossRef] [PubMed]

**6. **G. L. Hunter and E. R. Weeks, “The physics of the colloidal glass transition,” Rep. Prog. Phys. **75**, 066501 (2012). [CrossRef] [PubMed]

**7. **P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino, and D. A. Weitz, “Gelation of particles with short-range attraction,” Nature **453**, 499–503 (2008). [CrossRef] [PubMed]

**8. **T. Aste, M. Saadatfar, and T. J. Senden, “Geometrical structure of disordered sphere packings,” Phys. Rev. E **71**, 061302 (2005). [CrossRef]

**9. **M. Jerkins, M. Schröter, H. L. Swinney, T. J. Senden, M. Saadatfar, and T. Aste, “Onset of mechanical stability in random packings of frictional particles,” Phys. Rev. Lett. **101**, 018301 (2008). [CrossRef]

**10. **S. R. Liber, S. Borohovich, A. V. Butenko, A. B. Schofield, and E. Sloutskin, “Dense colloidal fluids form denser amorphous sediments,” Proc. Nat. Acad. Sci. U. S. A. **110**, 5769–5773 (2013). [CrossRef]

**11. **X. Cheng, J. H. McCoy, J. N. Israelachvili, and I. Cohen, “Imaging the microscopic structure of shear thinning and thickening colloidal suspensions,” Science **333**, 1276–1279 (2011). [CrossRef] [PubMed]

**12. **U. Gasser, “Crystallization in three- and two-dimensional colloidal suspensions,” J. Phys. Condens. Mat. **21**, 203101 (2009). [CrossRef]

**13. **J. C. Crocker and D. G. Grier, “Methods of digital video microscopy for colloidal studies,” J. Colloid Interface Sci. **179**, 298–310 (1996). [CrossRef]

**14. **P. J. Lu, P. A. Sims, H. Oki, J. B. Macarthur, and D. A. Weitz, “Target-locking acquisition with real-time confocal (TARC) microscopy,” Opt. Express **15**, 8702–8712 (2007). [CrossRef] [PubMed]

**15. **Y. Gao and M. L. Kilfoil, “Accurate detection and complete tracking of large populations of features in three dimensions,” Opt. Express **17**, 4685–4704 (2009). [CrossRef] [PubMed]

**16. **http://tacaswell.github.io/tracking/html/

**17. **http://github.com/peterlu/PLuTARC_centerfind2D

**18. **S. W. Smith, *The Scientist & Engineer’s Guide to Digital Signal Processing* (California Technical, 1997).

**19. **D. G. Grier and J. C. Crocker, personal communication.

**20. **J.-P. Hansen and I. R. McDonald, *Theory of Simple Liquids* (Elsevier, 2006).

**21. **D. Frenkel, R. J. Vos, C. G. de Kruif, and A. Vrij, “Structure factors of polydisperse systems of hard spheres: A comparison of Monte Carlo simulations and Percus-Yevick theory,” J. Chem. Phys. **84**, 4625–4630 (1986). [CrossRef]

**22. **D. R. Wilkinson and S. F. Edwards, “The use of stereology to determine the partial two-body correlation functions for hard sphere ensembles,” J. Phys. D Appl. Phys. **15**, 551–562 (1982). [CrossRef]

**23. **M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars, *Computational Geometry: Algorithms and Applications* (Springer, 2008).

**24. **C. P. Royall, M. E. Leunissen, A.-P. Hynninen, M. Dijkstra, and A. van Blaaderen, “Re-entrant melting and freezing in a model system of charged colloids,” J. Chem. Phys. **124**, 244706 (2006). [CrossRef] [PubMed]

**25. **A. P. Cohen, E. Janai, E. Mogilko, A. B. Schofield, and E. Sloutskin, “Fluid suspensions of colloidal ellipsoids: Direct structural measurements,” Phys. Rev. Lett. **107**, 238301 (2011). [CrossRef] [PubMed]

**26. **A. P. Cohen, E. Janai, D. C. Rapaport, A. B. Schofield, and E. Sloutskin, “Structure and interactions in fluids of prolate colloidal ellipsoids: Comparison between experiment, theory, and simulation,” J. Chem. Phys. **137**, 184505 (2012). [CrossRef] [PubMed]