Robust and highly performant ring detection algorithm for 3d particle tracking using 2d microscope imaging

Three-dimensional particle tracking is an essential tool in studying dynamics under the microscope, namely, fluid dynamics in microfluidic devices, bacteria taxis, cellular trafficking. The 3d position can be determined using 2d imaging alone by measuring the diffraction rings generated by an out-of-focus fluorescent particle, imaged on a single camera. Here I present a ring detection algorithm exhibiting a high detection rate, which is robust to the challenges arising from ring occlusion, inclusions and overlaps, and allows resolving particles even when near to each other. It is capable of real time analysis thanks to its high performance and low memory footprint. The proposed algorithm, an offspring of the circle Hough transform, addresses the need to efficiently trace the trajectories of many particles concurrently, when their number in not necessarily fixed, by solving a classification problem, and overcomes the challenges of finding local maxima in the complex parameter space which results from ring clusters and noise. Several algorithmic concepts introduced here can be advantageous in other cases, particularly when dealing with noisy and sparse data. The implementation is based on open-source and cross-platform software packages only, making it easy to distribute and modify. It is implemented in a microfluidic experiment allowing real-time multi-particle tracking at 70 Hz, achieving a detection rate which exceeds 94% and only 1% false-detection.


Introduction
In the process of tracking particles in a microfluidic experiment, I was confronted with the need to detect and localise particles in 3d, based on 0.7 megapixel images, acquired at 70Hz, containing roughly 50 particles per frame.However, my needs are not unique because the study of dynamics often relies on tracking objects under the microscope.Indeed, precise and robust particle tracking is essential in many fields, including studies of micro-Rheology [1,2], chaotic dissipative flows [3], feedback for micromanipulation [4], and other soft condensed matter physics and engineering problems.Moreover, microfluidic systems play a growing role as part of lab-on-a-chip apparatus in micro-chemistry [5,6], bioanalytics [7,8], and other bio-medical research and engineering applications [9].Yet, detailed characterisation of the flow and transport phenomena at the micro-scale is still a non-trivial task.In general, the motion is three dimensional and automated tracking is cumbersome from the perspective of both instrumentation and software.Setting up several viewing angles as done for large systems [10]becomes even more complicated in microscopic systems, while scanning through the third axis, e.g.confocal microscopy, clearly compromises temporal resolution and concurrency.
When tracking small light emitting objects, such as fluorescent particles under the microscope, the appearance of rings is often a sign of the object going out of focus.Normally this results in the loss of the tracked object, which is thereafter considered as a hindering background source.However, these rings carry information of the 3-dimensional position of the particle.This has been used for localising a single light scattering magnetic bead based on matching the radial intensity profile to an empirical set of reference images [4].An axial range of 10µm was demonstrated and a temporal resolution of 25Hz was achieved using the knowledge of the particle's previous position.In fact, for fluorescent particles the radius of the most visible ring of each particle precisely indicates its axial position -the radius follows a simple scaling with the particle distance from the focal plane (see Fig. 3 in Appendix B).A similar approach was recently described in [11], where the measurements were, once again, limited to a single particle in the observation volume, with an axial range of 3µm and temporal resolution of 10Hz.
The development of the method presented here was motivated by the study of pair dispersion in a chaotic flow [12], taking place in a microfluidic tube of 140µm -the observation volume is larger by more than three orders of magnitude with respect to those reported in Refs.[4] and [11].The experiments consist of tracking tracers advected by the flow, at seeding levels of several tens to hundreds in the observation window, where it is necessary that the particles are resolved even when nearby to each other.The typical flow rates dictate sampling rates of 70Hz whereas the statistical nature of the problem requires data acquisition over weeks.Using a standard epifluorescence microscope and the fact that the parameters of the most visible ring can be mapped to the 3d position of the tracer, the particle localisation problem is converted to a circle detection problem.The constraints set by the nature of the experiment require an image analysis algorithm that is robust not only to the noise of the image acquisition process, but to rings overlaps, inclusions and occlusions as well.The typical complexity of the images is exemplified in a sub-frame from our experiment presented in Fig. 1(a).In addition, the data flow is of about 180GB per hour, an overwhelming rate which demands the optimisation of the algorithm for realtime analysis.The key steps for achieving high-performance are introduced following the presentation of the main concepts which contribute to the robustness of the algorithm.
Imagine for the moment that you have successfully identified which of the pixels in the image reside on a ring.The issue of doing so will be addressed later on.Given this set of coordinates, it may seem straightforward to find the parameters of the circles which best fit them.However there is a missing piece of information here, that is, which sub-sets of pixels belong together to form a ring.Moreover, we do not know a priori how many rings there are in the image and there may be false detected coordinates which we would like to disregard.Therefore, we need some method to classify/cluster the coordinates into sub-sets, each sub-set matching a single ring.
A circle in a two-dimensional image is uniquely specified by three parameters.In this work two designate the centre of the circle and the third specifies its radius.The parameter space of all possible circles is therefore three-dimensional.One can detect circles in an image by mapping the image intensity field to the circle parameter space.Peaks in this parameter space imply a circle well represented in the image.One approach to achieve this mapping is via a discrete Radon transform, which for the purposes of this presentation translates to convolving the image with a mask of a ring [13].Since each candidate radius calls for a separate convolution, this results in visiting all the pixels in the image over and over.Recall that the outermost ring is sufficient for 3d localisation of the imaged particles.Hence it is worth noting that even at moderate rings densities, the pixels lying on the outer-most ones consist a small fraction of the pixel population, less than 2% in my case.When there are more than a couple of potential radii this procedure would perform a plethora of useless computations [13].This fact directs to another approach, which may seem equivalent yet explicitly exploits this information sparsity -the circle Hough transform [14]: each pixel votes for all the candidate points in the parameter space of which it may be part.In this way, every pixel is visited once, potentially reducing the computation time by orders of magnitude.The discrete version of the parameter space is commonly referred to as the array of accumulators.During the voting procedure each vote increments an accumulator by one.Alas, in the literature of Computer Vision and Pattern Recognition it is well known that the standard circle Hough transform is rather demanding both for large memory requirements, which grow with the radii range, as well as for its 3d nature which renders peak finding in the parameter space a difficult task to tackle [15,16].
One way to address these computational challenges is to resort to lower dimensionality circle Hough transforms, but these usually miss circles having nearby centres and are less robust, resulting in higher false positive and false negative errors rates [15].Another path is to randomly sub-sample the information content in the image, giving way to non-deterministic methods; see [16] and references therein.However, due to their random nature these methods suffer from inferior detection rates and accuracy when compared with the deterministic ones [16].
For these reasons I developed an algorithm which is an offspring of the full 3d-circle Hough transform, yet the local maxima detection issues are addressed and it shows high performance and a small memory footprint.

Results
The key steps of the algorithm are conceptually outlined as follows: (i) detect directed ridges; (ii) map the directed ridges to the parameter space of circles; (iii) detect local maxima via radius dependent smoothing and normalisation; (iv) classify the coordinates of the ridge pixels according to the peaks in the circle parameter space, and fit each sub-set to a circle, achieving sub-pixel accuracy.This outline is presented graphically in Fig. 2 for a small sub-frame containing two fluorescent particles.The first step is locating the pixels of interest.Like many other feature detection algorithms, the standard circle Hough transform relies on an edge detection step, where edges are the borders of dark and bright regions.As the images contain rings rather then filled circles, I chose to implement   an algorithm that detects ridges, thin curves which are brighter than their neighbourhood, rather than edges.This exhibits better consistency.First note that the ring of interest is thicker than the inner ones.This is advantageous as the image admits scale selection [17] -the inner rings can be suppressed using a Gaussian smoothing having the appropriate scale, approximately that of the most visible ring.Ridges are then found using a differential geometric descriptor [17], which defines ridge pixels using the following two properties: (i) Negative least principal curvature, k − < 0; (ii) k − is a local minimum along the direction of the associated principal direction X − .The principal curvatures are the eigenvalues of the Hessian, the matrix of the image second spatial derivatives.Here k − and X − denote the smaller principal curvature and the corresponding eigenvector.This is demonstrated for the two fluorescent particles in Fig. 2(a), where this sub-image is analysed for directed ridges, represented by arrows overlaid on the k − field as background of Fig. 2(b).Note that X − is collinear with the direction to the centre of the ring.I use this to significantly reduce the complexity of the voting procedure, in a similar way to the gradient directed circle Hough transform [18] -each directed ridge pixel votes for all candidate points in the 3d parameter space of circles, provided that the coordinates of the circle centre is within the r min to r max range, directed along X − .
As votes from the ridge pixels accumulate, each ring in the image transforms into two mirroring coaxial cones, aligned along the radius axis, having a joint apex.This procedure results in a discrete scalar function over a 3d box.This is demonstrated in Fig. 2(c).The coordinates of the apexes, which are the local maxima of this function, are the candidate circle parameters.In practice, there are many sources which render the resulting circle parameter space very noisy: The raw image is a discrete representation of the intensity field, the image acquisition process itself is not noiseless, and the image complexity mentioned above, all may result in errors in the detected ridge position and direction, as well as false ridge detection and false negatives.Note that the deviation from ideal voting due to the error in determining the ridge direction grows linearly with the radius.Therefore, each equi-radius level of the ring parameter space is smoothed using Gaussian weights, whose width is proportional to the radius.Next we note that finding local maxima in the 3d parameter space requires the comparison of accumulators in different equi-radius levels, which asks for some normalisation as larger rings are expected to receive more votes.This leads to a natural normalisation by 1/r, following which an ideal ring is expected to receive 2π votes.Fig. 2(d ) shows how this procedure simplifies the parameter space.Local maxima can then be located by nearest-neighbours comparison.
The local maxima identified in the parameter space induce a classification on the ridge coordinates: for each peak an annulus mask is formed and a best fitting circle is found for all ridge coordinates within the annulus.Ridge pixels which are not covered by any annulus mask are not fitted for.This is desired as these usually result from non-circular features in the image or noise.See the example in Fig. 2(e).In this way sub-pixel precision is achieved.
Examining the robustness of the algorithm on experimental data reveals a detection rate that exceeds 94% with only 1% false-detection; for further details see Robustness assessment in the Methods section.
As was mentioned above, it is desired for our purposes to have the images processed in real-time.Here I briefly outline the key ideas behind the optimisation of the algorithm, the full details are available in the open-source code itself [19].The first key point for the algorithm optimisation is the splitting of the voting procedure -the votes are recorded for each ridge pixel as it is detected, such that ridge detection and votes collection are done in one-pass.
The population of the parameter space is performed at a separate stage, which leads to the second key point.Instead of holding an array for the full parameter space, only two sub-spaces are maintained, consisting of three consecutive equi-radius levels; the first for the raw parameter sub-space, the second for the smoothed and normalised one, where local maxima are searched for.The equi-radius levels are populated and processed one by one.Only the accumulators exceeding an integer vote threshold are regarded as hotspots and are mapped to the smoothed and normalised sub-space.Each time a radius-level is completed, regarded here as the "top" one, the hotspots in the level beneath, the "middle" one, are verified to exceed a preset floating point threshold, a fraction of 2π.Those which do are searched for local maxima by a nearest neighbour comparison within a 3×3×3 voxels box.Once this search is completed, the "bottom" level is no longer needed and a cyclic permutation takes place where the "bottom" level becomes the new "top".This allows memory recycling and avoids the need to initialise big arrays of zeros to represent the full 3d parameter space.
Registering modified array elements and undoing is the last key point.The radius-levels are required to be blank prior to their population.Re-calling the sparsity of the parameter space (see Figs. 2(c) and 2(d ) for example), going over all its elements is a waste of processing time.Instead, each time an array element is modified for the first time its indices are registered.Once the search for peaks is done, all modifications to the "bottom" radius-level are undone, preparing it for reuse.This lifts the need to clean the whole array.
The combination of memory recycling with registering modified array elements and undoing reduces the computation time and in my case results in a nearly ten times less memory consumption.In fact, the size of the arrays representing the parameter space kept in memory is now fixed and no longer grows with the radii range.For preliminary testing purposes, I first implemented the algorithm outlined in the beginning of the Results section (as well as in Fig. 2) using convolutions and other array based operations.The final implementation, inspired by the circle Hough transform and including the above optimisations, is more than 50 times faster.This is attributable to the reduction in the number of operations required once the sparsity of the data is taken advantage of.Further details and explanations can be found in the Methods section and Appendix A .

Discussion
In comparison with other existing methods for 3d particle tracking, the method presented here is advantageous when it comes to long measurements, temporal resolution and concurrency, as well as real-time applications.The confocal scanning microscope requires scanning the volume of interest.Therefore it is slower and cannot provide instantaneous information of the whole volume.Unlike Holographic microscopy [20,21], the proposed method does not pose long and heavy computational demands which is restrictive for real-time applications or when large datasets are required for statistics.
One could expect the optical method discussed here to produce patterns which are symmetric about the focal plane.When this applies, it may result in an ambiguity with respect to whether the particle is above or below focus.Our optical arrangement (see the Methods section) shows clear diffraction rings only on one side.Furthermore, as particles approach focus, the radius of the outer-most ring becomes too small to resolve.For these reasons the focal plane is placed outside the volume of interest (as reflected in Fig. 3 in Appendix B).Optical astigmatism offers a mean for discriminating between the two sides of the optical axis [22,23,24].The introduction of a cylindrical lens results in the deformation of a circular spot into an ellipsoidal one as a fluorescent particle goes further away from focus, with the ellipse major axis of a particle above focus aligned perpendicular to a one below.In [22], the axial range was limited to a couple of microns above and below focus; in [23,24] it was restricted to less than a micron.Within these ranges the tracer image can be approximated by an elliptical gaussian pattern.However, extending the range generates elliptical rings as well [22, Fig. 1].This requires dealing with two species of patterns, spots and rings.Moreover, deforming circular rings into elliptical ones, the dimensionality of the parameter space increases, and so does the technical complexity of the image analysis.Therefore the advantage of the stronger signal, by working closer to focus on both its sides, is expected to have a heavy computational cost once the range is extended such that diffraction rings appear as well.The method presented here requires working away from focus.Rings visibility decreases as the fluorescence signal spreads over a larger area, thus setting the lower bound for the exposure time.Nevertheless, I have found that the fluorescence signal-to-noise ratio allowed tracking particles moving chaotically at speeds exceeding 400µm/s.
One circle detection alternative to the one proposed here is the EDCircles algorithm described in [25].It detects contiguous edge segments and employs the Helmholtz principle for controlling false detections.It is attractive for two main reasons -the mathematical theory of perception it is founded on [26], and the promising results presented in the manuscript [25] promoting it as a general purpose circle detector.Yet, for the purpose of this work the EDCircles algorithm did not perform as well as the one presented here (see Robustness assessment in the Methods section).
Many of the algorithmic concepts introduced above can be advantageous when implemented in other algorithms if the data is sparse, namely memory recycling as well as registering modified array elements and undoing.Though this algorithm development was motivated by the analysis of fluorescence microscopy images, it is more general and can be applied to other cases as well.The interpretation of the Hough transform as a classification/clustering algorithm has a wider potential than merely image analysis.To name a physical example is the case of particle jets emerging from several sources of unknown loci.A dataset consisting of the positions and momenta of the particles at a certain time is analogous to that of the directed ridges and the jet sources can be identified with the local maxima over the parameter space.
The method I introduced above is currently implemented in an experiment which requires long unsupervised measurements lasting for days at high temporal resolution, sampling a volume of interest which contains tens to hundreds of particles.Thanks to the fact that the whole volume of interest is sampled at once by a single 2d image, concurrency is achieved.The use of LED and the relatively short exposure times can be potentially exploited to avoid photo-damage and bleaching.It is demonstrated to be robust to the overlap, inclusion and occlusion of the ring pattern of the imaged particles.It features high performance admitting real-time applications.
This method paves the way for studies of 3d flows in microfluidic devices.Its robustness to the vicinity of particles to each other allows to study the dynamics of particle pairs [12], triples, etc.As such it also has a potential for biomedical research.A possible immediate application is the detailed characterisation of the transport induced by the presence of cells in confined flows, a phenomenon presented in [27].Together with the development of high signal tracers, labelling techniques and sensitive cameras, this method may be useful in other life sciences studies such as cellular trafficking [28], cell migration [29] and bacterial taxis [30].

Methods
Algorithm implementation I implemented this algorithm relying on freely available open-source and cross-platform software packages only.The source is available at [19] .Most of the heavy lifting is achieved using the Cython language [31].It has a Python-like syntax from which a C code is automatically generated and compiled.This allows the code to be short and easy to read while enjoying the performance of C. For example, this implementation exploits the Numpy/Cython strided direct data access [31,32] by fully sorting the votes.In the image pre-processing step the image is smoothed using a Gaussian convolution and the smoothed image spatial derivatives are calculated using a 5×5 2 nd order Sobel operator; these operations are done using OpenCV's Python bindings [33].
The equi-radius levels of the circle parameter space are smoothed using Gaussian weights exp − 1 2 x−x hotspot σ

2
, whose width is proportional to the radius σ(r) ∝ r.The explicit form, σ(r) = 0.05r + 0.25, was found empirically.The slope coefficient is interpreted as accounting for ∼ 0.1rad uncertainty in the direction of X − .

Experimental details
The imaging system consists of an inverted fluorescence microscope (IMT-2, Olympus), mounted with a Plan-Apochromat 20×/0.8NAobjective (Carl Zeiss) and a fluorescence filter cube; a Royal-Blue LED (Luxeonstar) served for the fluorophore excitation.A CCD (GX1920, Allied Vision Technologies) was mounted via zoom and 0.1× cmount adapters (Vario-Orthomate 543513 and 543431, Leitz), sampling at 70Hz, 968px × 728px, covering 810 × 610µm 2 laterally.The experiments were conducted in a microfluidic device, implemented in polydimethylsiloxane elastomer by soft lithography, consisting of a curvilinear tube (see grey broken line in Fig. 4(a)).The rectangular cross-section of the tube was measured to be of 140µm depth and 185µm width.The working fluid consisted of polyacrylamide in aqueous sugar (sucrose and sorbitol) syrup, seeded with fluorescent particles (1 micron 15702 Fluoresbrite R YG Carboxylate particles, PolySciences Inc.).The flow was driven by gravity.
For the empirical calibration, the same working fluid was sandwiched between two microscope glass slides.The separation distance between the slides was set to 161µm by micro-spheres (4316A PS NIST certified calibration and traceability, Duke Standards) serving as spacers.The microscope objective was translated in steps of 2µm to acquire images of the tracers at different off-focus distances ∆z.The microscope focus knob was manipulated by a computer controlling a stepper motor.During the rest stages of the objective, the ring radii of every detected tracer were averaged over 210 frames spanning 3s.Due to the high viscosity of the fluid, 1100 times larger than water viscosity, tracer motion due to diffusion is negligible during this time interval.The median of the estimated standard-deviations of the data presented in Fig. 3(a) is 0.03px and the maximal is 0.27px.In practice, to account for the uncertainties in finding the focal position and due to optical aberrations, 25 tracers dispersed throughout the observation volume were accounted for.Their curves were aligned via shifting ∆z by the larger root of a quadratic polynomial fit.Then, the conversion function r −1 (r) was obtained by inversion of the quadratic polynomial fit accounting for all the data together; see Fig. 3(b).The resulting root-mean-squarederror, (∆z − r −1 (r)) 2 = 1.97µm, and the maximal measured absolute error is 5.35µm; these reflect the uncertainty due to the empirical calibration procedure taken here.Finally, the out-of-focus distance of the objective ∆z has to be converted to the physical distance via multiplication by the ratio of the refractive indices, 1.58 in this case.The observed axial range exceeds 180µm.

Robustness assessment
In order to estimate the robustness of the algorithm, images from the chaotic flow experiment were analysed and cropped to a sub-region of 340px×370px.The analysis results of 600 such subframes were examined.These sub-frames contained 14 rings on average, out of which 67.7% were in a ring cluster (overlap and inclusion configurations).This examination shows an average of 6.8% False-Negative errors.In some sub-frames a ghost ring would appear accompanied by a strong distortion of the tracer image in its real position.This is attributed the microfluidic walls and observed only when a particle is very close to the wall.Excluding from the statistics two such tracers, the False-Negative error rate is reduced, corresponding to a detection rate of 94.7%.This examination does not show any significant sensitivity to rings overlap and inclusion.On the contrary, 95.5% of the rings in clusters were identified correctly while only 0.8% of the reported rings in clusters were non-existing particles.
To illustrate the high-quality of the above results, Fig. 4(a) was analysed by the on-line demo of the parameter-free EDCircles algorithm, provided by the authors of [25].The results are provided here for comparison.The EDCircles demo detected 35 circles and 3 ellipses.One of the ellipses is a false merging of two rings, which were correctly resolved by the algorithm presented here.Excluding these two, the EDCircles did not detect 19 out of the other 56 rings found by the method proposed in this manuscript.None of the particles correctly reported by the former was missed by the latter.
Performance assessment The performance assessment is based the analysis of 1500 full frames containing 50 particles on average; for a typical example see Fig. 4(a).The test was run on an i7-3820 CPU desktop, running Ubuntu linux operating system.A single process analysed at an average rate of 6.28Hz.To achieve higher performance as required for our experiments, I use the multiprocessing package of Python, exploiting the multi-core processors available on modern computers.The analysis rate scales linearly with the number of processes.No deterioration of the processing rate per core was noted (tested up to twice the core number with hyper-threading).Based on a producer-consumer model, one can even transparently distribute the workload among several computers if needed.This is partially attributable to the small memory footprint of the algorithm.
Precision assessment To estimate the precision of the presented localisation method, smoothing splines were applied to the reconstructed trajectories providing an estimator for the error variance σ 2 [34].The mean values are as follows: σ x = 0.134µm, σ y = 0.135µm and σ z = 0.434µm, for x and y denoting the lateral coordinates, and z the axial one.This axial uncertainty corresponds to 0.3% of the axial range covered by the particles.According to the data provided in Ref. [4] a value of 0.1% was achieved and a similar one in Ref. [11].The uncertainty estimates reported here account for noise which rise not only due to the image analysis and the multi-particle scenario, but also due to other sources, namely the motion of the particles and the linking procedure.The details are as follows.
From a 3m30s measurement, 2014 trajectories which span more than 1s were analysed (discarding shorter ones).Particle positions, converted to microns, were linked to reconstruct their trajectories; for a sub-sample see Fig. 4(b).The linking algorithm was adapted from the code accompanying Ref. [35], generalised to n-dimensions, the kinematic model was modified to account for accelerations, and a memory feature was introduced to account for occasional misses.Finally, a natural cubic smoothing spline was applied to smooth-out noise and obtain an estimate of particle velocity and acceleration [36,37].The smoothing parameter was automatically set using Vapnik's measure, using a code adapted from the Octave splines package [38].

Appendix A Detailed algorithm
As mentioned in the main text, the standard circle Hough transform is often avoided not only for its challenging local maximum detection in noisy 3d space but for its heavy memory requirements as well.The standard circle Hough transform requires a 3-dimensional array of accumulators.The coordinates of each array element are the parameters of a candidate circle.The value of the accumulator at these coordinates indicates how well this circle is represented in the image.Code optimisation for high-performance and small memory footprint is achieved following this scheme: 1. Image pre-processing step: the image is smoothed using a Gaussian convolution and the smoothed image spatial derivatives are calculated using a 5×5 2 nd order Sobel operator [33].Using these derivatives, the local least principal curvature k − is estimated as the smaller eigenvalue of the Hessian matrix.
2. One-pass ridge detection and votes collection: for each pixel in k − which is smaller than a pre-defined curvature threshold (the latter is no greater than zero), the corresponding X − is calculated.If this pixel is found to be a local minimum along the direction of X − , its coordinates are recorded in the ridge container.At this stage, its votes are collected as well, that is, the potential circles parameters to which it may belong.
3. Sort the votes stack according to the radii: this allows performing the parameter space incrementing procedure equi-radius level by level.In order to achieve higher performance, the votes are further sorted by the row index and then by the column index exploiting the numpy/cython strided direct data access [31,32].For this reason each circle parameter triple is represented as an integer using a bijection.

4.
Circle parameter space population and local maximum detection via radius-dependent smoothing and normalisation: This is done using two arrays representing a sub-space of the full circle parameter space.Each consists of 3 consecutive equi-radius levels; the first for the raw accumulators sub-space, the second for the smoothed and normalised one, where local maxima are to be searched for.There are two votes thresholding steps: an integer threshold for the raw accumulators and a floating point threshold, a fraction of 2π, for the smoothed and normalised array elements.In describing the procedure it is assumed that the r − 2 and r − 1 levels have already been populated in both sub-space triples and the r levels are blank, i.e. all zeros.As long as the votes drawn from the votes stack point to the same radius, the corresponding radius-level is populated by incrementing the indicated accumulator.Recall that the votes are fully sorted hence all votes pointing to a certain voxel will come out from the stack in a row.Every time a new circle parameter triple is encountered, its coordinates are recorded as modified.In case the previously incremented voxel has surpassed the 1 st votes threshold its coordinates are recorded as a hotspot -a circle candidate.Once there are no more votes for this r-level, it is mapped to the second subspace: for each hotspot voxel a spatial average is calculated, weighted by a Gaussian function, which width is linearly dependent on the radius; the value of the average is then normalised by 1/r.After mapping all the hostpots of the current r-level, a local maximum is searched for among the hotspots of the (r − 1)-level which pass the 2 nd votes threshold.This is done using a nearest neighbours comparison within a 3 × 3 × 3 voxels box.Array elements which are local maximum and exceed the threshold are registered as rings.Once all hotspots have been processed, all modifications to the (r − 2)-level are undone as its data are no longer needed.By this it is made ready to be regarded as the next r-level and a cyclic permutation among the levels takes place.In practice, this is performed by accessing the equi-radius levels using the modulo operation -the radius indices are calculated using r (mod 3).

Figure 1 :Figure 2 :
Figure 1: Snapshots from the experiment and a demonstration of the algorithm robustness: (a) typical image complexity is exemplified in an unprocessed sub-frame consisting of 1/9 part of the full frame, corresponding to lateral dimension of 215 × 315µm 2 .The axial range available for the particles is 140µm.(b) the corresponding analysis result; in red are the radii in pixels units.(c) & (d ) time sequences of sub-frames (400ms each).Red coloured particles in (c) demonstrate pair dispersion, in which the algorithm is required to resolve rings with similar parameters.The yellow particle in (d ) shows radius change corresponding to a downwards translation.Each sub-frame in (c) & (d ) images a box which lateral dimensions is 190 × 270µm 2

5 .Figure 4 :
Figure 4: 2d out-of-focus images → 3d trajectories: (a) A single raw full 2d frame imaging a volume of 810µm × 610µm × 140µm, including 60 tracers in one period of the curvilinear tube (which boundaries are denoted by the broken grey line).Smaller rings result from tracers being closer to the focal plane, which is positioned above the tube.(b) A subsample of 40 trajectories reconstructed from a time-lapse sequence of such frames, which spans 12 seconds of data acquisition.The colour coding indicates the speed ranging from 80µm/s (blue) to 400µm/s (red); the mean flow is rightwards.Corresponding axes are denoted by colours.The isometric view of the bottom panel can be obtained by three rotations, starting from the orientation of the upper panel: −90 • about the green axis, −45 • about the blue axis, and approximately 35 • about the new horizontal axis.