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.

 

Fig. 1 (a) Unprocessed confocal microscope image of a dense fluid at ϕ = 0.39 of colloidal spheres (σ =2.4 μm; 19 pixels). (b)–(d) Images filtered via convolution with several kernels, whose intensity maps are inset (not to scale): (b) square T_sq (28 pixels), (c) circular pillbox T_cir (28 pixels) and (d) circular gaussian T_gau (24 pixels). After filtering, the same particle-location algorithm is applied; final particle positions are marked with white dots.

Download Full Size | PDF

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.

 

Fig. 2 (a) g(r) derived from particle positions, such as in Fig. 1(b), using T_sq of size 2.3 μm (green squares). Theoretical g(r) calculated for a fluid of hard spheres at the same ϕ (solid grey curve) is completely different; changing the size of T_sq (blue and pink squares) does not improve the agreement. Inset: 2D density distribution function is highly anisotropic and shows an unphysical square symmetry. (b) By contrast, filtering with a rotationally-symmetric kernel, either T_cir (red circles) or T_gau (blue circles) leads to very close agreement for r > σ. Inset: 2D density distribution function is circularly symmetric.

Download Full Size | PDF

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.

 

Fig. 3 (a) P(θ) of disordered particles in the T_sq-filtered image in Fig. 1(b) has a cross-like shape, independent of sample orientation, either along the x̂-axis of the confocal scanner (χ = 0°, red dashed curve) or physically rotated by χ = 45° (black curve). P(θ) is insensitive to a 45°-image rotations before T_sq filtering (blue dotted curve), yet rotates with the image when this rotation follows T_sq filtering (green curve). (b) P(θ) after refining particle positions with fracshift [15] (purple curve) has the same artifact; anisotropy is less pronounced for lower ϕ =0.05 (orange circles). (c) P(θ) after filtration with rotationally-symmetric kernels, either T_cir (black circles) or T_gau (white circles), is correctly isotropic. (a)–(c) are at the same scale. (d) P(θ) for a crystalline sample, filtered with T_sq (red squares) vs. T_cir (blue circles) and T_gau (blue triangles), show smaller differences.

Download Full Size | PDF

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²) = exp(−x² − y²) = exp(−x²)exp(−y²); 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.

 

Fig. 4 (a) Confocal microscope image of a colloidal sediment prepared in a capillary with a centrifuge; effective gravity direction marked by the arrow. (b) P(θ) for this sediment demonstrates preferential bond orientation along the direction of gravity (θ = 0) and at integer multiples of 60° with respect to this direction; data points (blue circles) have uncertainty of ±0.015 and B-spline (solid curve) is included as a guide to the eye. Physically rotating the sample by χ = 30° on the microscope stage, then rotating the resulting P(θ) by −30° (red squares and curve), follows the original data closely.

Download Full Size | PDF