1.

Introduction

Superresolution localization microscopy, such as photoactivated localization microscopy,¹ stochastic optical reconstruction microscopy,² and fluorescence photoactivation localization microscopy,³ can achieve nanometer resolution using photoswitching to image fluorophores sequentially and localization algorithms to measure their positions precisely. The performance of localization algorithms is a key factor in determining superresolution imaging performance.^4,5 Most of these algorithms are based on fitting the point spread function (PSF) of a single fluorophore with a known mathematical model. Alternatives to fitting approaches include centroid calculation and the radial symmetry-based algorithm.⁶ Since localization microscopy is a 2-D wide-field imaging modality using a CCD camera as the detector, acquiring a 3-D superresolution image needs sophisticated PSF engineering to represent the depth information in the distinct characteristics of the microscope’s PSFs, such as astigmatism,⁷ double helix,⁸ Airy function,⁹ saddle-point and tetrapod,^10,11 which require appropriate fitting algorithms to retrieve such information.

Another approach to obtain depth information from a 2-D image is based on defocused imaging. The PSF in a defocused image is made of a central spot and multiple concentric rings due to Fresnel diffraction.¹² The radii of the rings are correlated to the defocus level, thus the depth ($z$) information. By locating the centroid and measuring the radius of the outermost ring, the 3-D position of a single fluorophore can be determined. Such defocused imaging has been used to track single particles with Ångström accuracy.¹³ Recently, we developed a temporal focusing two-photon microscope that can track single fluorophores at nanometer precision with a depth range of $100\mu \mathrm{m}$.¹⁴ The temporal focusing two-photon microscope is a widefield 2-D imaging modality. It can achieve $z$-sectioning capability (axial resolution on the order of a few micrometers) by stretching the femtosecond laser pulse in the temporal domain before the objective lens and compressing it to its shortest temporal width at the focal plane of the objective.^15,16 When a fluorophore is out of the focal plane, it forms a defocused image on the CCD camera. We implemented the same calculation method to obtain its 3-D position at 50-nm precision from such defocused images.

There are two obstacles preventing the wide adoption of defocused imaging for 3-D localization. First, to our best knowledge, there is no well-developed mathematical model for a defocused PSF,¹⁷ which makes it difficult to develop localization algorithms based on fitting a known PSF. Second, the signal-to-noise ratio (SNR) in a defocused image is substantially lower than that in an in-focus image since the emitted light is spread out and detected by many pixels of the CCD camera instead of only a few. Furthermore, the excitation efficiency is usually lower than that in in-focus imaging. Therefore, developing localization algorithms will greatly improve the capability of defocusing imaging in 3-D localization of particles.

A practical consideration in developing algorithms is that diffraction rings can exhibit an elliptical shape due to optical aberration such as astigmatism. In this article, we developed a particle localization algorithm based on radial symmetry and ellipse fitting to localize the center of defocused PSF, as well as both horizontal and vertical radii of the ellipse. Our method can localize the 3-D position of a fluorophore particle within 20-nm precision in three dimensions with a range of $40\mu \mathrm{m}$ in $z$-dimension using a single defocused 2-D image.

2.

Materials and Methods

2.2.

Particle Localization

The flowchart of our algorithm is shown in Fig. 1. The approach is to first preprocess the raw image to retain the region of interest (outermost ring) and crop out other parts. Manual segmentation of a large number of images is neither accurate nor practical. Two subroutines are developed for automatic preprocessing as shown in the first box (rough center localization) and the second box (rough radius estimation) of Fig. 1. After the first two steps, ellipse fitting subroutine is executed to calculate the lengths of the long and short axes as shown in the third box.

Fig. 1

Flowchart of proposed particle localization algorithm.

2.2.1.

Rough center localization

To utilize the twofold rotational symmetry in an ellipse, the first step in rough center estimation is to determine the symmetrical axis in four directions by calculating the correlation between the original image and a new image obtained from the original image after 180-deg rotation. For a given raw image $I$, four equally spaced radial directions with polar angles at 0 deg, 45 deg, 90 deg, and 135 deg are chosen. For each direction, a series of parallel lines $l_{1,2,3,\dots ,n}$ are tested as the centrosymmetric axis by calculating symmetric correlation coefficients $c_i(i=1\dots n)$ as

Eq. (1)

$$c_i=\frac{{\sum }_{x=1}^H{\sum }_{y=1}^W[I(x,y)-\overline{I}]\times [I_r(x,y)-{\overline{I}}_r]}{\sqrt{{[{\sum }_{x=1}^H{\sum }_{y=1}^W(I(x,y)-\overline{I})]}^2\times {[{\sum }_{x=1}^H{\sum }_{y=1}^W(I_r(x,y)-{\overline{I}}_r)]}^2}},$$

where $I$ and $I_r$ denote pixel intensities in the original image and the new image after 180-deg rotation, $\overline{I}$ is the average intensity, and $H$ and $W$ are the height and width of images, respectively.

The variation of $c_i$ versus $i$ in all four directions is shown in Figs. 2(a)–2(d), respectively, and the line with maximum correlation coefficient is the symmetry axis in that direction [red points in Figs. 2(a)–2(d)]. The symmetry axis in each direction is represented as

Eq. (2)

$$y=K_{d}x+C_d,d=0,45,90,135.$$

Fig. 2

Rough center localization. (a)–(d) The distribution of correlation in four directions, 0 deg, 45 deg, 90 deg, and 135 deg, respectively. Horizontal axis is the index of parallel line. The red dot in each figure marks the peak of $c_i$, and the numbers in brackets show index $i$ value and the maximal correlation coefficient. (e) The original image with four blue lines indicating the symmetrical axes in four directions. (f) The magnified view of (e). The red circle indicates the estimated center.

Ideally, these lines intersect exactly at the center. However, due to noise, this does not happen, as shown in the zoomed-in image near the center [Figs. 2(e) and 2(f)]. We determine the center based on least square model, which searches for a pixel $(x,y)$ that has the minimum sum of Euclidean distances $D_{\mathrm{sum}}$ to all four lines by the following equation:

Eq. (3)

$$D_{\mathrm{sum}}(x,y)=\sum _d\frac{|K_{d}x-y+C_d|}{\sqrt{1+K_d^2}}.$$

Apparently, $D_{\mathrm{sum}}$ achieves minimum at $(x_s,y_s)$ when the following two criteria are met simultaneously

Eq. (4)

$$\{\begin{array}{c}\frac{\partial D_{\mathrm{sum}}(x_s,y_s)}{\partial x_s}=0\\ \frac{\partial D_{\mathrm{sum}}(x_s,y_s)}{\partial y_s}=0\end{array}.$$

Therefore, $(x_s,y_s)$ are the rough center coordinates [red circles in Figs. 2(e) and 2(f)].

2.2.2.

Rough radius estimation

Based on the center found above, a histogram-based radius estimation algorithm is applied. Figure 3(a) gives a schematic demonstration of calculating the distance-to-center histogram. For all the pixels that have equal Euclidean distance ($D_0$) to the center $(x_s,y_s)$ marked with a red cross, the average pixel intensity $I_0$ was calculated and is plotted as a function of $D_0$ in Fig. 3(b). If the noise level is low, this histogram curve should have a peak as the average pixel intensity reaches its maximum value when the Euclidean distance-to-center equals the radius. However, high-frequency noises appear in this histogram, and this peak is not identifiable. To remove these noises without greatly distorting the signal, the Savitzky–Golay (SG) filter is applied to smooth the histogram by a 597-point quadratic polynomial.¹⁸ Compared with other methods such as derivative calculation and wavelet transform (WT), the SG filter performs a fast least-square-fit convolution procedure by setting only two parameters: the width of the smoothing window and the degree of the smoothing polynomial, which significantly saves computation costs with fewer manual parameter settings.¹⁹ Rough radius $R_0$ is calculated by finding the local maximum of the smoothed curve [red dot in Fig. 3(b)].

Fig. 3

Rough radius estimation. (a) Schematic graph interprets the histogram. Given a distance ($D_0$) to center ($x_s$, $y_s$), the corresponding averaged intensity is calculated by the average intensities of all pixels whose Euclidean distances to center equal $D_0$. (b) Blue and magenta curves represent the original histogram and S–G filtering result. The horizontal coordinate of the red solid dot indicates the radius.

2.2.3.

Ellipse fitting

With the estimated center $(x_s,y_s)$ and radius $R_0$ known, most of the background region can be excluded from calculation except the region of interest (the ring), which was determined as

Eq. (5)

$${(x-x_s)}^2+{(y-y_s)}^2\in [{(R_0-b)}^2,{(R_0+b)}^2].$$

The parameter $b$ is the half width of the ring and is empirically set as five pixels in this experiment. The equation of an ellipse is represented by an implicit second-order polynomial

Eq. (6)

$$E(P;X)=P^T·X=a{x}^2+bxy+c{y}^2+dx+ey+f=0,$$

where $P={[a,b,c,d,e,f]}^T$ and $X={[x^2,xy,y^2,x,y,1]}^T$. In a 2-D image, each pixel coordinate is represented as $(x_k,y_k)$, and the pixel intensity is $I_{\mathrm{k}}$ ($k=1\dots N$, $N$ is the total number of data points). $E(P;X_k)$ represents the algebraic distance between point $X_k$ and the ellipse $E(P;X)=0$. According to least-square criterion, the optimal ellipse fitting can be achieved by minimizing the sum of squared algebraic distance of all $N$ data points to the ellipse [Eq. (7)], and the optimal coefficient vector $P_0$ is achieved when such error is minimal as shown in Eq. (8)

Eq. (7)

$$D_a(P)=\sum _{k=1}^{N}E{(P;X_k)}^2,$$

Eq. (8)

$$P_0={\mathrm{argmin}}_{\mathrm{P}}{D}_a(P).$$

After the background is cropped out, there are still noises inside the rings. Intuitively pixels with higher intensity are more likely to be the real signal and pixels with lower intensity tend to be noise. Therefore, a weight factor $w_k$ is implemented to enhance signal and suppress noise

Eq. (9)

$$w_k=\mathrm{norm}(I_k^c),$$

where $c$ is a positive parameter to control the weight among pixels with low (noise) and high (signal) intensities. We empirically set $c=10$. Now the squared algebraic distance with weight factor from the ellipse to $N$ data points is represented as

Eq. (10)

$$D_a(P)=\sum _{k=1}^N{w}_k*E{(P;X_k)}^2.$$

When finding the polynomial coefficient vector $P$ by minimizing $D_a(P)$, some constraint is needed to avoid the trivial solution $P=0_6$ and several possible solutions representing the same ellipse. In general, constraints can be expressed in the matrix form of $P^{T}CP=1$, where $C$ is a $6\times 6$ constraint matrix. Here, the Bookstein constraint, one type of Euclidean-invariant constraint, is implemented by the following matrix equation:²⁰

Eq. (11)

$$P^T\left[\begin{array}{cccc}1& 0& 0& \\ 0& 1/2& 0& 0_{3\times 3}\\ 0& 0& 1& \\ & 0_{3\times 3}& & 0_{3\times 3}\end{array}\right]P=1.$$

Now the minimization of Eq. (10) can be solved by considering a rank-deficient generalized eigenvalue system with a Lagrange multiplier $\lambda$

Eq. (12)

$$S^{T}SP=\lambda CP,$$

where $S={[X_1,X_2,\dots ,X_N]}^T$. $P$ is solved by calculating generalized eigenvectors of Eq. (12) under the constraint of Eq. (11).²¹ Figure 4 shows an example of the fitting with $P={[0.0036,0.0004,0.0035,-0.4796,0.4329,30.3565]}^T$.

Fig. 4

Weighted least square ellipse fitting. (a) Original image, (b) ellipse fitting based on ring cropping of (a). The red solid curve in (b) is the optimal solution.

3.

Results and Discussion

The performance of this algorithm is first tested on simulated noisy images, and the algorithmic error is evaluated with a method developed in Ref. 6. Simulated rings with intensity obeying Gaussian distribution are considered ideal noise-free images ($I_{\text{ideal}}$), and the long and short axes are $L$ and $S$ with the center at $(x_s,y_s)$. Simulated noisy images ($I_{\text{noisy}}$) with artificially added Poisson noises to $I_{\text{ideal}}$ spanning a range of SNR are used to evaluate this algorithm. Figure 5 shows three simulated noisy images with SNR at 1.05, 1.5, and 2, respectively. The algorithmic error is defined as the difference between algorithm outputs from $I_{\text{noisy}}$ and real ring parameters from $I_{\text{ideal}}$. For the purpose of brevity, only the algorithmic errors of the horizontal axis ($\mathrm{\Delta }R$) and center coordinate $(\mathrm{\Delta }{x}_s,\mathrm{\Delta }{y}_s)$ are displayed in the figures.

Fig. 5

Simulated noisy images with different levels of SNR. (a) $\mathrm{SNR}=1.05$, (b) $\mathrm{SNR}=1.5$, and (c) $\mathrm{SNR}=2$.

In the first test, different levels of Poisson noises ranging from $\mathrm{SNR}=1.01$ to 2 at intervals of 0.1 are added. At each SNR level, 1000 test images are generated. The mean of absolute value of algorithmic error (red circle dots) and standard deviation (blue plus signs) of algorithmic error at different SNR levels are shown in Fig. 6(a). Furthermore, since the SNR of experimental raw images is around 1.05, we calculate the mean of absolute value of algorithmic error (red circle dots) and standard deviation (blue plus signs) of error in a fine range from $\mathrm{SNR}=1.01$ to 1.1 with intervals of 0.01 [Fig. 6(b)]. Data of 1000 images at the SNR level of 1.05 is shown in Fig. 6(c). The errors $\mathrm{\Delta }R$, $\mathrm{\Delta }{x}_s$, and $\mathrm{\Delta }{y}_s$ in these 1000 experiments are randomly distributed with a standard deviation of 9.69, 7.62, and 7.51 nm with corresponding parameters in $I_{\text{ideal}}$ as 0. These tests show that, under the current image SNR of around 1.05, we have $\sim 95\%$ confidence to expect the algorithmic error for $R$, $x_s$, and $y_s$ to be around 20, 15, and 15 nm, respectively.

Fig. 6

Evaluation of algorithmic error and mechanical instability. (a) Mean and standard deviation of $\mathrm{\Delta }R$, $\mathrm{\Delta }{x}_s$, and $\mathrm{\Delta }{y}_s$ versus SNR in large SNR range (from 1.01 to 2 at 0.1 intervals). (b) Mean and standard deviation of $\mathrm{\Delta }R$, $\mathrm{\Delta }{x}_s$, and $\mathrm{\Delta }{y}_s$ versus SNR in short SNR range (from 1.01 to 1.1 at 0.01 intervals). Note that at each SNR level in (a) and (b), such noise adding is repeated 1000 times. (c) The distribution of $\mathrm{\Delta }R$, $\mathrm{\Delta }{x}_s$, and $\mathrm{\Delta }{y}_s$ versus number of experiments when SNR is set as 1.05. The mean SNR of experimental raw data is 1.05.

This algorithm is further applied on real microscopic images. Figures 7(a)–7(c) show three experimental defocused images acquired by the temporal focusing two-photon microscope at $z=20$, 35, and $50\mu \mathrm{m}$, respectively. At most depths, the center lobe still has the highest intensity, and the outmost ring intensity increases as the level of defocus increases. Figures 7(d) and 7(e) show the distribution of SNR for all raw images in group #1 and group #2, respectively. The mean of SNR is 1.045 with standard deviation 0.032 in group #1, and the mean of SNR is 1.051 with standard deviation 0.008 in group #2. The elliptical long and short radii for both experimental groups are shown in Figs. 8(a) and 8(b), respectively. In a long range, the relationship of radius versus $z$ is not linear; however, in a short range, it shows better linearity. The ratio of long axis over short axis is about 1.02, indicating the aspect ratio of ellipse. Apparently the fluctuations of the calculated radii of rings are much larger than the means and standard deviations in test images (Fig. 6). The most possible reason is the mechanical instability of the microscope since it is not equipped with feedback control.

Fig. 7

Experimental images at (a) $z=20\mu \mathrm{m}$, (b) $z=35\mu \mathrm{m}$, and (c) $z=50\mu \mathrm{m}$, (d) and (e) The SNR level is between 1 and 1.15.

Fig. 8

Radius versus depth for (a) group #1 (large range) and (b) group #2 (small range). Red: long axis and blue: short axis.

To estimate the mechanical instability, the same specimen was kept stationary and imaged for 20 min with a total of 20 images taken at 1-min intervals. The fluctuation of calculated horizontal radius over time after subtracting average is shown in Fig. 9. Since at $\mathrm{SNR}=1.05$ the algorithmic errors for $R$, $x_s$, and $y_s$ are around 20, 15, and 15 nm, respectively, [Fig. 6(c)], this shows that the microscope has a fluctuation for $R$, $x_s$, and $y_s$ in round 50 nm. The mean and standard deviation of Figs. 6(c) and 9 are summarized in Table 1. These tests demonstrate that mechanical instability is the limiting factor in demining particle position precision in our current experiments.

Fig. 9

Mechanical stability measurement for $R$, $x_s$, and $y_s$.

Table 1

Mean and standard deviation of ΔR, Δxs, and Δys for simulated images [Fig. 6(c)] and raw experimental images (Fig. 8).

All calculations were run on a Lenovo™ Ideapad 320 laptop computer with an AMD A12-9720P processor. The average execution time of this ring localization algorithm for a single image is shorter than 1.3 s (Table 2).

Table 2

Mean and standard deviation of computation time for rough center localization (S1), rough radius estimation (S2), and ellipse fitting (S3).

In summary, we developed an advanced algorithm to optimize the precision of localization in defocused imaging superresolution microscopy. A previous particle tracking work based on the calculation of radial symmetry centers from in-focus superresolution images achieved near theoretical limits with orders-of-magnitude faster execution time over the Gaussian fitting approach.⁶ With improvement of mechanical stability, we expect that this method could achieve nanometer-level resolution under current SNR. The advantage of defocused imaging is the capability of imaging multiple fluorophores at different depths simultaneously. Such capability allows us to achieve 3-D superresolution imaging without PSD engineering. The challenge is that, when multiple particles are imaged simultaneously, their PSFs overlap. Further development is needed for segmentation of overlapping rings in such scenarios.