Supporting Information-Determination of the positions and orientations of concentrated rod-like colloids from 3 D microscopy data

Confocal microscopy in combination with real-space particle tracking has proven to be a powerful tool in scientific fields such as soft matter physics, materials science and cell biology. However, 3D tracking of anisotropic particles in concentrated phases remains not as optimized compared to algorithms for spherical particles. To address this problem, we developed a new particle-fitting algorithm that can extract the positions and orientations of fluorescent rod-like particles from three dimensional confocal microscopy data stacks, even when the fluorescent signals of the particles overlap considerably. We demonstrate that our algorithm correctly identifies all five coordinates of uniaxial particles in both a concentrated disordered phase and a liquid-crystalline smectic-B phase. Apart from confocal microscopy images, we also demonstrate that the algorithm can be used to identify nanorods in 3D electron tomography reconstructions. Lastly, we determined the accuracy of the algorithm using both simulated and experimental confocal microscopy data-stacks of diffusing silica rods in a dilute suspension. This novel particle-fitting algorithm allows for the study of structure and dynamics in both dilute and dense liquid-crystalline phases (such as nematic, smectic and crystalline phases) as well as the study of the glass transition of rod-like particles in three dimensions on the single particle level.


INTRODUCTION
Colloidal particles are applied throughout industry, for example in paints, personal care products, food, ceramics and pharmaceutics [1][2][3]. Additionally, they are also applied in recent commercially available products such as the electronic ink in e-readers [4]. As a result, the characterisation of the structure and dynamics of colloidal suspensions is important for many industrial applications. Furthermore, hard-sphere colloidal suspensions have proven to serve as a model system to investigate phenomena such as crystallization, the glass transition and flow induced behaviour on the single particle level [5][6][7][8][9][10][11][12]. In many of these studies, an image processing technique was applied based on the algorithm described by Crocker and Grier [13]. In their algorithm, spherical particles are located in 2D digital microscopy images using a local brightness maxima criterion. The position is refined by calculating the brightness-weighted centroid of a cluster of pixels. This method was extended to 3D either slice-by-slice [5,7] or by considering the full 3D image [9]. Crocker and Grier also reported a method to obtain the trajectories of the individual particles in time, known as particle tracking [13]. Since then, there have been numerous algorithms that locate or track spherical particles with increased accuracy or performance [14][15][16][17][18][19][20]. These extensions and alternatives are all based on processing images of spherical particles. However, due to recent progress in particle synthesis, well-defined (shape) anisotropic colloids are becoming widely available, see e.g. Refs. 21-26. These particles can often be observed directly with a (confocal) microscope and therefore enable quantitative measurement of not only their positional but also their rotational degrees of freedom. Therefore, a rapid increase in the number of algorithms that extract coordinates of anisotropic particles from microscopy images has taken place [23,[27][28][29][30][31]. Most of these algorithms are based on processing of 2D (bright-field) images of quasi-2D systems. Mohraz and Solomon, however, were one of the first to determine the 3D position and orientation of uniaxial ellipsoidal particles, i.e. all five coordinates, using 3D confocal microscopy and a novel anisotropic feature-finding algorithm [23]. Their algorithm identifies the points that are located on the central axis (or backbone) of a rod. These points are then grouped together by cluster analysis as individual rodbackbones, from which the centroid location and orientation are determined. This algorithm enabled the quantitative determination of the 3D translational and rotational motion of a dilute suspension of ellipsoids [32].
Quantitative 3D real-space study of concentrated phases of anisotropic particles is, however, much less progressed compared to studies on spherical colloids. Progress has been made for suspensions of ellipsoids, where nematic order was found using a centrifugal field [33] and local crystalline order with an external electric field [34]. In contrast with the system of ellipsoidal particles, it was recently shown by some of us that a system of fluorescent silica rod-like particles forms both nematic and smectic phases in equilibrium [21,22]. However, determination of all the 3D positions and orientations of the particles in these dense phases was not possible due to the significant overlap of fluorescent signals. In this paper we demonstrate a novel 3D image processing algorithm that is capable of quantifying fluorescent silica rods in concentrated (liquid-crystalline) phases. The algorithm is tailored to work even when the fluorescent sig-nals of the particles overlap considerably and a threshold method and subsequent clusters analysis alone does not suffice. The algorithm in principle also works for other uniaxial particles such as ellipsoids or dumbbells. This paper is organized as follows. First, we describe the basics of particle-locating algorithms. Second, we describe our algorithm in detail. Third, we demonstrate the performance of the algorithm with 3D image stacks of concentrated fluorescent silica rods. Then, we illustrate that our algorithm can also be applied to 3D electron tomography data of gold nanorods. Next, we evaluate the accuracy of the algorithm by measuring the translational and rotational motion of non-overlapping rod-like particles. Finally, we compare our results with recent progress in the field and give an outlook on further studies that the algorithm enables.

Locating particles in confocal microscopy data sets
The aim is to identify and locate (rod-shaped) particles in a set of real-space images (or snapshots) and to obtain the full configuration of the system. A specific configuration of a system of particles is given by a set of parameters, one for each degree of freedom of every particle. In the case of rods, these degrees of freedom for particle i are centre position r i , orientationû i and possibly length l i , diameter d i and brightness b i . If the length and diameter are known in advance they can be fixed, but if the particles vary in size they can also be left as free parameters. If the particles vary in brightness this can be added as an additional degree of freedom. Variations in brightness can be caused by the synthesis method, scattering or shading in the sample, but also by photo bleaching. In the case of fully symmetric, homogeneously dyed rods it is not possible to distinguish between the two ends of the rods. However, we also synthesised rods with a gradient in brightness, with one bright and one much darker end [35], of which the orientation could be fully determined. To keep the notation short we introduce p i = {û i , l i , d i , b i } which contains all the degrees of freedom except the position.
To obtain the configuration (r i and p i ) we need to elaborate on what is measured. In case of fluorescent confocal laser scanning microscopy we can assume that the imaging system is linear so that we can add intensities. The measured image intensity M (r) at position r can be written as the sum of the ideal (noiseless or averaged) images of the single particles, and RSP(r, p i ) is the image of a single particle placed in the origin, or rod spread function (RSP). The image of a single particle at the origin depends on all the internal degrees of freedom of the particle such as orientation, length, diameter and brightness, but also on the point spread function (PSF) of the imaging system. It is given by which is a convolution ( * ) of the dye distribution ρ dye (r, p i ) of particle i placed in the origin and the PSF. In a dilute sample this RSP(r, p i ) can be measured directly but it can also be calculated when the dye distribution is simple and the parameters of the optical systems are known. Different approaches to obtain the particle coordinates are possible. If all the parameters such as the PSF and RSP are known, the locating problem becomes in principle a deconvolution. However, the RSP and the PSF can be time consuming to determine accurately, and deconvolutions are sensitive to small changes in the kernel function [36]. This is unfortunate since e.g. polydispersity will introduce changes in the RSP which would make the deconvolution difficult. If the RSP is not known, there exist several other possible options. The first option is to assume that the overlap between the RSPs is not too severe and to determine centre-of-mass and orientation with methods that are insensitive to the details of the optical system. This is the method used by centroiding algorithms and is also the method used in this article.
Another option is to use a Bayesian method [37]. This method searches for the configuration that has the largest probability of having resulted in the observed image. This method has proven to work well for two-dimensional data sets [38]. It is, however, slow and complex and therefore not practical for large three-dimensional data sets.

Generation of test images
To test our algorithm we generated 8-bit confocal-like images from sets of computer-generated particle trajectories. Using the centres-of-mass r i and particle orientationsû i , we generated 3D stacks of xy-images of spherocylinders with aspect ratio l/d = 5, where l is the endto-end length of the particle and d the diameter. This was done by calculating the closest distances D to a line segment, representing the backbone of a particle. The distance from a point in the origin to a line segment from x 1 to x 2 with length l = |x 1 − x 2 | is given by where α = (û · x 1 ) andû = (x 1 − x 2 )/l the unit vector along the length of the line segment. If this distance was less than the diameter of the particle, the pixel was given a value of 0.95. This was then repeated for all particles. We approximated the effect of the PSF in our test images by convolving them with a Gaussian kernel with fixed standard deviation σ x /d = σ y /d = 0.3 and σ z /d = 0.3, 0.6 and 0.9 with d the diameter of the particle. The full-width-at-half-maximum (FWHM) of the Gaussian function, given by 2 √ 2 ln 2 σ i , is a direct measure of the resolution of the images. Besides variation of resolution, we also varied the amount of noise in the images. Although noise from modern detectors is essentially photon-limited, suggesting a Poisson distribution [39], we added noise to each pixel in our images with a simple Gaussian distribution with standard deviation σ n = 0.10 − 0.30. Because the amount of noise is known a priori, it is still straightforward to calculate the signal to noise ratio (SNR), which we define as SNR = (σ 2 g /σ 2 n − 1) 1/2 , with σ 2 g the variance of the constructed image and σ 2 n the variance of the noise [15]. Finally, we converted all our data, with pixel-values between 0 and 1, to 8-bit grayscale tiff images.

Our algorithm
To demonstrate the three-dimensional rod tracking algorithm we will first illustrate all steps of the algorithm with an artificially created set of images of a single rod, shown in Fig. 1. This will allow us to demonstrate clearly what is going on on a single pixel/voxel level. Later we will demonstrate how the algorithm fares with real colloidal suspensions. The following description is for three dimensions but most of the steps are straightforward to modify for two dimensions.
Reading In Figs. 1a-c we show three orthogonal slices through a generated 3D image that acts as the source image. The particle shown in the image has a diameter d = 13.0 pixels, and is blurred with a Gaussian kernel σ x /d = σ y /d = 0.3 and σ z /d = 0.9. Gaussian pixel noise of σ n = 0.1 was added to the image. The first step is to read in these source images. To avoid accumulating rounding errors and to allow the use of images of arbitrary bit depth we perform all image manipulations on floating point numbers between zero and one. Next, the image is rescaled to make sure the voxels are cubic, which is often not the case for confocal microscopy image stacks. The rescaling avoids having to account for different x, y and z scales in all following routines. To make sure no information is lost, this is done by enlarging the image using a bicubic interpolation. Care should be taken not to use overexposed images since this will result in a loss of information and an increase of positional error. See Ref. 15 for a more detailed description and the optimal shape of the intensity histogram. We generally choose the magnification such that the particles are approximately 10 pixels in diameter. Larger magnification results in a large file size without any additional benefit.
Filter The aim of the first filter step is to reduce image noise. We apply a Gaussian blur to the image, i.e. a convolution with a Gaussian kernel, that acts as a low pass filter. The optimal width of the function depends on the noise level in the images; a value between 1.5 and 3 pixels was found to give the best results for the images obtained in the present paper. A value that is too large will result in the loss of resolution and in missing particles, a value that is too small will result in additional, incorrectly identified, particles. To ensure a black background for the particles, a background value is subtracted from every pixel. This background value is assumed to be mostly the result of photon noise, but it can also originate from other sources such as fluorescence from the solvent or immersion fluid. Pixels that have a negative value after the background value has been subtracted, are set to zero. In most cases a background value between 0.01 and 0.1 is used. This value should be chosen such that approximately half the empty pixels (not containing a particle) of the image are zero. We also save a copy of the image that has not been filtered. This allows us to perform the final fitting step on the original image. An example of a computer-generated image that has been filtered is shown Fig. 1d.
Well separated particles When the intensity distributions of the individual particles do not overlap significantly we apply what we call a threshold method. This threshold method works as follows. A typical value for the threshold is between 0.4 and 0.7 and can be determined by plotting a histogram or by a quick test on a single image in a program like Photoshop, Gimp or ImageJ. The next step is to group all connected pixels above the threshold value into sets, as described in the next section. This methods works when these sets of pixels belong each to a single particle and each particle only corresponds to a single set of pixels. In Fig. 1e an example is shown of the threshold method applied to a single particle. All pixels above the threshold are marked in yellow. The particle coordinates can be obtained by applying a fit to these sets of pixels, as described later in this section. When this threshold method works, it is preferred over more complex methods since it is both robust and accurate.
General case When a threshold does not successfully separate the image into regions belonging to single particles another method has to be used. The first step of this method is similar to the Crocker and Grier algorithm and is aimed at providing the final fitting steps with a good initial starting point.
In this step, we roughly locate the line segment starting from one end of the rod and ending at the other end, called the backbone of the particle. To locate the backbone, we look at all voxels brighter than a predetermined cut-off value. A good value for this is in general between The different stages of identification of the position and orientation of a single rod-shaped particle. (a,b,c) Orthogonal slices through a computer generated 3D image stack. The particle has a diameter of 13.0 pixels and is blurred with a Gaussian kernel with width σx/d = σy/d = 0.3, and σz/d = 0.9. Pixel noise has been added by adding Gaussian noise with σn = 0.1. (d) The same image after the filter-step. (e) After a threshold step, the pixels above the threshold are marked in yellow. (f) After a backbone step, the pixels identified as backbone pixels are marked in yellow. (g,h,i) The rod as it is located, viewed from the xy, yz and xz plane. (j) The histogram of the average intensity along the rod length after smoothing and background removal. The dashed vertical lines mark the fitted end-points of the rod. 0.1 and 0.5 depending on the intensity fluctuations between the rods. For these bright pixels we then check whether they are part of a backbone. To do this we first note that all local maxima should be part of the backbone. To check if the brightness of the pixel is a local maximum we compare its intensity to that of all pixels within a distance r bb . If none of these pixels are brighter the pixel is a local maximum. To find the parts of the backbone that are not on a local maximum, we look at the distribution of brighter pixels around the pixel in question. If the pixel is part of the backbone they should be on a ridge. Backbone pixels can have brighter pixels to one side or two sides but all these brighter pixels should be more or less on a line through the pixel in question. So to check if the pixel is part of a backbone we need to check if the pixels brighter than the pixel in question are on a straight line. To do this we fit a line to these bright pixels and sum the squared residuals χ, the squared distance between the brighter pixels and the line. If these bright pixels are part of the backbone of a rod this number will be low since the pixels will form an almost perfect line while on other places they will not form a line and the residuals will be much higher. We found that r bb = 3 pixels and a maximum value χ max = 80 work well for all our data. This step depends on the initial filtering and on the thickness of the rod in pixels. Fig. 1f shows the pixels that have been identified as backbone pixels in yellow.
After having identified the backbone pixels, we group them into connected clusters. Due to noise there can be small gaps between the backbone pixels of a rod, so we use the same search range r bb as before to identify neighbouring pixels. This should work as long as the diameter of a rod is larger than r bb .
We now have groups of pixels most likely belonging to a single rod. To continue, we fit (least square) a straight line to these pixels using a singular value decomposition and the algorithm described in Ref. 40. The coordinates resulting from this fit are accurate, but still have a strong pixel bias since they only fit to a few backbone pixels. To eliminate this bias and to obtain more accurate results, we use these coordinates, lengths and orientations as a starting point to fit the real image again.
Fitting The fitting steps work best when applied to the unfiltered image. The Gaussian blur filter will result in an additional overlap of the RSPs which can result in a decreased accuracy. The fitting is done in three steps; first the centre of mass of each group of pixels is computed, then the orientation is fitted and finally the length is fitted. The position is taken from the centre of mass, weighted with the pixel intensity, of the pixels within half a diameter from the previous fit. The orientation is obtained by fitting a straight line to these pixels where the fit is weighted with the intensity of the pixels using the same least square fitting algorithm as for the backbones. The length is obtained by calculating the average intensity of pixels along the rod length, see Fig. 1j. The histogram that is obtained from this is smoothed with a Gaussian kernel to avoid noise. The end points are then obtained by determining where the histogram value drops below I end I max , where I max is the maximum intensity value in the smoothed histogram and I end is a parameter that can be set manually. Usually a value of I end = 0.6 -0.8 was found to give good results, see the (blue) dashed lines in Fig. 1j. To obtain sub-pixel accuracy we fit a straight line to the 2 pixels above and 2 pixels below the point where the histogram crosses this value. To determine which pixels to take into account in the generation of the histogram and the other fits, we use the pixels within one radius of the central line segment of the previous fit. Therefore, the result of the fit might improve when the step is repeated. The fitting algorithm normally converges in one or two steps. If this is not the case there is something wrong with the data or one of the parameters. Figs. 1g-i show the same orthogonal sections as Figs. 1a-c with the backbone of the rod highlighted in yellow and the outline of the rod (resulting from the fit) highlighted in magenta.
Filtering The final step is to filter out particles that are found more than once, particles that do not contain enough intensity or sometimes particles that are not long enough. Ideally not much filtering is required.

3D particle tracking
To study particle dynamics, we applied our algorithm to time-series of 3D image-stacks. We first identified the positions and orientations of the rods in each 3D stack separately. Then, we obtained the particle trajectories using standard IDL-based routines [13]. To uniquely track the tip of the (up-down indistinguishable) rods, it is required that the angular displacements between successive frames [û(t + 1) −û(t)] 2 < 2. Therefore, care was taken that displacements with [û(t + 1) −û(t)] 2 > 2 were negligible. We then calculated the mean squared displacement (MSD) and the mean squared angular displacement (MSAD). We fitted the MSD to the expression with D t the rotationally averaged translational diffusion coefficient and t the error in measurement of each of the coordinates of the particle [41]. For the MSAD we used the expression [42,43] with D r the average rotational diffusion coefficient and r the measurement error in the determination ofû(t). For short times, equation (5) reduces to To estimate the sedimentation velocity at infinite dilution, assuming complete decoupling of rotations, translations and sedimentation [44], we use the Svedberg equation [45] v with v p the volume of the particle, g the gravitational acceleration, ρ p the mass density of the particle and ρ s the mass density of the solvent.

Expressions for the diffusion coefficients
To test the validity of our experimental measurements of the diffusion coefficients, we compared them to analytical expressions for hard cylinders at infinite dilution, as proposed by Tirado, Martinez and de la Torre [46], with η the solvent viscosity, p = l/d the aspect ratio of the particle and δ i a correction term for the finite aspect ratio of the cylinders, given by [46] δ ⊥ = 0.839 + 0.185/p + 0.233/p 2 , Experimental methods

Dense sediments of silica rods
For the preparation of dense samples of silica rods, two different batches of particles were used. The first batch consisted of rods with length l = 2.37 µm (δ = 10%) and diameter d = 640 nm (δ = 7.5%), with δ the polydispersity (standard deviation over the mean) [21]. A transmission electron microscopy (TEM) image of these particles is shown in Fig. 2a. The particles contained a non-fluorescent core, a 30 nm fluorescein isothiocyanate (FITC) labelled shell, and a 190 nm non-fluorescent outer shell. For the second batch of silica rods, with length l = 2.6 µm (8.5%) and diameter d = 630 nm (6.3%), rhodamine isothiocyanate (RITC) dye was added during synthesis, which resulted in an intensity gradient of dye molecules along the major axis of the particle [35]. The particles were coated with a 175 nm non-fluorescent outer shell. Particle suspensions were prepared by dispersing the rods in an index-matching mixture (n 21 D = 1.45) of either dimethylsulfoxide (DMSO) and ultrapure water (Millipore system) or glycerol and ultrapure water. The particles were first dispersed in DMSO or glycerol, after which water was added until the suspension was indexmatched by eye. This resulted in mixtures of 91 wt% DMSO in water and 85 wt% glycerol in water.
Next, sample cells were constructed with standard microscopy slides and No. 1.0-1.5 glass coverslips (Menzel-Gläzer). After the cells were filled with the suspension, they were sealed with UV-glue (Norland No. 68). The suspensions were imaged with a confocal microscope (Leica SP2 or Leica SP8) using a 63x/1.4 or 100x/1.4 oilimmersion confocal objective (Leica). We corrected the 3D images for distortion of the axial (z ) distances due to the refractive index mismatch between sample (n 21 D = 1.45) and immersion oil (n 21 D = 1.51), which resulted in an increase of axial distances of 5% [47]. Fig. 2b shows a 3D confocal microscopy image of a single rod suspended in 85 wt% glycerol in water. In Figs. 2c-e, three orthogonal slices through this 3D volume are shown. The larger width of the PSF in the axial (z) direction is clearly visible. Notice that the pixel size in x,y (50 nm) is smaller than in z (78 nm). Figs. 2f-h show the same rod after rescaling to cubic pixels, filtering and particle-fitting. In Fig. 2i, we show the intensity histograms of two rods that were oriented parallel to the xy image plane of the confocal microscope. The continuous (red) line shows the intensity histogram of a single uniformly dyed rod and the dashed (green) line that of a gradient-dyed rod [35].

Freely diffusing silica rods
For the experimental measurements on a dilute suspension of silica rods we used particles with length l = 3.3 µm (δ = 10%) and diameter d = 550 nm (δ = 11%), as measured with TEM. The particles were fluorescently labelled with a 30 nm (FITC) shell. The particles were dispersed in an index matching mixture of 85 wt% glycerol in water. The density of the solvent mixture was ρ = 1.222 g/ml [48] and the viscosity η = 92 cP (22 • C), as measured with an SV10 viscometer (A&D Company). This mixture not only matches the refractive index of the particles (n 25 D = 1.45), the high viscosity slows down the particle dynamics enough to measure their short-time self-diffusion in 3D. Because the density of this mixture is significantly lower than the density of the particles ρ = 1.9 g/ml [35], sedimentation cannot be avoided. We assume, however, complete decoupling between translational motion, rotational motion and sedimentation [44]. A fused quartz capillary (Vitrocom) was filled with a dilute suspension (volume fraction φ < 1%) of the fluorescent silica rods. The suspension was imaged with a confocal microscope (Leica SP8) equipped with a fast 12 kHz resonant scanner and hybrid detector. Images with 8-bit pixel-depth were acquired using a white light laser with a selected wavelength of 488 nm. A confocal glycerol immersion objective 63x/1.3 (Leica) was used, which is optimized for refractive index n D = 1.45. If we assume a Poisson distribution of the noise, we can easily estimate the signal to noise ratio (SNR) of a single image because of the photon counting mode of the hybrid detector. We use the definition SNR = √ n p with n p the number of detected photons in the brightest part of the image [49]. To avoid hydrodynamic interactions with the wall, particles The scale bar is 800 nm. (f-h) The particle after rescaling, filtering and fitting. The magenta outline indicates the final fit from which the position and orientation is computed. (i) The (normalized) intensity histograms along the major axis of two differently dyed particles that are oriented in the xy plane, obtained from confocal microscopy images. The (red) solid line is from a uniformly dyed silica rod. The (green) dashed line is from a silica rod with a gradient in dye distribution.
were imaged 20 µm deep into the sample. We recorded 800 repeats of 3D image stacks consisting of 512 x 261 x 66 pixels with voxel size 144 x 144 x 331 nm. The time to record a single 3D volume was τ = 1.80 s. During this time, the particles are expected to translate on average √ 2 D t τ = 110 nm in each direction and rotate only √ 4 D r τ = 0.1 rad.

AuNRs@SiO2 & 3D electron tomography
For the fabrication of a spherical cluster of nanorods, we first synthesized gold nanorods following the method described in Ref. 50. Next, the gold rods were coated with a layer of mesoporous silica (AuNRs@SiO 2 ) [51], which resulted in particles with length l = 119 nm and diameter d = 68 nm, as measured with TEM. Afterwards, clusters were fabricated via an emulsification process [52,53]. Brightfield TEM tilt series of an 11-particle NR-cluster were acquired by tilting the sample over a range of -65 • to 65 • and recording images every 2 • . Images were taken on a Tecnai 20 (FEI) transmission electron microscope, operating at 200 kV with an LaB 6 electron source, in bright field mode. Tomographic reconstructions of the images were made with the iMOD software package using the simultaneous iterative reconstruction technique (SIRT) [54,55]. After reconstruction, the data stack was filtered using a low frequency Fourier filter (iMOD) and inverted to ensure light particles on a dark background to enable individual particle identification.

Determination of 3D particle positions & orientations in dense suspensions
To test our 3D particle-fitting algorithm we identified the fluorescent particles in a concentrated suspension of silica rods, as shown in Fig. 3. The particles were uniformly dyed, had a length l = 2.37 µm (10%), diameter d = 640 nm (7.5 %) aspect ratio l/d = 3.7 and were dispersed in a 85 wt% glycerol in water mixture. Small regions of hexagonally stacked particles existed in the sample (Fig. 3c), however there was no long-ranged order in the sample and particles seemed jammed or arrested in different orientations. Fig. 3a shows that 5.2 µm deep in the sample, the fluorescent signals of the particles did not overlap significantly in xy, despite the high particle concentration. This is due to the 190 nm non-fluorescent outer shell of the particles which was deliberately grown around the particles during synthesis to resolve them individually, even when they were lying side-by-side. However, the orthogonal slices in Figs. 3b-c show that particle signals did overlap in the z -direction, even after noise filtering. Nevertheless, by visual inspection of the (magenta) particle outlines in Figs. 3d-f we conclude that the algorithm correctly identified the orientations and positions of the particles, despite the high particle concentration. Fig. 3g shows a computer generated reconstruction of the sample, with colours indicating the 3D orientation of the particles. Fig. 4 shows a second example of the performance of our fitting-algorithm in a concentrated suspension. The rods in this sample had length l = 2.6 µm (8.5 %), di- ameter d = 630 nm (6.3 %) and were dispersed in an index-matching mixture of DMSO/water. After the particles had been left to sediment for several days, they ordered into smectic layers, more or less parallel to the xy-plane, as can be seen from Fig. 4a (12.5 µm deep in the sample). It can also be seen that the particles had an intensity gradient along their major axis and that there was significant overlap of the fluorescent signals in the xy-image (Fig. 4a). As expected, it was even more difficult to resolve individual particles in the z-direction (Figs. 4b-c), however, it is clear from the hexagonal pattern in Fig. 4b that the particles formed a smectic-B phase. The magenta outlines in Fig. 4d-f show the result of the particle fitting. By visual inspection of the outlines in the complete image-stack (containing 1699 particles), we conclude that > 98% of the particles had been correctly identified by the algorithm. In Fig. 4g we show a 3D reconstruction of a part of the image-stack, which clearly shows 3D orientational order, smectic layering, and transverse (red and blue) particles.

Determination of 3D positions and orientations of gold nanorods
Although our algorithm was written for analysis of confocal microscopy images, it is also applicable to other 3D image-stacks of uniaxial symmetric particles.
As an example, we show results of the identification of gold nanorods (AuNRs) from a 3D transmission electron microscopy (TEM) tomographic reconstruction in Fig. 5. The TEM micrograph in Fig. 5a shows the gold nanorods (in black), that were coated with a layer of mesoporous silica (dark grey). Fig. 5b-c  orthogonal sections through the 3D reconstruction of the cluster. The images were inverted to enable particle identification with our algorithm. Fig. 5d-e show the same orthogonal sections after filtering, with identified particles outlined in red. Finally, Fig. 5f shows the 3D reconstruction, with color-coding of the 3D orientation of the rods. The algorithm had identified all 11 particles and the reconstruction clearly shows that there was some degree of orientational ordering inside the cluster.
We are aware that a substantial amount of information on the 3D structure of the nanoparticles can be measured directly (and manually) from the 3D tomogram itself. Our image-processing algorithm however can determine unambiguously the 3D positions and orientations of the particles and can therefore be useful for the quantification of (larger) nanoparticle assemblies and should in principle also work on other types of samples, e.g. self-assembled clusters of nano-dumbbells [56].
Testing the accuracy of the algorithm for non-overlapping particle signals In this section, we assess the accuracy of our algorithm in more detail. We focus on the fitting accuracy of the algorithm when applied to images containing particle signals that are well separated. Although this situation is much less demanding compared to partially overlapping signals, care has to be taken when fitting this type of data as well. The main reason is that the (fluorescent) diameter of typical rod-like particles used in our experiments (d f l ∼ 300 nm) is comparable to the resolution of a typical confocal microscope (200−300 nm in the lateral and 500 − 700 nm in the axial direction [47,57]). Additionally, the PSF itself is anisotropic, which can result in a (strongly) distorted particle shape. Things become progressively worse when there is a refractive-index mismatch between the sample and immersion fluid, which deteriorates the PSF, introduces an intensity fall-off with height and distorts axial distances [47,58].
We therefore determined the accuracy of our algorithm using two approaches. In the first approach we investigated both the effect of the theoretically approximated PSF and the effect of noise on particle tracking accuracy using computer generated data. The second approach consisted of an experimental measurement of the translational and rotational diffusion of a dilute suspension of silica rods.
Details on the construction of the test-images and variation of the theoretically approximated PSF and noise can be found in the Supplementary Information. The results are summarized in Table I, which shows that for our worst-case scenario of a z-resolution of 636 nm and signal to noise ratio of 1.7, we obtain for the error in the determination of the main-axis of the rod r = 0.07 rad, which corresponds to a small measurement error of 4.1 • , see also Fig. S1. Additionally, we did not find any significant pixel-bias in either the position or the orientation (see Fig. S2). For the error in the positional measurement, we found t /d ∼ 0.05, which indicates sub-pixel accuracy. Finally, we measured the diffusive motion in a dilute suspension of fluorescent silica rods experimentally, which provides a real-life test of the accuracy of our algorithm. The rods that were used had length l = 3.3 µm (δ = 10%), diameter d = 550 nm (δ = 11%) and aspect ratio l/d = 6.0. From the number of photons in the brightest part of the image we estimated the signal to noise ratio to be SNR ≈ 3, which is in the range stated in Table I. The tracking results, averaged over 8 particles, are shown in Fig. 6. A typical translational trajectory of 12 min is shown in Fig. 6a. From a fit to the average linear displacements (of all 8 particles) in the z-direction, we estimate the sedimentation speed to be v sed = 0.331 ± 0.005 µm/min. This value is slightly higher but comparable to the value of v sed = 0.28 µm/min that we obtained from equation (7). For further analysis we subtracted the average linear displacements from the trajectories. Fig. 6b shows a rotational trajectory of 12 min for a single particle. In Fig. 6c we show the probability distribution of the norm of the displacement |∆r|, for three different timesteps ∆t. In Fig. 6d we show the same distribution for the norm of the displacements of the unit vector |∆û|. The solid black lines in Figs. 6c,d are fits proportional to |α| 2 exp(−|α| 2 ) with α = ∆r, ∆û respectively.
To extract the translational diffusion coefficient, we calculated the rotationally averaged mean squared displacement ∆r 2 , as can be seen in Fig. 6e. For ∆t > 10 s we found that ∆r 2 ∼ t 0.97 indicating diffusive behaviour. The statistical error in the individual measurement points is smaller than the symbol size. Fitting the data with equation (4), we obtain the shorttime rotationally averaged translation diffusion coefficient D t = (3.06 ± 0.01) × 10 −3 µm 2 /s and static error t = 45, 46 and 59 nm in the x, y and z-direction respectively, which confirms that we can locate the particles with sub-pixel accuracy. The value for D t is in strong agreement with the theoretical value obtained from equation (10) which is D t = 3.2 × 10 −3 µm 2 /s. Finally, we calculated the mean squared angular displacement ∆û 2 , as shown in Fig. 6f. This time we obtained ∆û 2 ∼ t 0.92 and for the short-time rotational diffusion coefficient D r = (1.32 ± 0.02) × 10 −3 rad 2 /s. This is in good agreement with the theoretical value D r = 1.5 × 10 −3 rad 2 /s, obtained from equation (11). For the corresponding rotational relaxation time we found τ r = 1/(2 D r ) = 3.8 × 10 2 s, which confirms that we measured in the short-time diffusion regime. From the fit we also obtained r = 0.07 rad, which corresponds to a small angular uncertainty of ∼ 4 • .

DISCUSSION
In this paper we demonstrated a new image-processing algorithm that is capable of extracting the positions and orientations of fluorescent rod-like particles in both dilute and concentrated suspensions. Although the algorithm was written for three dimensions, most steps are straightforward to modify for two dimensions [59]. Mohraz and Solomon [23] were the first, as far as we know, to describe an algorithm that can detect the position and orientation of ellipsoidal particles in 3D confocal microscopy images and this work follows a similar approach. The algorithm of Mohraz and Solomon groups clusters of pixels together to form backbones but does not use additional fitting steps, which we found necessary to correctly iden- tify particles when there is significant overlap of particle signals. The difference in particle geometry (ellipsoids versus rods) combined with the small (fluorescent) particle diameter in our study might be the reason why we find that using only a maximum threshold and cluster analysis is not sufficient to identify rods in concentrated suspensions, even when the rods have a large (> 150 nm) nonfluorescent shell and a considerable electric double layer (∼ 100 nm) [21,22]. The rod-like particles used in this study have a repulsive interaction potential and therefore form dense smectic-like phases, which we now can identify on the single-particle level in the bulk. The algorithm also enables the study of glassy phases of anisotropic particles in three dimensions, which is promising since all current real-space glass-transition studies of anisotropic particles so far are either 2D [60][61][62] or tracer-host [63]. Finally, we would like to mention that the algorithm is also applicable to study the dynamics of (concentrated) 'active colloids' (e.g. self-propelled particles and bacteria), a field that is rapidly emerging [64].
Since the typical fluorescent diameter of the rod-like particles is around 300 nm, deconvolution of the imagestacks before particle fitting can be useful when particles are difficult to resolve individually. The necessary higher (Nyquist) sampling rate, however, is not always practical or even not possible for faster moving particles. Additionally, deconvolutions are sensitive to small changes in the rod-spread-function (RSP), introduced by e.g. polydispersity. A clear improvement of the algorithm, therefore, is to fit the particles with the RSP, analogous to the fitting of the sphere-spread-function (SSF) reported by Jenkins et al. [15], which is work currently ongoing. With this type of extension of the current algorithm it should also be possible to accurately measure in-situ particle polydispersity, which is known to have a large effect on e.g. the liquid-crystalline phase behaviour [65].
By measurement of freely diffusing rods, we acquired additional information on the accuracy of our algorithm. Although the motion is analysed in the lab-frame and therefore translational and rotational motion should be coupled [28], we did not observe such behaviour since the friction anisotropy in our 3D measurement is small D /D ⊥ = 1.3 and because we averaged over an ensemble of particles and over many initial orientations. We found that the error in locating the rods ( t = 45, 46 and 59 nm in the x, y and z-direction respectively) confirms sub-pixel accuracy and agrees roughly with the criterion for spherical particles that t ∼ M/N , with M the pixelsize and N the diameter of the particle in pixels [13]. The value for the (short-time) rotational diffusion coefficient D r = (1.32 ± 0.02) × 10 −3 rad 2 /s, is one order of magnitude larger than previously accessible with 3D confocal microscopy [32] which is, however, due to the equipment rather than the image-processing. The error in the determination of the orientation of the rod ( r = 0.07 rad) is in the range of the values obtained via simulated test-images, shown in Table I. The rule-of-thumb that r ∼ 1/P a with P a the half-length of the rods in pixels [32] seems to hold quite well in our case, since 1/P a = 0.08 in our measurements. CONCLUSION We developed an algorithm that extracts the positions and orientations of rod-like particles from 3D confocal microscopy images. The algorithm is tailored to a system of fluorescently labelled silica rods and can identify these particles even in the bulk of 3D concentrated phases where the fluorescent signals of the particles overlap considerably. This allowed us to determine the 3D positions and orientations of particles in a concentrated disordered phase and in a liquid-crystalline smectic-B phase. The algorithm also works on electron tomography reconstructions of gold nanorods, which enables the 3D quantification of (large) nano-particle assemblies. It is also expected to work on other uniaxial particles such as ellipsoids or dumbbells. We determined the accuracy of the algorithm for varying z-resolution and noise levels from generated 3D test-images. Despite the (anisotropic) distortion of the theoretically approximated point spread function (PSF) and the low signal to noise ratio (SNR), the error in the determination of the orientation of the particles remained small. These results confirmed that we can accurately track rod-like particles with (fluorescent) diameters down to 300 nm. With our algorithm and a fast confocal microscope we determined the translational and rotational motion of a dilute suspension of sedimenting silica rods. We demonstrated that the measured diffusive motion was in good agreement with theory (neglecting sedimentation) and that we can track the particles with sub-pixel resolution.
This novel algorithm therefore allows for studies of structure and dynamics on the particle level of dense liquid-crystalline phase behaviour (such as nematic, smectic and crystalline phases), but also allows for studies of the glass transition of anisotropic rod-like particles in three dimensions. Of course, the algorithm will also be applicable to dilute suspensions or in cases where rodlike particles are used as tracers, such as in biophysical or micro-rheology studies.