AI Video Guide AI Picture GuideAI One SentenceAI Short Abstract

【AIGC One Sentence Reading】：Our upsampled PSF method enables high-accuracy 3D superresolution imaging with sparse sampling, enhancing imaging throughput and field of view in SMLM.AI Voice 

【AIGC Short Abstract】：We propose an upsampled PSF inverse modeling method for SMLM, allowing high-accuracy 3D superresolution imaging with sparse sampling. This approach enhances imaging throughput, reduces data volume, and expands the field of view, all while maintaining precise single-molecule localization.AI Voice 

Note: This section is automatically generated by AI . The website and platform operators shall not be liable for any commercial or legal consequences arising from your use of AI generated content on this website. Please be aware of this.

1. INTRODUCTION

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

Figure 1.Concept of upsampled PSF modeling. First, the z-stack data of beads are collected (top left), from which ROIs are extracted to construct an experimental PSF library at different positions (top center). Next, the initial global parameters (upsampled PSF represented by a 3D matrix) and local parameters (3D position, photon count, and background) are used to generate a fine PSF library. Then, the fine PSF library is pixel-integrated to create a coarse PSF library (top right). The loss function, based on MLE, is calculated between the coarse PSF library and the experimental PSF library (top center). Finally, through backpropagation (L-BFGS-B optimization algorithm), both global and local parameters are updated iteratively, ultimately yielding the final upsampled PSF. We also localize the experimental bead data with both the upsampled PSF (bottom center) and averaged coarse PSF, revealing that only the upsampled PSF could achieve unbiased and optimal localization results (bottom left).

2. RESULTS

A. Upsampled PSF Inverse Modeling

B. Localization Algorithm for Large-Pixel Data

C. Investigating Optimal Pixel Size for Large-Pixel Imaging

Figure 2.Comparison of the CRLB and RMSE for PSFs with different pixel sizes. (a)–(c) Comparison of the theoretical CRLB for $x$, $y$, and $z$ for PSFs with different pixel sizes. (d)–(f) Comparison of the RMSE in the $x$, $y$, and $z$ directions localized with PSFs of different pixel sizes. $R_{up}$ represents the percentage degradation in CRLB or RMSE averaged over the range from $-600$ to 600 nm, relative to the values obtained with a pixel size of 110 nm. Each molecule contributes a total of 2000 photons, with a background level of $2\times {10}^{-3}\text{photons}/{\mathrm{nm}}^2$. Localization accuracy at each $z$ position is computed from RMSE of the localization results from 1000 single molecules.

D. Validation Using Simulated Data

Figure 3.Validation of the upsampled PSF modeling using simulated data from astigmatic PSF. (a) The 330 nm pixel size PSF (top) and 110 nm pixel size astigmatic PSF (bottom) estimated from 300 simulated data sets with a pixel size of 330 nm. (b)–(d) Comparison of the impact of the estimated 330 and 110 nm pixel size astigmatic PSF on the localization accuracy of the 330 nm pixel size simulated data along the $x$, $y$, and $z$ directions, respectively. Each molecule contributes a total of 10,000 photons, with a background level of $2\times {10}^{-3}\text{ photons}/{\mathrm{nm}}^2$. The $x/y/z$ bias as a function of $z$ is defined as in Appendix F Eqs. (F1) and (F2). Solid lines and shaded areas in localization plots indicate mean and standard deviation of bias, respectively.

Figure 4.Nup96-AF647 in U2OS cells reconstructed using the upsampled PSF model. (a) Single-molecule data of Nup96-AF647 with a pixel size of 127 nm were $2\times 2$ binned to produce undersampled data with a pixel size of 254 nm. Localization of the undersampled data using the upsampled PSF yielded a superresolved image. The red boxes indicate the single-molecule image with different pixel sizes. In the superresolution image, the white dashed box in the lower left corner corresponds to a magnified view of the region marked by the solid box. (b) $XZ$ view of the selected region (white dashed lines) in (a) reconstructed from single-molecule data using PSFs of different pixel sizes. (c) Superresolution image of 321 nm pixel size single-molecule data of Nup96-AF647 analyzed with 321 and 107 nm pixel size PSFs, respectively. The $XZ$ views of the selected region [white dashed lines in the left image of (c)] were compared.

E. Validation Using Experimental Data

3. DISCUSSION

APPENDIX A: METHODS AND EXPERIMENTAL DETAILS

SMLM imaging was performed at room temperature using a custom-built microscope setup (Fig. 11 in Appendix G). The sample was illuminated via a laser box coupled to a single-mode fiber. The excitation laser was reflected by a dichroic mirror (Di01-R405/488/561/635, Semrock) and focused on the back focal plane of a high NA objective (NA 1.5, UPLAPO 100XOHR, Olympus) for superresolution imaging. Fluorescence emission was collected by the objective and filtered through a quad-band emission filter (NF03-405/488/561/635E, Semrock). Following the tube lens (SWTLU-C, $f=180\mathrm{mm}$, Olympus), a bandpass filter (ET685/70m, Chroma) was inserted to remove residual laser light. Subsequently, the fluorescence passed through a 4f system composed of lenses L3 (AC254-150-A-ML, $f=150\mathrm{mm}$, Thorlabs) and L4 (AC254-030-A-ML, $f=30\mathrm{mm}$, Thorlabs). Finally, images were acquired using an sCMOS camera (ORCA-Flash4.0 V3, Hamamatsu) with a pixel size of $6.5\mathrm{\mu m}\times 6.5\mathrm{\mu m}$. Additionally, a plano–convex cylindrical lens (LJ1516L1-C, Thorlabs) was positioned in front of the camera to introduce astigmatism for 3D imaging.

Cell Culture. U2OS cells (Nup96-SNAP, catalog No. 300444, Cell Line Services) were cultured in DMEM supplemented with 10% (volume fraction) fetal bovine serum (FBS), $1\times$ penicillin-streptomycin (PS), and $1\times$ MEM nonessential amino acids (NEAAs; catalog No. 11140-050, Gibco). The cells were maintained at 37°C in a humidified incubator with 5% ${\mathrm{CO}}_2$ and were passaged every 2–3 days. Prior to seeding, 25 mm high-precision round glass coverslips (No. 1.5H, catalog No. CG15XH, Thorlabs) were thoroughly cleaned by sequential sonication in 1 mol/L potassium hydroxide (KOH), Milli-Q water, and ethanol, followed by 30 min of UV sterilization. For superresolution imaging, the U2OS cells were seeded onto the cleaned coverslips and cultured for 2 days until reaching 80%–90% confluency. Routine mycoplasma tests confirmed the absence of contamination throughout the study.

Fluorescence bead sample. 100 nm TetraSpeck beads (Thermo Fisher Scientific, catalog No.  T7279) were used. A clean 25 mm coverslip was incubated with 40 μL of 100 mmol/L $\mathrm{Mg}{\mathrm{Cl}}_2$ and 360 μL of 1:1000 diluted bead solution for 5 min. The coverslip was then thoroughly rinsed with Milli-Q water 3 times. After washing, the coverslip was placed in a custom sample holder, and 1 mL of Milli-Q water was added to maintain hydration.

Nup96-SNAP-AF647 in U2OS cells. To label Nup96, U2OS-Nup96-SNAP cells were prepared following established protocols. Cells were initially prefixed in 2.4% paraformaldehyde (PFA) for 30 s, permeabilized with 0.4% Triton X-100 for 3 min, and then fixed in 2.4% PFA for an additional 30 min. Following fixation, the cells were quenched with 0.1 mol/L ${\mathrm{NH}}_4\mathrm{Cl}$ for 5 min and washed twice with PBS. To minimize nonspecific binding, cells were blocked for 30 min using Image-iT FX Signal Enhancer (Invitrogen, catalog No. I36933). For labeling, cells were incubated for 2 h in a dye solution containing 1 μmol/L SNAP-tag ligand BG-AF647 (New England Biolabs, catalog No. S9136S), 1 mmol/L dithiothreitol (neoFroxx, catalog No. 1111GR005), and 0.5% bovine serum albumin (BSA) in PBS. After incubation, cells were washed 3 times with PBS for 5 min each to remove unbound dye. Finally, cells were postfixed with 4% PFA for 10 min, washed 3 times with PBS, and stored at 4°C until imaging.

Collection of bead data on single-objective SMLM systems. The fluorescence bead sample was prepared as described above. When $\mathrm{bin}=2$, we typically consider that collecting data from more than 10 bead stacks is sufficient for robust upsampled PSF estimation; when $\mathrm{bin}=3$, collecting data from over 30 bead stacks is considered adequate. Each bead stack was acquired by moving the sample stage from $-1$ to 1 μm, with a step size of 50 nm. One frame per $z$ position was collected, and each frame could include multiple beads, with approximately 16,000 photons and 80 background photons per bead. If the photon count per bead is low, the number of beads can be increased accordingly. Our method, based on Poisson distribution estimation, is adaptable to data with different photon counts (SNR). When using the demo on our GitHub to fit the upsampled PSF, a localization bias average below 5 nm indicates reliable results.

Imaging of Nup96-AF647 on single-objective SMLM systems. Nup96-SNAP-AF647 labeled U2OS cells were imaged with the Hamamatsu sCMOS camera over an FOV covering $56\mathrm{pixels}\times 56\mathrm{pixels}$. The camera was operated under the rolling shutter readout mode with an exposure time of 20 ms. 150,000 frames were acquired. The position of the cylindrical lens before camera was adjusted so that $\sim 80\mathrm{nm}$ astigmatism aberration was introduced to the system.

APPENDIX B: UPSAMPLED PSF MODELING IN THE SPATIAL DOMAIN

The initial position for each bead is $$x_{\mathrm{init}}=y_{\mathrm{init}}=z_{\mathrm{init}}=0.$$(B1)

The background value for each bead is estimated as $${\mathrm{bg}}_i=\min (D_i(x,y,z)\otimes G(x,y,z,{\sigma }_x={\sigma }_y={\sigma }_z=2))/{\mathrm{bin}}^2,$$(B2)where $D_i(x,y,z)$ represents the cropped 3D image stack of bead $i$, while $G$ denotes a 3D Gaussian kernel with standard deviations ${\sigma }_x$, ${\sigma }_y$, and ${\sigma }_z$ in the respective $x$, $y$, and $z$ directions. bin factor refers to the ratio between the pixel size of the data and the pixel size of the upsampled PSF. The symbol $\otimes$ denotes the convolution operation. The minimum value is taken over all voxels within the cropped and filtered region. The initial background for each bead is set to the median value of ${\mathrm{bg}}_i$ across all bead data, $${\mathrm{bg}}_{\mathrm{init}}=\underset{i}{\text{median}}(b_i).$$(B3)

The initial photon value for each bead is estimated as $${\mathrm{Np}}_{\mathrm{init}}^i=\underset{z}{\mathrm{avg}}(\sum _{x,y}{D}_i(x,y,z)-b_i).$$(B4)

The initial PSF model is a 3D array with each element set to $1/(N_s^2\cdot {\mathrm{bin}}^2)$ with $N_s$ the ROI size of bead data. For example, if the ROI size is $N_s\times N_s=21\times 21$ pixels and the sampling factor bin is 4, the initial value is 0.00014. This assumes that the sum of each axial slice of the PSF model is 1.

The upsampled PSF model at a given bead location ($x_i$, $y_i$, $z_i$) can be calculated as follows [22]: $$U_i^{\mathrm{up}}(x-x_i,y-y_i,z-z_i)={\mathcal{F}}_{3D}^{-1}({\mathcal{F}}_{3D}({\mathrm{Np}}_i·{\mathrm{PSF}}_{\mathrm{up}}(x,y,z)+{\mathrm{bg}}_i)·g(q)·e^{i{\phi }_{\text{shift}}}),$$(B5)where the shifting phase is calculated from $${\phi }_{\text{shift}}=2\pi (q_x{x}_i+q_y{y}_i+q_z{z}_i),$$(B6)where $q_x$, $q_y$, $q_z$ are the Cartesian coordinates of the frequency space as given by $$q_x=\frac{x}{L_x},\hspace{1em}{q}_y=\frac{y}{L_y},\hspace{1em}{q}_z=\frac{z}{L_z},$$(B7)where $L_x$, $L_y$, and $L_z$ are length of the upsampled PSF model in $x$, $y$, and $z$ directions. Considering the bead size, we model the bead as a sphere and $g(q)$ is the analytical Fourier transform of a solid sphere of radius $r_0$, $$g(q)=\frac{J_{\frac{3}{2}}(2\pi q{r}_0)}{{(q{r}_0)}^{\frac{3}{2}}}{r}_0^3,$$(B8)where $J$ is the Bessel function of the first kind and $q$ is the spherical coordinate in the frequency space and is calculated from $$q=\sqrt{q_x^2+q_y^2+q_z^2}.$$(B9)

To ensure that the variation of $U_i^{\mathrm{up}}$ is not affected, $g(q)$ is normalized to its maximum value. To match $U_i^{\mathrm{up}}$ with the pixel size of the data, we further apply a binning operation to $U_i^{\mathrm{up}}$, $$U_i(k,l)=\sum _{i=1}^{\mathrm{bin}}\sum _{j=1}^{\mathrm{bin}}{U}_i^{\mathrm{up}}(\mathrm{bin}·k+i,\mathrm{bin}·l+j),$$(B10)where $U_i$ is the binned (merged) forward model, which can be directly compared with the measured bead data $M_i$. $k$ and $l$ denote the horizontal and vertical pixel indices of the forward model.

We utilized the L-BFGS-B algorithm from the SciPy optimization package for our optimization process. To account for the varying scales of gradients across different variable types, we applied a tailored scaling factor to each type. This adjustment ensured uniform update rates for all variables throughout each iteration. The original variables in the forward model were substituted with transformed variables, ensuring consistent and balanced updates during the optimization procedure. For example, $$\begin{gathered} {\mathrm{Np}}_i=\frac{{\mathrm{Np}}_i}{w_{{\mathrm{Np}}_i}}{w}_{{\mathrm{Np}}_i}={\mathrm{Np}}_{w,i}{w}_{\mathrm{Np}}, \\ {\mathrm{bg}}_i=\frac{{\mathrm{bg}}_i}{w_{{\mathrm{bg}}_i}}{w}_{{\mathrm{bg}}_i}={\mathrm{bg}}_{w,i}{w}_{\mathrm{bg}}, \\ {\mathrm{PSF}}_{\mathrm{up}}=\frac{{\mathrm{PSF}}_{\mathrm{up}}}{w_{{\mathrm{PSF}}_{\mathrm{up}}}{\mathrm{bin}}^2}{w}_{{\mathrm{PSF}}_{\mathrm{up}}}={\mathrm{PSF}}_w{w}_{\mathrm{PSF}}, \end{gathered}$$(B11)where ${\mathrm{PSF}}_w$, ${\mathrm{Np}}_{w,i}$, ${\mathrm{bg}}_{w,i}$ are scaled variables, with respect to which the gradient will be calculated, and $w_{\mathrm{PSF}}$, $w_{\mathrm{Np}}$, $w_{\mathrm{bg}}$ are scaling factors for the upsampled PSF model, photons, and background. The scaling factors can be determined based on the gradient values of their respective parameters, ensuring that all gradient values are within the same order of magnitude. Proper variable scaling is crucial in PSF learning, as it ensures that the optimization process converges to the global minimum and facilitates the efficient optimization of all variables.

The loss function for upsampled PSF estimation is $$\begin{gathered} \mathrm{loss}=\mathrm{LL}+a_{\text{drift}}{f}_{\text{drift}}+\beta (a_{\mathrm{PSF},\min }{f}_{\mathrm{PSF},\min }+a_{\mathrm{bg},\min }{f}_{\mathrm{bg},\min } \\ +a_{\mathrm{Np},\min }{f}_{\mathrm{Np},\min }+a_{\mathrm{norm}}{f}_{\mathrm{norm}}), \end{gathered}$$(B12)where $\mathrm{LL}$ is the log-likelihood function, assuming that the converted pixel values follow a Poisson distribution, $$\mathrm{LL}=\underset{i}{\mathrm{avg}}(U_i-D_i-D_i\log (U_i)+D_i\log (D_i)).$$(B13)

Note that $D_i$ and $D_i\log (D_i)$ are normalization factors for the likelihood probabilities, aimed at improving the stability of the fitting process.

When sample drift is considered, $f_{\text{drift}}$ is equal to the $L^1$ norm of all drift rates over the bead data; this term is to constrain the estimated drift rates to be close to zero, to avoid adding an arbitrary constant drift rate. We define the bead’s lateral position at each axial slice as $$\begin{gathered} x_{i,z}=g_{x,i}z,\hspace{1em} \\ {y}_{i,z}=g_{y,i}z, \end{gathered}$$(B14)where $z=1,2,\dots ,N_z$, and $g_{x,i}$ and $g_{y,i}$ are additional learning variables that define the drift rates along $x$ and $y$ for each bead stack. To apply the $z$-dependent lateral drift, the obtained forward model is shifted slice-wise through a 2D Fourier transform, $$\begin{gathered} U_{\text{drift},i}(x-x_{i,z},y-y_{i,z},z) \\ ={\mathcal{F}}_{2D}^{-1}({\mathcal{F}}_{2D}(U_i(x,y,z))e^{i2\pi (k_x{x}_{i,z}+k_y{y}_{i,z})}). \end{gathered}$$(B15)

The next three terms serve to constrain the values of the upsampled PSF model, photons, and background to be positive, $$\begin{gathered} f_{\mathrm{PSF},\min }=\sum _{x,y,z}\min {({\mathrm{PSF}}_{\mathrm{up}},0)}^2, \\ {f}_{\mathrm{Np},\min }=\sum _i\min {({\mathrm{Np}}_i,0)}^2, \\ {f}_{\mathrm{bg},\min }=\sum _i\min {({\mathrm{bg}}_i,0)}^2. \end{gathered}$$(B16)

Here min denotes an element-wise minimum, comparing each value in an array with zero. As default we set $a_{\mathrm{PSF},\min }=100$, $a_{\mathrm{Np},\min }=1$, $a_{\mathrm{bg},\min }=1$.

The final term, $f_{\mathrm{norm}}$, in the loss function is introduced to enforce the sum of each axial slice of the PSF model to remain constant. This constraint is based on the principle of energy conservation, ensuring that the total photon count of a single emitter remains invariant across different axial positions. Thus, we define $$f_{\mathrm{norm}}=\underset{z}{\mathrm{avg}}{\left(\sum _{x,y}U-\sum _{x,y,z}\frac{U}{N_z}\right)}^2.$$(B17)

This term is often utilized in conjunction with the estimation of $z$-dependent photon fluctuations to account for photobleaching and variations in illumination intensity during data acquisition. By default, we set $a_{\mathrm{norm}}=0$.

APPENDIX C: LOCALIZATION METHODS USING UPSAMPLED PSF

Before performing cubic spline interpolation on the estimated upsampled PSF, a convolution operation must first be applied. This is because the experimental data are obtained through integration over large pixel sizes, while the estimated upsampled PSF is based on integration over small pixel sizes. The difference between the two is determined by the bin factor, which represents the ratio between the large and small pixel sizes. To align the two, a convolution is performed on the upsampled PSF using a kernel of size corresponding to the bin factor, with all elements set to 1, and a stride of 1. The detailed convolution procedure is as follows:, $$\begin{gathered} {\mathrm{PSF}}_{\mathrm{con}}(k,l)={\mathrm{PSF}}_{\mathrm{raw}}\otimes \mathrm{ones}(\mathrm{bin},\mathrm{bin}) \\ =\sum _{i=1}^{\mathrm{bin}}\sum _{j=1}^{\mathrm{bin}}{\mathrm{PSF}}_{\mathrm{raw}}(sk+i,sl+j), \end{gathered}$$(C1)where ${\mathrm{PSF}}_{\mathrm{con}}$ represents the convolved upsampled PSF, $\otimes$ denotes the convolution operation, and $\mathrm{ones}(\mathrm{bin},\mathrm{bin})$ is a matrix of size equal to the bin factor, with all elements set to 1. $s$ denotes the convolution stride, fixed at 1 in this instance. The bin factor refers to the ratio between the pixel size of the data and the pixel size of the upsampled PSF.

For the localization step using the ${\mathrm{PSF}}_{\mathrm{con}}$, the data were analyzed with the cubic spline-fitting method. The cubic spline interpolation of a given PSF model is [19,20] $$f_{i,j,k}(x,y,z)=\sum _{m=0}^3\sum _{n=0}^3\sum _{p=0}^3{a}_{i,j,k,m,n,p}{\left(\frac{x-x_i}{\frac{\mathrm{\Delta }x}{\mathrm{bin}}}\right)}^m{\left(\frac{y-y_j}{\frac{\mathrm{\Delta }y}{\mathrm{bin}}}\right)}^n{\left(\frac{z-z_k}{\mathrm{\Delta }z}\right)}^p,$$(C2)where $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ are the $x$ and $y$ pixel sizes of the camera, $\mathrm{\Delta }z$ is the axial step size of the z-stack, $a_{i,j,k,m,n,p}$ are the spline coefficients, and $x_i$, $y_j$, $z_k$ are the start position of each voxel $(i,j,k)$ in the upsampled PSF. $x$, $y$ are the positions corresponding to the camera pixels, while $z$ represents the position in the z-stack. After building the spline PSF model, MLE with Poisson statistics was used to localize beads or single molecules with the objective function given by [21] $${\chi }_{\mathrm{mle}}^2=2(\sum _{k,j}(U_{k,j}-M_{k,j})-\sum _{k,j,M_{k,j}>0}{M}_{k,j}\ln \frac{U_{k,j}}{M_{k,j}}),$$(C3)where $U_{k,j}$ and $M_{k,j}$ are the expected photon number and measured photon number in the pixel $(k,j)$, respectively. We used a modified Levenberg–Marquardt (L-M) algorithm [28] to minimize ${\chi }_{\mathrm{mle}}^2$ for the parameter estimation.

APPENDIX D: VECTORIAL PSF MODEL CALCULATION

To accurately model the image formation process in a microscope equipped with a high NA objective, a vectorial PSF model is employed, accounting for the refractive index mismatch between the medium-cover slip interface and the cover slip-immersion medium interface. Given that fluorescent probes are typically flexibly attached to the target molecules and are capable of free rotation, we assume an isotropic emitter PSF model for our analysis. The vectorial PSF can be expressed as [25,26,29,30] $$\begin{gathered} \mathrm{PSF}(x-x_i,y-y_i,z-z_i)=\sum _{\begin{array}{c}m=x,yn=p_x,p_y,p_z\end{array}}|{\mathcal{F}}_{\mathrm{czt}}(h(k_x,k_y)w_{mn}{e}^{i2\pi (k_x(x-x_i)+k_y(y-y_i)-k_z(z-z_i))}){|}^2, \\ {k}_z=\sqrt{k^2-k_x^2-k_y^2}, \end{gathered}$$(D1)where $h(k_x,k_y)$ is the pupil function and can be written as $$h(k_x,k_y)=T_{a}A(k_x,k_y)e^{i\mathrm{\Phi }(k_x,k_y)},$$(D2)where $A$ and $\mathrm{\Phi }$ are the magnitude and phase components of the pupil function, and each of them is a 2D array. $T_a$ is the apodization factor of the objective lens and is equal to $$T_a=\frac{\sqrt{\cos {\theta }_{\mathrm{imm}}}}{\cos {\theta }_{\mathrm{med}}},$$(D3)where ${\theta }_{\mathrm{med}}$ and ${\theta }_{\mathrm{imm}}$ represent the angles of the optical rays in the sample medium and the immersion medium, respectively, and are constrained by the NA of the objective lens. The wave vector $k$ has Cartesian components $k_x$, $k_y$, $k_z$, with its magnitude $k=\frac{n_{\mathrm{imm}}}{\lambda }$, where $n_{\mathrm{imm}}$ is the refractive index of the immersion medium, and $\lambda$ is the central wavelength corresponding to the emission filter. The term $e^{-i2\pi k_z(z-z_i)}$ accounts for the defocus phase in the propagation.

We employed the chirp Z-transform based on Bluestein’s algorithm (denoted as ${\mathcal{F}}_{\mathrm{czt}}$) to compute the 2D Fourier transform of the pupil function. The advantage of using Bluestein’s algorithm is that the pupil function size becomes independent of the camera’s pixel size, allowing a pupil size of $64\mathrm{pixels}\times 64\mathrm{pixels}$ to offer sufficient sampling for accurate representation.

In this context, $w_{mn}$ denotes the m-component of the electric field at the pupil plane generated by the $n$-component of the dipole moment of the fluorophore. The calculation of $w_{mn}$ is as follows: $$\begin{gathered} w_{xn}=P_n\cos \phi -S_n\sin \phi , \\ \hspace{1em}{w}_{yn}=P_n\sin \phi -S_n\cos \phi , \end{gathered}$$(D4)where $\phi$ is the angular component in the polar coordinate of the frequency space. $P_n$ and $S_n$ are electric field components in p and s polarizations relative to the incident plane at the sample space, $$\begin{gathered} P_{p_x}=T_p\cos {\theta }_1\cos \phi , \\ {P}_{p_y}=T_p\cos {\theta }_1\sin \phi , \\ {P}_{p_z}=-T_p\sin {\theta }_1, \\ {S}_{p_x}=-T_s\sin \phi , \\ {S}_{p_y}=T_s\cos \phi , \\ {S}_{p_z}=0, \end{gathered}$$(D5)where $p_x$, $p_y$, $p_z$ are the Cartesian components of the dipole moments, and $T_p$ and $T_s$ are the total transmission coefficients of $\mathrm{p}$- and $\mathrm{s}$-polarized light. The fluorescence light starting from the dipole emitter propagates through the sample medium, the coverslip, and the immersion medium. In this three-layer system, we ignore the multiple reflections at the two interfaces; then $T_p$ and $T_s$ can be calculated from $$\begin{gathered} T_p={\tau }_{P13}{\tau }_{P23}, \\ {T}_s={\tau }_{S13}{\tau }_{S23}, \end{gathered}$$(D6)where ${\tau }_{Pij}$ and ${\tau }_{Sij}$ are the Fresnel transmission coefficients of $\mathrm{s}$- and $\mathrm{p}$-polarized light traveling from medium $i$ to medium $j$, $$\begin{gathered} {\tau }_{Pij}=\frac{2n_i\cos {\theta }_i}{n_i\cos {\theta }_j+n_j\cos {\theta }_i}, \\ {\tau }_{Sij}=\frac{2n_i\cos {\theta }_i}{n_i\cos {\theta }_i+n_j\cos {\theta }_j}, \end{gathered}$$(D7)where $n_i$ and ${\theta }_i$ are the refractive index and the light propagation angle in medium $i$, and the subscripts 1, 2, 3 denote the sample medium, the coverslip, and the immersion medium, respectively.

In order to realistically account for the effects of fluorescence dipole emission, pixelation, particle size, dispersion, and chromatic aberration, the PSF model is rescaled using the OTF as $${\mathrm{PSF}}_{\mathrm{otf}}={\mathcal{F}}_{3D}^{-1}({\mathcal{F}}_{3D}(\mathrm{PSF})·e^{-i2{({\sigma }_x{q}_x\pi )}^2-i2{({\sigma }_y{q}_y\pi )}^2}),$$(D8)where ${\sigma }_x$ and ${\sigma }_y$ are the standard deviation (in pixel unit) of a 2D Gaussian kernel in real space. By including the photon count and background of each bead stack, the final forward model is $$U=\mathrm{Np}·{\mathrm{PSF}}_{\mathrm{otf}}+\mathrm{bg},$$(D9)where $\mathrm{Np}$ and $\mathrm{bg}$ represent the photon count and background, respectively.

APPENDIX E: ANALYTICAL CRLB CALCULATION FOR PIXELATED PSF

The CRLB is calculated from the diagonal elements of the inverse of the Fisher information matrix $I(\mathrm{\Theta })$, which measures the amount of information that an observation (PSF) carries about the estimated parameters $\mathrm{\Theta }$, $$\mathrm{CRLB}{(\mathrm{\Theta })}_{ij}=\mathrm{diag}(I{(\mathrm{\Theta })}_{ij}^{-1}),$$(E1)where diag refers to extracting the diagonal elements of a matrix. The Fisher information matrix is defined as $$\begin{gathered} I{(\mathrm{\Theta })}_{ij}=\sum _j\sum _k\frac{1}{U_k}\frac{\partial U_{k,l}}{\partial {\mathrm{\Theta }}_i}\frac{\partial U_{k,l}}{\partial {\mathrm{\Theta }}_j}, \\ {U}_{k,l}={(\mathrm{Np}·{\mathrm{PSF}}_{\text{pixelated}}(k,l)+\mathrm{bg})}_k, \end{gathered}$$(E2)where $\mathrm{\Theta }$ is a set of parameters being estimated, and $k$ and $l$ denote the horizontal and vertical pixel indices of the pixelated PSF, respectively. The parameters $\mathrm{\Theta }$ for pixelated PSF are $x$, $y$, $z$, photon $\mathrm{Np}$, and background $\mathrm{bg}$ for each emitter. The $U_{k,l}$ is the forward model as defined in Eq. (B10).

The pixelated PSF and its first-order derivatives with respect to the parameters are obtained by integrating the upsampled PSF and its corresponding derivatives. The first-order derivatives of the pixelated forward model $U_k$ with respect to the parameters $\mathrm{\Theta }$ are given by $$\begin{gathered} {\mathrm{PSF}}_{\text{pixelated}}(k,l)=\sum _{i=1}^{\mathrm{bin}}\sum _{j=1}^{\mathrm{bin}}{\mathrm{PSF}}_{\text{upsampled}}(\mathrm{bin}·k+i,\mathrm{bin}·l+j), \\ \frac{\partial U_{k,j}}{\partial {\mathrm{\Theta }}_{x,y,z}}=\mathrm{Np}\sum _{i=1}^{\mathrm{bin}}\sum _{j=1}^{\mathrm{bin}}\sum _{\begin{array}{c}m=x,yn=p_x,p_y,p_z\end{array}}\mathrm{Re}\left(E_{mn}^{*}\frac{\partial E_{mn}}{\partial {\mathrm{\Theta }}_{x,y,z}}\right), \\ \frac{\partial U_k}{\partial \mathrm{Np}}={\mathrm{PSF}}_{\text{pixelated}}(k,l), \\ \frac{\partial U_k}{\partial \mathrm{bg}}=1, \\ {E}_{\mathit{mn}}={\mathcal{F}}_{\mathrm{czt}}(h(k_x,k_y)w_{\mathit{mn}}{e}^{i2\pi (k_{z\mathrm{med}}{z}_i-k_z{z}_s)}{e}^{i2\pi (k_x{x}_i+k_y{y}_i)}), \\ \frac{\partial E_{\mathit{mn}}}{\partial {\mathrm{\Theta }}_{x,y,z}}={\mathcal{F}}_{\mathrm{czt}}(i{k}_{x,y,z\mathrm{med}}h(k_x,k_y)w_{\mathit{mn}}{e}^{i2\pi (k_{z\mathrm{med}}{z}_i-k_z{z}_s)}{e}^{i2\pi (k_x{x}_i+k_y{y}_i)}), \end{gathered}$$(E3)where bin factor refers to the ratio between the pixel size of the data and the pixel size of the upsampled PSF.

Due to the significant variation in CRLB at different positions within the large-pixel vector PSF (Fig. 8 in Appendix G), 121 points were uniformly selected within a single pixel. The CRLB was calculated for each point, and the average of these values was taken to obtain the final mean CRLB, $${\mathrm{CRLB}}_{\mathrm{mean}}^{\frac{1}{2}}=\sum _{i=1}^{\mathrm{num}}{\mathrm{CRLB}}^{\frac{1}{2}}(x_i,y_i)/\mathrm{num},$$(E4)where $x_i$ and $y_i$ represent different position coordinates within a single pixel. The num represents the number of different positions sampled within a single pixel.

The CRLB calculated in the current work is slightly better than that in previous work using a numerical approach, especially for larger pixel conditions. We found that blurring of the theoretical PSF with a Gaussian function [Eq. (D8)] to mimic experimental PSF will reduce the impact of the pixel size on the theoretical CRLB, which was not considered by Huang et al. (Fig. 12 in Appendix G). We also showed that our upsampled spline PSF fitter could achieve a localization accuracy approaching the calculated CRLB (Fig. 8), showing that both the upsampled PSF modeling and upsampled spline PSF fitter could achieve the optimal resolution.

The calculation of ${\mathrm{CRLB}}_{3D}$ is as follows [26]: $${\mathrm{CRLB}}_{3D}=\frac{1}{N_z}\sum _{i=1}^{N_z}({\mathrm{CRLB}}_{\mathrm{mean},x}+{\mathrm{CRLB}}_{\mathrm{mean},y}+{\mathrm{CRLB}}_{\mathrm{mean},z}).$$(E5)

Similarly, the calculation of ${\mathrm{RMSE}}_{3D}$ is as follows: $${\mathrm{RMSE}}_{3D}=\frac{1}{N_z}\sum _{i=1}^{N_z}{({\mathrm{RMSE}}_x^2+{\mathrm{RMSE}}_y^2+{\mathrm{RMSE}}_z^2)}^{\frac{1}{2}},$$(E6)where $N_z$ represents the number of z-positions within the axial range.

APPENDIX F: LOCALIZATION TEST

To evaluate the accuracy of the learned upsampled PSF model, we performed localization tests using the same bead data that were utilized during training. The accuracy was assessed by measuring the localization bias in the $x$, $y$, and $z$ dimensions, with the localization bias in $z$ calculated from $$z_{\mathrm{bias},i,z}=z_{i,z}-z_{\mathrm{GT},i,z},$$(F1)where $Z_{\mathrm{GT}}$ represents the ground-truth position, which corresponds to the stage position. Similarly, the localization bias in $x/y$ is calculated from $$\begin{gathered} x_{\mathrm{bias},i,z}=x_{i,z}-\underset{z}{\text{median}}{x}_{i,z}, \\ {y}_{\mathrm{bias},i,z}=y_{i,z}-\underset{z}{\text{median}}{y}_{i,z}, \end{gathered}$$(F2)where $x_{i,z}$ and $y_{i,z}$ are estimated using MLE for the lateral localization of each 2D bead image, and the subscripts $i$ and $z$ refer to the indices of the beads and the axial slices within each bead stack, respectively.

APPENDIX G: SUPPLEMENTARY FIGURES

Figs. 5–12 are figures supplementary to the main text. More details are shown in their figure captions.

Figure 5.Comparison of the shape and intensity distribution of 330 nm pixelated PSFs at different positions. (a) PSFs with a pixel size of 330 nm at different $x$ positions showing the shape of PSF is spatially variant. $x_c$ is the $x$ position of the PSF, with 0 nm locating at the center of the middle pixel. (b) Relative intensity distribution at the white dashed line in (a) for PSFs at different $x$ positions.

Figure 6.Validation of spline fitting for different sampling rates. (a) To match the pixel size of the upsampled PSF with the simulated or experimental data, a convolution operation is performed on the estimated raw upsampled PSF (Raw) to generate the pixel integrated upsampled PSF (Convolved), using a convolution kernel size of $\mathrm{bin}\times \mathrm{bin}$ with all values set to 1 [Eq. (B15)]. The convolved PSF is used for localizing experimental data. The bin factor represents the ratio between the pixel size of the simulated or experimental data and the pixel size of the upsampled PSF, with its value constrained to integer values. (b) The spline fitter returns the $z$-position results by fitting the same simulated data using upsampled PSFs at different sampling rates. (c) Localization results of 110, 55, and 27.5 nm pixel size PSF for 110 nm pixel size single-molecule data. The simulated data consist of 51 axial positions evenly distributed along $-600$ to 600 nm with 1000 single molecules at each axial position. Each molecule contributes a total of 2000 photons, with a background level of $2\times {10}^{-3}\text{photons}/{\mathrm{nm}}^2$.

Figure 7.CRLB and RMSE improved by using upsampled PSF model and novel spline PSF fitter. (a) Comparison of RMSE for 330 nm pixel data, with and without pixel integration [convolution in Fig. 6(a)]. RMSE_con shows results after convolution, while RMSE_raw shows results without convolution. Only with convolution (RMSE_con) does the localization achieve the theoretical CRLB limit. (b) Comparison of RMSE for localizing 330 nm data with 110 nm upsampled PSF and 330 nm PSF using our novel fitting algorithm. The results show that only the upsampled PSF achieves the theoretical CRLB limit. Simulation parameters match those in Fig. 6.

Figure 8.Comparison of localization accuracy at different lateral positions. Localization accuracy of 110 nm pixel size PSFs on 330 nm single-molecule data at different lateral positions, (0 nm, 0 nm), (0 nm, 115 nm), (115 nm, 0 nm), and (115 nm, 115 nm) for (a), (b), (c), and (d), respectively. Here, the center of the central pixel is set as coordinate (0 nm, 0 nm). We used a vector PSF model (Appendix D) to generate 1000 simulated single molecules at different $x$ and $y$ positions at each axial position. Twenty-five axial positions uniformly distributed along $-600$ to 600 nm were evaluated. Each molecule contributes a total of 2000 photons, with a background level of $2\times {10}^{-3}\text{ photons}/{\mathrm{nm}}^2$. Localization accuracy at each $z$ position is computed from the RMSE of the localization results from 1000 single molecules.

Figure 9.Validation of the upsampled PSF modeling using simulated data from tetrapod PSF. (a) The 330 nm pixel size PSF (top) and 110 nm pixel size PSF (bottom) estimated from 300 simulated tetrapod PSF data sets with a pixel size of 330 nm. (b)–(d) Comparison of the impact of the estimated 330 and 110 nm pixel size tetrapod PSF on the localization accuracy of the 330 nm pixel size simulated data along the $x$, $y$, and $z$ directions, respectively. Each molecule contributes a total of 10,000 photons, with a background level of $2\times {10}^{-3}\text{ photons}/{\mathrm{nm}}^2$. The $x/y/z$ bias as a function of z is defined as in Appendix F, Eqs. (F1) and (F2). Solid lines and shaded areas in localization plots indicate mean and standard deviation of bias, respectively.

Figure 10.Validation of the upsampled PSF modeling using experimental data from astigmatic PSF. (a) 321 nm pixel size experimental bead data (top) and the estimated 321 and 107 nm pixel size astigmatic PSF (bottom). (b)–(d) Comparison of the impact of the estimated 321 and 107 nm pixel size astigmatic PSF on the localization accuracy of the 321 nm pixel size experimental data along the $x$, $y$, and $z$ directions, respectively. The $x/y/z$ bias as a function of $z$ is defined as in Appendix F, Eqs. (F1) and (F2). Solid lines and shaded areas in localization plots indicate mean and standard deviation of bias over $50$ beads, respectively.

Figure 11.Detailed layout of the optical setup. M, mirror; DM, dichroic mirror; L, lens; TS, translation stage; FC, fiber coupler; fiber, single-mode fiber; BFP, back focal plane; FW, filter wheel; TBL, tube lens; AP, aperture; QPD, quadrant photodiode. The excitation lasers are first reflected by dichroic mirror DM1 and coupled into a single-mode fiber through the fiber coupler FC. Before being reflected by the main dichroic mirror DM2 to enter the objective for sample illumination, the beam is collimated and reshaped by a pair of lenses (L1 and L2) and a slit at AP1. In the imaging path, the fluorescence collected by the objective passes through the dichroic mirror DM2 and is filtered by the filter wheel FW. It is then focused using a tube lens TBL. Subsequently, the fluorescence passes through a 4f system composed of lenses L3 ($\text{focal}\text{length}=150\mathrm{nm}$) and L4 ($\text{focal}\text{length}=30\mathrm{nm}$), before being detected by the camera. Additionally, a beam excited by a 785 nm laser, reflected off the coverslip, is detected by the quadrant photodiode (QPD), providing feedback control to the Z-stage for focus locking.

Figure 12.Comparison of the impact of PSF Gaussian blurring on the theoretical CRLB. The $\sigma$ corresponds to the ${\sigma }_x$ and ${\sigma }_y$ terms in Appendix D, Eq. (D8). When $\sigma =0$, Gaussian blurring is not applied. The pixel size of the data used is 330 nm. The simulation parameters are the same as those in Fig. 5.