## Abstract

Volumetric imaging and 3D particle tracking are becoming increasingly common and have a variety of microscopy applications including *in situ* fluorescent imaging, *in-vitro* single-molecule characterization, and analysis of colloidal systems. While recent interest has generated discussion of optimal schemes for localizing diffraction-limited fluorescent puncta, there have been relatively few published routines for tracking particles imaged with bright-field illumination. To address this, we outline a simple, look-up-table based 3D tracking strategy, which can be adapted to most commercially available wide-field microscopes, and present two image processing algorithms that together yield high-precision localization and return estimates of statistical accuracy. Under bright-field illumination, a particle’s depth can be determined based on the size and shape of its diffractive pattern due to Mie scattering. Contrary to typical “super-resolution” fluorescence tracking routines, which typically fit a diffraction-limited spot to a model point-spread-function, the lateral (*XY*) tracking routine relies on symmetry to locate a particle without prior knowledge of the form of the particle. At low noise levels (signal:noise > 1000), the symmetry routine estimates particle positions with accuracy better than 0.01 pixel. Depth localization is accomplished by matching images of particles to those in a pre-recorded look-up-table. The routine presented here optimally interpolates between LUT entries with better than 0.05 step accuracy. Both routines are tolerant of high levels of image noise, yielding sub-pixel/step accuracy with signal-to-noise ratios as small as 1, and, by design, return confidence intervals indicating the expected accuracy of each calculated position. The included implementations operate extremely quickly and are amenable to real-time analysis at frame rates exceeding several hundred frames per second.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Recently, intense effort has been devoted to developing 3D localization microscopy techniques for applications that include *in situ* single-molecule fluorescence imaging, characterization of colloidal systems, and *in-vitro* single-molecule assays such as tethered-molecule force spectroscopy [1–5]. While scanning techniques such as confocal and light-sheet microscopy suffice for imaging stationary or slow-moving systems (such as fixed tissue samples), highly dynamic systems require frame rates only accessible with wide-field microscopy. In particular, force spectroscopy applications such as Optical Traps (OT), Magnetic Tweezers (MT), or Acoustic Force Spectroscopy (AFS) often require sampling rates as high as 1 kHz, so that high-speed wide-field or non-image-based tracking techniques are the only viable options [6–9]. Moreover, studies employing those techniques typically require nanometer-scale localization accuracy. One of the biggest barriers to setting up an assay utilizing 3D localization, particularly those incorporating force spectroscopy, is that there are few turn-key solutions available commercially. Consequently, new assays often require building or retrofitting a custom microscope. Moreover, while excitement over super-resolution fluorescence microscopy has driven a surge in manuscripts detailing optimal routines for tracking diffraction-limited spots, there have been relatively few algorithms published discussing how to track particles imaged under bright-field illumination.

This manuscript presents a set of straightforward, computationally efficient image analysis routines for 3D tracking of bright-field illuminated particles in real-time, with nanometer-scale precision. The scheme can be conceptualized in two parts. First particles are localized in 2D (within the image plane). Next, their appearance is used to determine their depth along the optical axis of the microscope. Unlike typical fluorescence tracking routines, which require prior knowledge about the image of a particle (e.g. the point-spread-function of the optical system), the routine presented here relies on image symmetry to determine lateral particle position. This is accomplished by extending the method presented by Parthasarathy [10]; in particular, we show how one can determine the uncertainty of particle position based on the noise characteristics of each processed image. Furthermore, this allows us to determine the optimal way to weight the image data for maximum accuracy. By fitting captured images to a measured, depth-dependent look-up-table, particle position along the optical axis can be computed quickly without requiring detailed knowledge of the characteristics of the condenser optics. We extend the look-up-table method, initially proposed in [11], Gosse, *et al*., to dramatically improve the resolution of axial positions. We additionally provide the uncertainty of the axial position for each measurement. Our methods build on the computational efficiency of [10–12] and thus can execute in real time at rates up to 1 kHz. Likewise our methods retain the hardware advantages of [11,12], Gosse *et al*. and van Loenhout, *et al*., in that they are applicable to bright-field microscopes modified only to include automated focus control.

#### 1.1. 3D imaging via diffractive scattering

As opposed to 3D scanning modalities (e.g. light-sheet, confocal), which record and distinguish voxels, the approach used herein discriminates axial position based on the appearance of a nano- or micro-particle when imaged under bright-field illumination. Mie scattering creates concentric interference fringes surrounding a particle. The periodicity and extent of these fringes depend on the wave-front of the incident light from the condenser, the particle size, and its distance with respect to the focal plane of the microscope. With the first two parameters held constant, the fringe pattern of a particle changes smoothly as it moves axially (*Z*). Consequently, the *Z*-position can be determined either through analytical fitting using digital holographic microscopy or by matching fringe patterns to a look-up Table [11,13,14].

While analytical fitting approaches do not require look-up-tables and therefore can be implemented without a prior calibration procedure, they are less versatile. For instance, diffraction patterns are objective-dependent and vary with optical alignments. Consequently, careful alignment of the objective and condenser in bright-field illumination is essential [13]. Furthermore, digital holographic routines are computationally intensive, and are ill-suited to situations where real-time readout is beneficial (particularly true for force-spectroscopy techniques). In contrast, look-up table (LUT) methods are computationally inexpensive and they do not rely on exacting microscope alignments provided the LUTs are recorded with the same settings.

The only major requirement in retrofitting a stock microscope for bright-field LUT-based depth tracking is the need for automated focus control to conveniently record LUTs. The LUT directly dictates the accuracy of depth determination. Various commercially available, motor-based, focusing systems allow coarse (micron-scale) *Z*-tracking, while nanometer accuracy can be obtained using piezoelectric translators. There are a number of high-precision piezoelectric objective scanners and microscope stage inserts that provide 10-20 nanometer accuracy and closed-loop repeatability. The data presented here was acquired using a Physik-Instrumente P-721.CDQ objective scanner and E-665 piezo controller (Physik Instrumente, Auburn, MA). Other scanners [6,15] have been used in similar implementations.

#### 1.2. Established tracking routines

Realizations of LUT-tracking vary, although, typical implementations adopt the algorithm described by Croquette *et al.* in their seminal publication detailing a magnetic tweezer [11]. As described, prior to the start of an experiment, a piezoelectrically-controlled objective scanner is used to vary the focus and capture a series of images of stationary particles similar in size and shape to the particles that will be tracked. Then, during the course of an experiment with a stationary objective, images of moving particles with surrounding fringe patterns are matched to the entry in the *Z*-stack that they most closely resemble, yielding the approximate focal depth. See Figs. 1(a) and (b) for illustration.

Before fringe patterns can be matched to the LUT to track lateral (*XY*) motion, centroids of particles must be identified. In fluorescence and dark-field microscopy, sub-pixel *XY* localization is typically accomplished by calculating the center of mass (intensity; COM, a.k.a. first-order central image moment) of bright puncta or fitting the puncta with an appropriate point-spread-function, via maximum-likelihood estimation (MLE) [16]. However, for diffractive scattering, neither approach is appropriate. For bright-field images where diffractive fringes oscillate from dark to light on a gray background, COM (which weights pixels by intensity) does not make sense. The intensity and size of the fringe pattern can be predicted using Mie theory; therefore, fitting a particle’s image using MLE amounts to performing analytical holographic fitting (which can be slow and necessitates precise optical alignment, as discussed above). An alternative to COM and MLE is to localize particles based by considering their origin of symmetry within an image. The origin of symmetry is often determined using Discrete Fourier Transform (DFT)-based cross-correlation [6,12,17,18] which involves discriminating the sub-pixel maximum of an image convolved with its mirror image (see appendix for a complete mathematical description). For large images this can be time consuming. *M* × *N* pixel images require approximately 2*MN* log_{2} *MN* operations. Computation time can be reduced by employing a 1D convolution scheme. If the axes of symmetry of a particle are already known to within a few pixels (via COM or peak detection, for example), convolution and parabolic fitting routines can be applied to extract the *x-* and *y*-coordinates from the XY bands centered on the particle. Unfortunately, noise can significantly impact the 1-D approach, because off-axes information is discarded. A hybrid approach presented by van Loenhout *et al.* aims to mitigate this loss using a quadrant interpolation (QI) scheme to construct 1D intensity profiles that roughly correspond to the azimuthal averages of each quadrant, which are then cross-correlated to determine the approximate error in the initial COM-estimate [12]. Nevertheless, in all three cases the convolution routines suffer from two major shortcomings: First, subpixel localization cannot be determined without either some form of interpolation or image-upscaling and second, convolution routines do not provide an obvious mechanism for determining localization confidence limits.

#### 1.3. Features of our 3D LUT tracking algorithm

Force spectroscopy assays routinely require localization resolution on the order of 10-20 nm. With a 100× objective and pixel sizes ranging from 16 µm, for sensitive, scientific-grade cameras, down to 4 µm, for high-speed, industrial vision cameras, 10 nm resolution translates to localizing a particle to within 0.06-0.25 pixels. Consequently, routines for maxima detection or kernel convolution that have no better than pixel accuracy are inadequate. Similar requirements apply to axial-localization. For example, the objective scanner used for the experiments described below has a minimum step size of 20 nm, meaning discrete matching is insufficient for achieving <20 nm accuracy. Also, due to vibration, mechanical drift, and imaging noise, it is virtually impossible to record LUT steps of exactly equal intervals, see Fig. 1(c) (inset). Thus, to overcome the limitation of discrete sampling and measurement noise, some form of interpolation is needed. We achieve this by adapting the LUT algorithm to use splines for fast interpolation with accuracy limited only by the level of image noise.

Additionally, particle tracking analyses can be improved by including precision estimates, which enables systematic identification of low quality or missing data and the propagation of localization error through subsequent calculations. For example, characterization of both optical trap stiffness and force applied via Magnetic Tweezing involve tracking the Brownian motion of a confined particle [8,19]. For stiff traps (or high tension) particles become very confined, in which case localization errors greater than tens of nanometers (several tenths of a pixel, as determined by the resolution estimates above) can result in wildly inaccurate measurements. As such, localization precision estimates can be used to determine the limits of stiffness/force sensitivity. In this manuscript we demonstrate how to estimate precision for both *XY* and *Z* localization algorithms.

Finally, for high frame rate applications (such as that required for the fluctuation analysis used in force-spectroscopy), capturing and storing images to disk for offline processing is problematic due to both the data storage requirements and the delays which prohibit real-time decision-making during experiments. Thus, fast (100-1000 frame/sec or better) and parallel-processing algorithms are advantageous. Building upon pre-existing, efficient algorithms, we present a high-speed, 3D particle tracking routine which achieves excellent *XY* resolution even in the presence of noise, sub-step *Z* precision, and numerical precision estimates.

## 2. Methods

#### 2.1. Radial symmetry detection

We start by describing our routine for determining a particle’s lateral position. As an alternative to image convolution, the image of a particle possesses an origin of radial symmetry that can be determined with on the order of 3*MN* operations using radial symmetry detection [10], the details of which are reprised in Eqs. (1)–(5), below. The radial symmetry routine leverages the fact that the gradient vectors of a radially symmetric image will always be directed along lines through the origin, illustrated in Fig. 2(a). A linear least-squares fit of the point at which all of the gradient lines in an image intersect yields the origin of symmetry; we extend the radial symmetry routine by showing that the mathematical-machinery of least-squares fit can be used to calculate (per-image) error estimates.

Without noise, the origin of a radially symmetric image lies at the intersection of all gradient-lines. For a noisy image, the origin (*x _{0}, y_{0}*) is the point which minimizes the distance to each line. Consequently, the origin of symmetry can be determined by linear least-squares where the error function being minimized is given by

*x*) specifies the coordinate of the

_{k}, y_{k}*k*-th pixel, and $({\Delta}_{x}{I}_{k},{\Delta}_{y}{I}_{k})=\stackrel{\rightharpoonup}{\nabla}{I}_{k}$are the components of the gradient at each pixel. The gradient is numerically-computed via Eqs. (28) and (29). In matrix-form this amounts to computingwhere$x={[{x}_{0},{y}_{0}]}^{T}$, and

Image noise has an overwhelming influence on gradient-direction for pixels with low-magnitude gradients, i.e. intensities similar to their neighbors, see cyan “x” in Fig. 2(b). Consequently, noise biases the solution to Eq. (1) towards the center of the image.

The poor susceptibility to noise of the unweighted algorithm is due to two factors. First, pixels far away from the true symmetric origin collectively contribute greatly to the error function given in Eq. (1). Provided that the particle center can be roughly estimated by COM via Eq. (21), for example, then an exclusion filter can be applied [10]. Second, because the effect of image noise is most pronounced in areas where the gradient is shallow, those pixels inordinately bias the calculation. A solution is to also weight the least-squares fit as a function of the magnitude of the image gradient, using weighting matrix

**W**is the diagonal matrix defined in Eq. (4). In the original implementation, Parthasarthy somewhat arbitrarily chose${W}_{ii}=\frac{{\left|\nabla I\right|}^{2}}{{r}_{i}}$; however, as discussed in the next section, this choice can be further optimized.

#### 2.2. Localization error estimation and optimal weighting function

Building on the method outlined in Eqs. (1)–(5), we can extend the routine to also return (per-measurement) localization error estimates derived from the image noise. Additionally, by statistically characterizing error we can determine the optimal weighting function, Eq. (4).

For an un-weighted two-variable least-squares fit of $N$data points, the average residual of the estimate, $x$in Eq. (2), is given by

and the standard error of $x$is given by the matrixIn this context, the residual, defined as $R=(b-Ax)\text{,}$corresponds to**defined in Eq. (1). For the weighted least-squares, the average residual becomes**

*d***SE**describe the orientation and scale of the major and minor axis of an ellipse which characterizes the fit precision. The larger of the two eigenvalues can be interpreted as the radius of a circle representing the “worst-case” precision estimate for both

*x*and

*y*directions. In subsequent analysis and plots when refereeing to a

*scalar*“weighted standard error” we mean the leading eigenvalue of

**SE**.

A good choice for the weighting matrix, **W** in Eq. (4), can be determined analytically and further optimized by numerical evaluation. If the entire image (or region of interest) is assumed to have gradient vectors, which (in the absence of noise) point to the symmetric origin, then the propagation of error due to Gaussian-distributed image noise yields that the least-squares minimizer, ${d}_{k}$in Eq. (1), is expected to be distributed as

*m*and gradient-exponent

*n*can be determined by finding the values which minimize the weighted standard error, Eq. (9). For a given image morphology, the optimal values can be calculated by applying the algorithm to several representative images and evaluating which exponent values yield the lowest weighted standard error.

We generated patterns representative of bright-field imaged microspheres by first imaging a 2.8 µm diameter microsphere adsorbed to a glass slide using a 63×, 1.40 NA objective, brightly illuminated to maximize signal to noise. Next, the symmetry algorithm with weight exponent *n* = 2 in Eq. (4) and distance-dependence exponent *m* = 0 was used to determine the approximate center of symmetry of the image. Using the estimated center, $\tilde{x}$and$\tilde{y}$, each pixel was assigned a radial coordinate${r}_{ij}=\sqrt{{({x}_{i}-\tilde{x})}^{2}+{({y}_{j}+\tilde{y})}^{2}}$. Finally, pixel intensities versus radial coordinate were fit with a smoothing spline (using MATLAB’s csaps function), yielding a continuous diffraction profile *I*(*r*). Using *I*(*r*), images of arbitrarily located particles could be quickly generated by simply evaluating the function over a grid defined by the radial coordinate *r _{ij}*. Image noise was generated by adding a Gaussian-distributed random array$\eta [i,j]$:

*m*, and gradient-exponents,

*n*, for varying levels of noise. See Figs. 2(c) and (d) and Fig. 3. The routine performs best with a choice of

*n*≈5 and

*m*ranging from 0 – 1.5. As shown in Fig. 3, with

*n*= 5, the distance-exponent has little effect on accuracy, for

*m*<1.5. A map of the effective pixel weight with ${W}_{ii}=|\nabla {I}_{i}{|}^{5}$for an example image is shown in Fig. 2(e). Choosing

*m*= 0 simplifies the algorithm by eliminating the need to compute an initial centroid guess; although, if fringe patterns of neighboring particles “spill-over” into the region of interest then an initial guess can be useful for applying a size exclusion filter prior calculating the symmetric origin.

To summarize our method of *XY* localization, we use Eqs. (2), (3) and (11), with *m* = 0 and *n* = 5, to locate particles, and Eqs. (8) and (9) to determine the uncertainty of the position.

#### 2.3. Z-Localization through pattern lookup

Axial-tracking is accomplished by matching images of particles with the best-fitting entry in a previously recorded lookup table. For bright-field imaging, a common approach used by the MT and AFS community is to use a piezoelectrically actuated objective scanner to record a series of regularly spaced images of a stationary particle in order to assemble a lookup table of diffraction patterns of particles as a function of relative distance from the objective. During measurements, the LUT is used to match the measured radial profile $\widehat{I}[{r}_{i}]$with the most similar entry in the LUT:

In the simplest implementation, the discrete step ${z}_{j}^{*}$, is accepted as the position, limiting the depth-resolution to the sampling interval. To achieve sub-step resolution, the approach initially presented by Croquette *et al.*, and adopted by many others, is to fit the residual in the neighborhood of the j-th step, ${R}^{2}[{z}_{j}^{*}]$, to a quadratic function and use the vertex as an estimator of the sub-step position [11,15,17,18]. Unfortunately, this scheme has two major drawbacks.

First, it does not account for image noise and positional error during LUT capture. Due to vibration, thermal drift, and limited precision in focusing hardware, it is very unlikely that a LUT can be sampled at consistent discrete intervals. Moreover, camera noise means the captured intensity is also expected to deviate from the “true” value. Consequently, each step of a LUT should be expected to fluctuate around its ideal value

This fluctuation is illustrated in Fig. 1(c) (inset). The example LUT was constructed by recording the azimuthally averaged bright-field diffraction profile of a bead imaged over a series of roughly 20 nm-spaced focal steps. The intensity profile (along the*R*-direction of the LUT) is relatively smooth for each step due to the fact that image-noise is azimuthally averaged-out; however, plotting a single radial position as a function of

*z*reveals a significant amount of scatter, which if not accounted for directly translates to scatter in the least-squares residuals in Eq. (14).

Second, this scheme relies on a quadratic fit in order to find the sub-step location. This is problematic because there is no guarantee that the residual versus *Z* is quadratic near the minimum. Moreover, noise also affects the reliability of the quadratic fit. If only three points (${R}^{2}[{z}_{*}]$, ${R}^{2}[{z}_{*}-1]$, and ${R}^{2}[{z}_{*}+1]$) are fit, the resulting parabola will pass through each point, meaning the vertex is susceptible to noise-induced bias. On the other hand, more points will not necessarily follow a parabola and the resulting vertex would not correspond to the best *Z*-estimate.

#### 2.4. Sub-step z-localization using non-linear least-squares

The calculation performed in Eq. (14) is essentially a discrete-version of a non-linear least-squares (NLS) optimization problem. If the LUT, $I[{z}_{j},{r}_{i}]$, is replaced with a set of continuous functions, ${f}_{{r}_{i}}(z)$, which interpolate between the sampled points

#### 2.5. Spline interpolation of LUT

The optimal choice of continuous fitting-functions, Eq. (16), used to describe the LUT requires some consideration. In the absence of a physically motivated equation, the choice of fitting functions seems somewhat arbitrary. Rather than attempt to fit each oscillatory profile function,${f}_{r}(z)$, to a single polynomial, trigonometric series, or any other global function, we use piecewise cubic splines functions. Spline functions have the distinct advantage of being simple to compute while still closely matching the fitted data over the entire range. Generic splines tend to over-fit data, resulting in oscillations between knot-points, which would be both non-physical and problematic for numerical minimization. To overcome this problem, we use Reinsch-type smoothing splines, which are designed to simultaneously minimize fitting error and total curvature, Eq. (34) [21]. Standard implementations of Reinch-type smoothing splines create a knot-point at each unique value of the independent-variable passed during the fit. For LUTs with a large number of samples this creates an excessive number of knots, which increases the computational burden in finding the knot corresponding to the lowest residual. To circumvent this problem, we use a knot-reduction routine to replace the set of smoothing-splines with cubic Hermite splines defined by a minimal set of knots, see Appendix for details. Reduction to Hermite splines was chosen because the Gauss-Newton method only requires that the first derivative be continuous and calculable, which is guaranteed to be true for the entire domain of a Hermite spline. As illustrated in Fig. 1(c) (inset), the knot-reduced smoothing-splines do a good job of capturing the slow oscillation of the data while simultaneously averaging noise.

#### 2.6. Implementation and error estimation of non-linear least-squares

The Gauss-Newton optimization requires an initial guess,${z}_{0}$, to start the calculation. Because the optimal *Z*-location,$\widehat{z}$, is expected to be close to *Z*-value of the row in the LUT with the smallest residual, Eq. (14), that value, ${z}_{j}{}^{*}$, serves as a good starting point. Since splines are simply a set of polynomials, it is straightforward and computationally efficient to evaluate the function and its first-derivative at any point in its domain. As such, the Newton step, Eq. (32), can be computed quickly without numerical approximation.

Similar to the error estimate for the 1-D linear least-squares problem, once the optimal *z*-value is determined, the standard error of the estimate is given by

In summary, the LUT tracking routine starts with discretely sampling the depth-dependency of the radial profile of a particle, Eq. (15). Then, interpolation of the discrete LUT using Reinsch smoothing splines, Eq. (34), simplified by recursive bisection (see appendix, section 6), yields a continuous LUT defined by a set of functions, Eq. (16). During measurements, the images of particles are compared in real-time to the continuous LUT using a Gauss-Newton solver, Eqs. (17) and (30)–(33), which returns the best-fitting *Z*-position. The statistical confidence of that fit is given by Eqs. (18) and (19).

## 3. Results and discussion

#### 3.1. XY noise sensitivity comparison

Localization algorithm performance was assessed using simulated images following the scheme described in Section 2.3. Generated images corresponded to diffraction profiles sampled approximately 1 µm below the apparent focal plane of the particle (as determined by the narrowest “waist” of the diffraction pattern). That position was chosen because it lies within a typical calibration range used in tethered particle experiments [6,22]. Using the simulated images, we directly compared the noise-sensitivity and tracking performance of the radial symmetry routine against the performance of the convolution-based routines discussed in the introduction and outlined in the Appendix [12]. The results of each test are summarized in Fig. 3. Detailed plots of the results for the bright-field tests can be found in Figs. 4 and 5. For each of the comparisons, particles were simulated at different positions within an image and with varying levels of additive noise.

Under moderate to low-noise conditions (SNR>2) the radial symmetry routine outperforms the convolution and Quadrant Interpolation (QI) routines in terms of reliability and versatility. The discrete Fourier transform used in the convolution routine introduces a stepping artifact which results in an oscillatory localization error every ½-pixel increments. (See Fig. 5.) While this artifact is less severe for the QI routine, it does not outperform the radial routine (Fig. 6). Moreover, the convolution and QI routines do not work as particles approach the edge of the field of view. The convolution routines both require the particle’s fringe pattern be completely captured within the image (Fig. 4). Even worse, the QI routine has very poor noise-resiliency once the particle is away from the center of the image by more than a few pixels (Fig. 6). In contrast, for low-to-moderate noise levels the radial symmetry routine has nearly constant performance across the entire field of view, Fig. 4. To be fair, the shortcomings of the convolution and QI routines can be partially overcome by adding an additional processing step wherein image COM is used to roughly approximate the centroid, around-which a new region of interest can be defined. Nevertheless, because the radial symmetry routine uses gradient vectors to determine centroid location, it performs well even when particles lie partially outside the field of view.

#### 3.2. Sensitivity to uneven background and pixel scale

With the choice of Gaussian weighting exponent, *n* = 5 in Eq. (11), the radial symmetry routine is relatively insensitive to shading. To illustrate that point, shading, which in real optical systems could be due to misaligned optics or dust, was simulated by multiplying generated particle images by a 2D wave-like equation defined by

*A*specifies the amplitude of shading, and $\lambda $is a factor proportional to the apparent particle diameter that sets the length-scale of the shading. At shading amplitudes ranging from 0 to 1 and wavelengths varying from two-diameters (~280 px) up to 50-diameters (7000 pixel), locations were computed for 100 particle diffraction patterns simulated at random pixel coordinates in a 250 × 250 grid. Example images and results are summarized in Fig. 7. Even in the situation where the shading varies over distances just twice as large as the particle with intensity variations as large as fringes of the particle (i.e. $\lambda =2\cdot Dia.$and

*A*= 1), the average measurement error was less than one pixel.

#### 3.3. XY localization sensitivity to magnification and focal plane

The symmetry routine also performs well over a range of particle sizes and for beads located on various focal planes. To test the sensitivity to magnification, images were generated as discussed previously, rescaled to varying degrees, and distorted with additive noise defined in Eq. (12). Even with SNR = 1, the localization error remains below 0.06 pixel over effective magnifications ranging from 21× to 63×, see Fig. 8(a).

Similarly, sensitivity to particle focal position was tested by generating diffraction patterns using the fringe pattern stored in the LUT illustrated in Fig. 1(b). The pattern was modified with additive noise and particles were localized. Tracking accuracy is relatively constant across the entire 10 µm range defined in the LUT.

#### 3.4. Gauss-Newton LUT performance and noise sensitivity

*Z*-tracking performance was tested using a similar scheme. A bright-field look up table was recorded by capturing a 10 µm thick *Z*-stack of the 2.8 µm bead used in the *XY* performance comparison. Objective position was controlled using a PI P-721.CDQ piezo objective scanner (Physik Instrumente). The *Z*-stack step size was roughly 20 nm, and three images were captured at each step. The objective scanner includes a close-loop position detector, which reports the actual displacement. Closed-loop measurements reveal differences from the target position due to errors in the voltage output of the piezo controller and thermal drift. The LUT was assembled by indexing each image with the closed-loop readout at the moment of capture. For each image in the *Z*-stack, the particle center was identified using the radial symmetry routine and the radial profile was obtained by computing an intensity versus radius histogram with 1-pixel wide bins. Profiles were normalized by the value of the largest radial bin. Next, the radial profile stack was fit using smoothing splines, as discussed above. The spline LUT was used to generate idealized, noise-free profiles corresponding to any position within the LUT range. Generated profiles were degraded by additive noise. Finally, the degraded profiles were fit using the Gauss-Newton method. Fitting accuracy of our routine was compared to the accuracy of a simple, discrete, closest-match approach, where axial position is determined by which *Z*-step in the discrete LUT yields the lowest residual for Eq. (8), and the quadratic fitting scheme described by in [11], Gosse *et al*. Results are summarized in Fig. 9, and the details are shown in Figs. 10 and 11.

Without noise, the accuracy of the discrete fit is limited to half the step size (10 nm). For moderate noise (SNR> = 5), the Quadratic fitting routine achieves errors that are one third lower than those of discrete-matching. However, as noise increases, the residual, Eq. (14), near the best-fitting step is not well described by a parabola and the performance of the quadratic routine declines significantly. The Gauss-Newton routine, which is agnostic with regard to the shape of the residual curve, outperforms both routines at all noise levels and converges to zero error in the absence of noise (SNR = ∞).

#### 3.5. Sensitivity to particle polydispersity

Up to this point, the LUT routine accuracy was tested using simulated image profiles drawn from the LUT against which they were compared. For real-world applications it would be advantageous to use a single LUT, stored during instrument calibration or immediately prior to experiments, to measure multiple particles. Unfortunately, such a scheme only works if the particles being measured are nearly identical. This is particularly true for bright-field imaging where factors including size, symmetry, and index of refraction all affect the diffraction pattern. To test how particle polydispersity affects LUT fit accuracy, commercially available 1 µm (CP02N Bangs Laboratories, Fischers, IN) and 2.8 µm (M280 Dynabeads, ThermoFisher Scientific) microspheres were non-specifically affixed to coverslips and imaged. The particles have a manufacturer-characterized size range of 1.04 ± 0.05 µm and 2.80 ± 0.08 µm, respectively.

For each sample, a roughly 2 µm deep *z*-stack was captured with 20 nm step spacing. Three images were captured at each step. In each image, particles were located by radial symmetry after which the radial profile was calculated by radially averaging around the local symmetric center for radii ranging from 5 to 100 pixels from the center. To enable comparison of profile patterns from different particles, all profiles were normalized by the average intensity at *R* = 100 pixels. For each particle, profiles from the odd *Z*-step indices (e.g. *Z* = 0.00, 0.04, 0.08 µm…) were fit with smoothing splines, yielding the LUTs. Next, the captured profiles from the even *z*-steps (e.g. *z* = 0.02, 0.06 µm…) were fit to the splines using the Gauss-Newton LUT algorithm.

Results for every particle-LUT pair are summarized in Fig. 12. Error is characterized as the difference between the actual objective position and the position estimated by the LUT algorithm. To account for tilt across the microscope field of view (which would manifest as one particle appearing vertically shifted, for all *Z*-values, when compared to a different particle’s LUT), the average difference between actual and fit positions for a given pair was subtracted from each result. The fit accuracy for thirty-six 1.04 µm particles and nine 2.8 µm particles were compared. For both nominal diameters, fit accuracy for a particle’s image compared with its own LUT was limited to approximately 5-10 nm, corresponding to the repeatability of the P-721.CDQ piezo objective scanner. For 1.04 µm particles, accuracy when comparing a particle’s image with a different LUT was limited to ± 0.1 µm; inter-particle comparisons for the 2.8 µm beads had an accuracy of ± 0.2 µm. These results indicate that even though the particles in each sample are very close in size, subtle differences in their respective diffraction patterns yield unique LUTs which cannot be used interchangeably. The inter-particle discrepancy becomes more pronounced for larger beads, this may be due to the fact that the nominal particle diameter of the smaller beads is relatively close to the wavelength of the illumination source (1.04 µm vs. 640 nm light), while the larger beads are more than 4 wavelengths in diameter. Consequently, high-accuracy localization is best achieved using a LUT captured for each measured particle.

## 4. Conclusions

Look-up-table based particle tracking is an inexpensive and easy-to-implement technique to retrofit nearly any wide-field microscope to locate particles in 3D with nanometer accuracy. The approach presented here includes a sub-pixel *XY* localization routine based on radial symmetry detection which, other than assuming that a particle is radially symmetric, does not require knowledge about how pixel intensity varies as a function of radius (as would be required by traditional PSF-fitting routines). We demonstrate that the routine works well for particles imaged via bright-field, even in the presence of significant image noise. Importantly, the routine also returns traditional confidence limits corresponding to standard error of the estimated position, providing a direct measure of the reliability of each computed location. At moderate to low noise levels (SNR > 2), the radial symmetry routine performs exceptionally well, yielding lateral error < 0.01 pixels. In addition to the lateral (*XY*) algorithms, we also present a noise-robust, look-up-table routine based on a Gauss-Newton, non-linear least-squares solver. By using smoothing splines to interpolate between points of a finitely sampled data set, the LUT routine returns sub-step accuracy limited only by measurement noise. The only requirement is that the values of the LUT vary smoothly between steps. Otherwise it does not depend on the structure of the stored image data, and it works well for both bright-field and fluorescent image stacks.

## 5 Appendix

#### 5.1. Intensity center of mass

For a grayscale image, the intensity center of mass (COM) is estimated by

#### 5.2. Convolution based symmetry detection

Provided the symmetry axes are parallel to the cardinal image dimensions (*x* and *y* axes), the symmetric origin can be estimated by finding the peak of the cross-correlation of the image with its mirror image [11]. In one dimension this corresponds to finding

At the expense of increasing the computational complexity or the inverse-FT, the accuracy can be improved via sinc-interpolation. For example, quarter-pixel accuracy can be achieved by padding the FT with *N*-zeros in the *x*-direction and *M*-zeros in the *y*-direction, increasing the computation complexity to roughly $4MN{\mathrm{log}}_{2}2MN$operations. An efficient alternative to sinc-interpolation is to locally fit $C[x,y]$with a paraboloid and use the vertex as a sub-integer estimate of $[x*,y*]$. However, the true shape of $C[x,y]$depends on the image and is not guaranteed to be quadratic near its maximum.

For an *M*-tall by *N*-wide image, the origin of symmetry can be determined by convolution. Assuming 1-indexing, the symmetric center $({x}_{c},{y}_{c})$is determined by

*N*and

*M*respectively.

If sinc-interpolation is used to increase the resolution of the fit, the Fourier transform is zero-padded and the scaling-relationship in Eq. (24) changes slightly. Quarter-pixel accuracy corresponds to padding the FT with *N*-zeros in the *x*-direction and *M*-zeros in the *y*-direction, and then applying the inverse-FT to compute the convolution. The origin is then given by

#### 5.3. Smoothed image gradient

To minimize the impact of additional noise resulting from finite differencing, we calculate the image gradient along ± 45°, with $\widehat{u}=\frac{1}{\sqrt{2}}(\widehat{x}+\widehat{y})$and $\widehat{\nu}=\frac{1}{\sqrt{2}}(\widehat{y}-\widehat{x})$and apply a 3 × 3 smoothing filter, so

#### 5.4. Gauss-Newton method

Starting from an initial guess, ${z}_{0}$, the solution,$\widehat{z}$, that minimizes the objective function in Eq. (16) is iteratively found by computing the Newton step

where ${{\rm H}}_{\varphi}$is the Hessian of $\varphi $, and then updating the guessThe gradient of the objective can be written in terms of the Jacobian, ${J}_{i}(z)=\frac{\partial {R}_{i}(z)}{\partial z}$,For small residuals (i.e. good guesses) the Hessian converges to ${H}_{\varphi}=J{(z)}^{T}J(z)$. The equation defining the Newton step reduces toIteration is terminated when the magnitude of the residual is less than a specified tolerance, $R{(z)}^{T}R(z)\le \text{TOL}\text{.}$#### 5.5. Smoothing splines

Reinsch-type smoothing splines are defined by minimizing the function

*p*is a user-specified parameter which determines the relative weighting of the two factors. We chose$p={(1+{(\Delta z/({Z}_{\mathrm{max}}-{Z}_{\mathrm{min}}))}^{3})}^{-1}$, where$\Delta z$is the z-step spacing, relative to the range of measured axial positions, of the LUT, because it qualitatively balances between over-fitting and over-smoothing.

#### 5.6. Hermite spline knot reduction

Smoothing splines are simplified by recursive bisection, similar to the strategy outlined by Dung & Tjahjowidodo [23]. First, the densely-knotted target spline is evaluated at 5 evenly spaced points spanning its range. The resulting values are fit with a piecewise cubic Hermite spline and divided in half, creating two segments. For each segment, the 2-norm residual difference is calculated between each knot in the original spline and resulting simplified Hermite spline. If the average residual is above a given threshold, the segment is sub-divided and the process repeated.

## Funding

National Institute of General Medical Sciences, National Institutes of Health (NIH) (R01 GM084070).

## Acknowledgments

We thank Cees Dekker and Jacob Kerssemakers for providing us with a MATLAB implementation of their Quadrant Interpolation tracking algorithm, enabling a direct performance comparison of their routines with our own. We also thank James Nagy for helpful discussions. All of the algorithms presented here, excluding the QI routine, are open source and are included in Code 1 [24]. As of publication, the radial symmetry and spline LUT routines can also be found at https://github.com/dkovari/ExtrasToolbox.

## References

**1. **T. Plénat, C. Tardin, P. Rousseau, and L. Salomé, “High-throughput single-molecule analysis of DNA-protein interactions by tethered particle motion,” Nucleic Acids Res. **40**(12), e89 (2012). [CrossRef] [PubMed]

**2. **P. Daldrop, H. Brutzer, A. Huhle, D. J. Kauert, and R. Seidel, “Extending the range for force calibration in magnetic tweezers,” Biophys. J. **108**(10), 2550–2561 (2015). [CrossRef] [PubMed]

**3. **I. De Vlaminck and C. Dekker, “Recent Advances In Magnetic Tweezers,” Annu. Rev. Biophys. **41**(1), 453–472 (2012). [CrossRef] [PubMed]

**4. **D. T. Kovari, Y. Yan, L. Finzi, and D. Dunlap, “Tethered Particle Motion: An Easy Technique for Probing DNA Topology and Interactions with Transcription Factors,” in *Single Molecule Analysis. Methods in Molecular Biology* (Springer, 2018) pp. 317–340.

**5. **M. Bierbaum, B. D. Leahy, A. A. Alemi, I. Cohen, and J. P. Sethna, “Light microscopy at maximal precision,” Phys. Rev. X **7**(4), 1–10 (2017). [CrossRef]

**6. **Y. Seol and K. C. Neuman, “Magnetic tweezers for single-molecule manipulation,” Methods Mol. Biol. **783**, 265–293 (2011). [CrossRef] [PubMed]

**7. **G. Sitters, D. Kamsma, G. Thalhammer, M. Ritsch-Marte, E. J. G. Peterman, and G. J. L. Wuite, “Acoustic force spectroscopy,” Nat. Methods **12**(1), 47–50 (2015). [CrossRef] [PubMed]

**8. **B. M. Lansdorp and O. A. Saleh, “Power spectrum and Allan variance methods for calibrating single-molecule video-tracking instruments,” Rev. Sci. Instrum. **83**(2), 025115 (2012). [CrossRef] [PubMed]

**9. **B. M. Lansdorp, S. J. Tabrizi, A. Dittmore, and O. A. Saleh, “A high-speed magnetic tweezer beyond 10,000 frames per second,” Rev. Sci. Instrum. **84**(4), 044301 (2013). [CrossRef] [PubMed]

**10. **R. Parthasarathy, “Rapid, accurate particle tracking by calculation of radial symmetry centers,” Nat. Methods **9**(7), 724–726 (2012). [CrossRef] [PubMed]

**11. **C. Gosse and V. Croquette, “Magnetic tweezers: micromanipulation and force measurement at the molecular level,” Biophys. J. **82**(6), 3314–3329 (2002). [CrossRef] [PubMed]

**12. **M. T. J. van Loenhout, J. W. J. Kerssemakers, I. De Vlaminck, and C. Dekker, “Non-bias-limited tracking of spherical particles, enabling nanometer resolution at low magnification,” Biophys. J. **102**(10), 2362–2371 (2012). [CrossRef] [PubMed]

**13. **S.-H. Lee, Y. Roichman, G.-R. Yi, S.-H. Kim, S.-M. Yang, A. van Blaaderen, P. van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Opt. Express **15**(26), 18275–18282 (2007). [CrossRef] [PubMed]

**14. **J. C. Crocker and D. G. Grier, “Methods of Digital Video Microscopy for Colloidal Studies,” J. Colloid Interface Sci. **179**(1), 298–310 (1996). [CrossRef]

**15. **I. D. Vilfan, J. Lipfert, D. A. Koster, S. G. Lemay, and N. H. Dekker, “Magnetic Tweezers for Single-Molecule Experiments,” in *Handbook of Single-Molecule Biophysics* (Springer US, 2009), pp. 371–395.

**16. **A. von Diezmann, Y. Shechtman, and W. E. Moerner, “Three-Dimensional Localization of Single Molecules for Super-Resolution Imaging and Single-Particle Tracking,” Chem. Rev. **117**(11), 7244–7275 (2017). [CrossRef] [PubMed]

**17. **T. Lionnet, J. F. Allemand, A. Revyakin, T. R. Strick, O. A. Saleh, D. Bensimon, and V. Croquette, “Single-molecule studies using magnetic traps,” Cold Spring Harb. Protoc. **2012**(1), 34–49 (2012). [CrossRef] [PubMed]

**18. **N. Ribeck and O. A. Saleh, “Multiplexed single-molecule measurements with magnetic tweezers,” Rev. Sci. Instrum. **79**(9), 094301 (2008). [CrossRef] [PubMed]

**19. **S. F. Nørrelykke and H. Flyvbjerg, “Power spectrum analysis with least-squares fitting: Amplitude bias and its elimination, with application to optical tweezers and atomic force microscope cantilevers,” Rev. Sci. Instrum. **81**(7), 075103 (2010). [CrossRef] [PubMed]

**20. **M. T. Heath, *Scientific Computing: An Introductory Survey*, 2nd ed. (McGraw-Hill, 2002).

**21. **C. H. Reinsch, “Smoothing by spline functions, II,” Numer. Math. **16**(5), 451–454 (1971). [CrossRef]

**22. **Z. Vörös, Y. Yan, D. T. Kovari, L. Finzi, and D. Dunlap, “Proteins mediating DNA loops effectively block transcription,” Protein Sci. **26**(7), 1427–1438 (2017). [CrossRef] [PubMed]

**23. **V. T. Dung and T. Tjahjowidodo, “A direct method to solve optimal knots of B-spline curves: An application for non-uniform B-spline curves fitting,” PLoS One **12**(3), e0173857 (2017). [CrossRef] [PubMed]

**24. **D. T. Kovari, “Radial Symmetry and Spline Look-up-table Inversion” figshare (2019) [retrieved 28 Jun. 2019], https://figshare.com/articles/Radial_Symmetry_LUT_Tracking_Code/8317136.