Structured illumination microscopy with unknown patterns and a statistical prior

Structured illumination microscopy (SIM) improves resolution by down-modulating high-frequency information of an object to fit within the passband of the optical system. Generally, the reconstruction process requires prior knowledge of the illumination patterns, which implies a well-calibrated and aberration-free system. Here, we propose a new \textit{algorithmic self-calibration} strategy for SIM that does not need to know the exact patterns {\it a priori}, but only their covariance. The algorithm, termed PE-SIMS, includes a Pattern-Estimation (PE) step requiring the uniformity of the sum of the illumination patterns and a SIM reconstruction procedure using a Statistical prior (SIMS). Additionally, we perform a pixel reassignment process (SIMS-PR) to enhance the reconstruction quality. We achieve 2$\times$ better resolution than a conventional widefield microscope, while remaining insensitive to aberration-induced pattern distortion and robust against parameter tuning.


Introduction
The Abbe diffraction limit was considered to be the fundamental limit for spatial resolution of an optical microscope for more than a hundred years. In the last decade, novel techniques have circumvented this limit in order to achieve super-resolution [1][2][3][4][5][6][7]. Structured illumination microscopy (SIM) [1][2][3][4], for example, uses illumination by multiple structured patterns to down-modulate high spatial frequency information of the object into the low-frequency region, which can then pass through the bandwidth of the microscope's optical transfer function (OTF) and be captured by the sensor. The reconstruction algorithm for SIM combines demodulation process which brings the high spatial frequency information back to its original position and synthetic aperture that extends the support of the effective OTF. Various structured patterns have been used to realize SIM: periodic gratings [1][2][3][4], a single focal spot (confocal microscope) [8,9], multifocal spots [10][11][12][13] and random speckles [13][14][15][16][17][18][19][20][21][22]. When the illumination patterns themselves are diffraction-limited, linear SIM is restricted to 2× the bandwidth of a widefield microscope [4], allowing up to ∼ 2.4× resolution enhancement (metrics explained in Sec. 3).
In practice, structured illumination systems are sensitive to aberrations and experimental errors. To avoid reconstruction artifacts that degrade resolution, the patterns that are projected onto the sample must be known accurately. Periodic grating patterns can be parameterized by their contrast, period and phase angle, which may be estimated in the post-processing [23][24][25][26]. For multifocal patterns, the location of each focal spot is required [10]. For random speckle patterns, the relative shifts of the patterns are needed [18,19]. Even with careful calibration and high-quality optics, distortions caused by the sample may degrade the result.
To alleviate some of the experimental challenges, blind SIM was proposed, enabling SIM reconstruction without many priors [16,17,21,22,27,28]. The only assumption is that the sum of all illumination patterns is uniform. Optimizationbased algorithms have been adopted, including iterative least squares with positivity and equality constraints [16,21,27], joint support recovery [17] and 1 sparsity constraints [22]. However, these algorithms are sensitive to parameter tuning and may show low contrast in reconstructing high spatial frequencies [16]. Another algorithm, speckle super-resolution optical fluctuation imaging (S-SOFI) realizes SOFI [29] by first projecting random speckle patterns onto the object, and then using the statistical properties of the speckle patterns as a prior to reconstruct a high-resolution image [20]. S-SOFI is experimentally simple and robust; however it only achieves a 1.6× resolution enhancement instead of 2.4× for conventional SIM techniques (as compared to a widefield microscope).
In this paper, we propose a new reconstruction algorithm for SIM that is applicable to any illumination patterns. Our method, termed pattern estimation structured illumination microscopy with a statistical prior (PE-SIMS), is as robust and insensi-tive to parameter tuning as S-SOFI, and achieves better resolution enhancement (up to 2×). Like blind SIM, the patterns need not be known (except for a requirement on the covariance of the patterns). We demonstrate our method using simulated and experimental results with both speckle and multifocal patterns. We discuss pattern design strategies to reduce the amount of data required and demonstrate an extension that uses pixel reassignment [30][31][32][33][34] to improve the reconstruction quality.

Theory and Method
Our algorithm takes in a SIM dataset consisting of multiple images captured under different structured illumination patterns (e.g. random speckles, multifocal spots). We reconstruct the super-resolved image in two parts. The first part is an iterative optimization procedure for estimating each illumination pattern based on an approximated object. The second part reconstructs the high-resolution image using the estimated patterns and the measured images, along with a statistical prior. Before introducing these two parts, we start by defining the SIM forward model. Figure 1: Example experimental setup for structured illumination microscopy (SIM) using a deformable mirror device (DMD) to capture low-resolution images of the object modulated by different illumination patterns. Our IPE-SIMS algorithm reconstructs both the super-resoloved image and the unknown arbitrary illumination patterns.

Forward model of structured illumination microscopy
A representative experimental setup is shown in Fig. 1. A DMD spatial light modulator (SLM) is used to project patterns onto the object through an objective lens.
The measured intensity for the -th captured image is the product of the object's fluorescence distribution o(r) with the illumination pattern p (r), where r = (x, y) denotes the lateral position coordinates. This product is then convolved with the system's incoherent detection-side point spread function (PSF), h det (r):

Part 1: Pattern estimation
The first part of our inverse algorithm is to estimate the illumination patterns. To do so, we start with an low-resolution approximation of the object. Then, we use this object and our measured images to iteratively estimate the patterns (see Fig. 2). Figure 2: The first part of our algorithm, Pattern Estimation (PE), iteratively estimates the illumination patterns from an approximated object given by the deconvolved widefield image.

Part 1a: Approximate widefield image
If we already knew the object o(r), it would be straightforward to estimate the pattern for each measured image by dividing out the object from each of the measurements. However, the object o(r) is unknown. Hence, we start by making a rough estimate of the object. We first take the mean of all the measured images: where · is the mean operation with respect to , and p 0 = p (r) is approximately a constant over the entire field of view. The resulting image will be equivalent to the low-resolution widefield image if the sum of all illumination patterns is approximately uniform.

Part 1b: Deconvolve widefield image
Since the widefield image represents the convolution of the object with its PSF, we can perform a deconvolution operation to estimate the low-resolution object: where F and F −1 denote the Fourier transform and its inverse, respectively,· denotes the Fourier transform of a certain function, u = (u x , u y ) are the lateral spatial frequency coordinates and β is a small Tikhonov regularization constant. Note that this object estimate has diffraction-limited resolution and will be used only for estimating the illumination patterns.

Part 1c: Pattern estimation
We then use the low-resolution object estimate o est (r) to recover each of the illumination patterns. Since each image is simply the product of the illumination and object, we could divide each image by the estimated object to get the pattern. However, we instead solve the problem as an optimization procedure in order to impose the correct Fourier support constraint and avoid reconstruction artifacts. The -th pattern estimate is the solution to the following problem where λ illu is the wavelength of the excitation light. The first term of the cost function, f diff (p ), in Eq. (4) is the least square error (residual) between the measured intensity and the predicted intensity based on our current estimate. The second term enforces a frequency support constraint for the illumination pattern via an indicator function I C . This is important to reduce artifacts in the pattern estimation because a normal division between the measured image and estimated object will create errors outside of this frequency support. In our epi-illumination geometry, the constraint is that the frequency content of each illumination pattern be confined within the OTF defined by the objective's NA. We implement a proximal gradient descent algorithm [35], summarized in Subroutine 1. Proximal gradient descent is designed to solve convex optimization problems like ours that have two cost function terms: one being a differentiable cost function term (e.g. the residual) and the other being a constraint or regularization term (usually nondifferentiable). When the constraint is defined by an indicator function, as in Eq. (4), the method is also known as a projected gradient method.
To implement, we first compute the gradient of the differentiable cost function term with respect to p (r) (5) where k denotes evaluation of the gradient using the pattern at the k-th iteration. We define the projection operation Π C to force the information outside of the OTF to be zero at each iteration. To reduce high-frequency artifacts, the following soft-edge filter is used where h illu (r) is the system's illumination-side PSF, and δ determines the amount of high-frequency information that is suppressed in the pattern estimation step. We repeat this process of updates and projections until convergence (typically ∼50 iterations to estimate each pattern). The convergence speed for proximal gradient descent is on the order of O(1/K) [35], indicating that the residual between the current and optimal cost functions is inversely proportional to the number of iterations K. To accelerate convergence, one extra step is conducted in Subroutine 1 to include the information of the previous estimate [36,37]. The convergence rate for this accelerated proximal gradient method, O(1/K 2 ) [37], is significantly faster than the normal proximal gradient method.

Part 2: SIM with a statistical prior
Once we have recovered the illumination patterns, the second part of the algorithm is to reconstruct a high-resolution image from the measured dataset I (r) and the estimated patterns p (r). We call this part of the algorithm Structured Illumination Microscopy with a Statistical prior (SIMS), summarized in Fig. 3. There are four steps, which are explained below. We will also describe how the statistical prior is used and why this procedure gives better resolution. Figure 3: The second part of our algorithm, termed structured illumination microscopy with a statistical prior (SIMS), estimates the high-resolution object from the measured images and the estimated illumination patterns obtained in Part 1.

Part 2a: Calculate the pattern-intensity covariance
Consider the case where the pattern p(r) is a random variable at position r and the measured intensity I(r) is also a random variable at position r. The -th image is thus the -th sample function for these random variables (one event out of the sample space). Covariance is a measure of how much two random variables change together. Since the intensity I(r) is the blurred version of the product between random patterns p(r) and deterministic object o(r) (Eq. (1)), the covariance between the pattern and the intensity should give high similarity wherever the object o(r) has signal and thus allow us to find the object underneath the random-pattern illumination [14,[38][39][40][41]. We calculate this covariance image I cov (r) as where ∆I (r) = I (r) − I (r) , and ∆p (r) = p (r) − p (r) . Regardless of which illumination pattern is imposed, the covariance image always gives an estimate of the object. However, the resolution of the reconstructed object may be different for different pattern statistics. We can quantify this by taking a closer look at the expression on the right-hand side of Eq. (7). The multiplication of detection PSF and covariance between p(r) and p(r ) acts as the PSF of the covariance image, which thus determines resolution. If the patterns are perfectly spatially correlated, the pattern-pattern covariance is a constant, and the pattern-intensity covariance image is a normal widefield image with PSF of h(r). If the patterns are perfectly spatially uncorrelated, the pattern-pattern covariance is |∆p (r)| 2 δ(r − r ), which, for a constant variance, results in the PSF being a delta function and the object being reconstructed with perfect resolution. In practice, this is not achievable, since the illumination is bandlimited and thus cannot be perfectly uncorrelated. In the general case, to find the resolution (PSF) of the covariance image, we need to calculate the spatial covariance of the patterns, which is the subject of Part 2b, below.

Part 2b: Calculate pattern-pattern covariance
To calculate the spatial covariance of the projected patterns, we first consider the pattern formation model. In our experiments, for example, we use a DMD to create random patterns at the sample plane. Assuming that the projected DMD pattern is sparse enough to avoid interference cross-terms, we can express our pattern under the incoherent model as where t (r) is the -th pattern on the DMD. With this model, the pattern-pattern covariance is where we have used an assumption that the DMD pattern values at position r 1 and r 2 are perfectly uncorrelated: with γ t being a constant that maintains unit consistency. This assumption is valid because the effective DMD pixel size is small compared to the FWHM of the optical system and we can control ∆t (r) to create an uncorrelated pattern. In the experiment, each position of t (r) is an independent and identically distributed random variable. When the number of patterns is large enough, the variance ∆t 2 (r 1 ) approaches the same constant for all the positions. We can then combine γ t and the variance into a single constant α t . Ideally, we can assume h illu (r) ≈ h det (r) when λ illu ≈ λ det , where λ det is the wavelength of the fluorescent emission detection light, and theoretically calculate the pattern-pattern covariance. We can also estimate h illu h illu (r) by numerically evaluating Eq. (9) using our estimated patterns, which accounts for possible aberrations in the illumination optics.

Part 2c: PSF deconvolution of the covariance image
The pattern-pattern covariance derived in Part 2b is related to the PSF of the patternintensity covariance calculated in Part 2a. Hence, we can plug the pattern-pattern covariance into Eq. (7) and write the covariance image as Importantly, the effective PSF for this correlation image is now [(h illu h illu ) · h det ](r), and the corresponding effective OTF is [|h illu | 2 ⊗h det ](u). Since both |h illu | 2 andh det have approximately the same Fourier support as the widefield OTF, the convolution between them covers around 2× the support of the widefield OTF, as in conventional SIM. Given the effective PSF, we implement a standard deconvolution to improve contrast at high spatial frequencies: where H(u) = [|h illu | 2 ⊗h det ](u) and ξ is a small regularization parameter.

Part 2d: Shading correction operation
When the number of images is not large enough to give uniform variance of the patterns at each pixel ( ∆t 2 (r ) from Eq. (9)), low-frequency shading artifacts will occur. Even if we assume the mean of the pattern to be flat in Eq. (2), the variance can still be non-uniform. These can be seen in the deconvolved covariance image in Fig. 3. To resolve this, we can estimate and correct for the variance across the image using our previously estimated projected patterns. Since the projected pattern p (r) is the blurred version of the pattern on the DMD, by ignoring the high-frequency component of the DMD pattern, we can approximate the variance of the DMD pattern by We divide out the spatially-varying variance α t in Eq. (11) from the deconvolved SIMS image, where is a regularizer and I SIMS (r) is the output from our SIMS reconstruction (Part 2c). This result of this step is our final reconstruction of the high-resolution object function.

Parameter Tuning and Algorithm Runtime
Our SIMS algorithm involves 4 regularizers: β, δ, ξ, and , described in Eq. (3), Eq. (6), Eq. (12), and Eq. (14), respectively. Each is decoupled from the others and acts similarly to a typical Tikhonov regularizer, so tuning may be done independently. Generally, we want the regularizers to be as small as possible, while still avoiding noise amplification. The procedure to tune the regularization parameters heuristically is summarized as follows. First, we check if the widefield images are well-deconvolved by finding the smallest β to give the image with best resolution but without obvious noise amplification, then we move on to check the deconvolved covariance image by tuning the SIMS regularizer ξ and the smooth-edge filter regularizer δ using the same principle, and finally we check the final reconstruction by using the smallest shading correction regularizer with enough shading correction but without evident noise amplification. Additionally, the negative values in all of the deconvolved images are set to zero since the fluorescent density is always positive.
The algorithm is implemented in MATLAB and run on an Intel i7 2.8 GHz CPU computer with 16 G DDR3 RAM under OS X operating system. To reconstruct an image with size of 200 × 200 pixels and 400 measurements, this computer takes about 200 seconds. The bottleneck of the algorithm is on the pattern estimation step. The estimation of each pattern takes around 0.5 second.

Definition of resolution
Before introducing and comparing any SIM algorithms, we want to first define the resolution criterion considered in this paper. Resolution of a microscopic image is usually defined by measuring the minimal resolvable distance between two points. Consider a widefield image with detected wavelength λ and numerical aperture N A; the Abbe resolution criterion is then 0.5λ/N A, the full width at half maximum (FWHM) of the widefield PSF. As two points get closer to each other, the contrast between them decreases. Under the separation set by Abbe's limit, two infinitely small points observed under widefield microscope will give an overlapped two-point image with a dip at the center with the contrast equal to 0.01. Hence, the Abbe resolution criterion can be thought of as setting the minimum acceptable contrast between two points at 0.01. We can therefore define the resolution of a microscope or a reconstruction algorithm by measuring the smallest resolvable fine features that have contrast between them of at least 0.01.

Comparison of algorithms
Given this definition of resolution, we quantify the resolution for various algorithms in Fig. 4. The Siemens star test target (o(r, θ) = 1 + cos 40θ in polar coordinates) has varying spatial frequencies along the radius. The resolution of different imaging methods is quantified by reading the minimal resolved period when the contrast reaches 0.01. The effective modulation transfer function (MTF) of each method is shown in Fig. 4b, measured as the contrast of the reconstructed Siemens star image at different radii. Our simulations use a SIM dataset with random patterns, so that we may compare against the previously proposed reconstruction algorithms of blind SIM [16] and S-SOFI [20]. We create N img = 400 speckle-illuminated images from shifted random patterns on the DMD, with shifts of 0.6 FWHM of the PSF across 20 × 20 steps in the x and y directions, respectively. In each pattern, only 10% of the DMD pixels are turned on. This noise-free situation allows us to compare the ideal achieved resolution for the different algorithms. Figure 4a shows the widefield, deconvolved widefield, confocal, and deconvolved confocal images of the Siemens star, as compared to blind SIM [16], S-SOFI [20] and our algorithm. At the bottom, we show the measured effective MTF for each algorithm. In terms of visual effect, S-SOFI [20] gives the least artifacts. To compare resolution, we use our definition of the minimal resolved separation when the contrast drops to 0.01 and summarize the results in Table 1. The enhancement metric gives the ratio resolution improvement over widefield imaging. S-SOFI resolves features down to 1.67 × smaller than the widefield microscope, which is close to the claimed 1.6× in [20], and Blind SIM achieves 1.84× improvement but lower contrast for high-frequencies, which is consistent with [16]. Our PE-SIMS and PE-SIMS-PR (PE-SIMS with pixel reassignment algorithm [30][31][32][33][34] described in Appendix B) algorithms give better resolution compared to other methods. We resolve features down to 1.84× and 2×, respectively, close to the limit set by the deconvolved confocal image. Hence, our method performs the best of the blind algorithms.
Ideally, if we know all the patterns and our spatial modulation covers the full Fourier bandwidth of the objective, we could reconstruct out to 4N A/λ in Fourier space, achieving enhancement of 2.42×, as in the case of deconvovled confocal image or periodic SIM with known patterns. The blind algorithms, however, deal with an ill-posed problem (measure N img images and solve N img + 1 images) that can only become well-posed through appropriate constraints. If the prior for these algorithms are not accurate enough, they may solve a different problem even if the problem becomes well-posed. This is why algorithms with different prior assumptions give different resolution performance for the same dataset, as we saw in Table 1.

Experimental Results
Our experimental setup is shown in Fig. 1. A laser beam (Thorlabs, CPS532, 4.5 mW) is expanded to impinge onto a reflective DMD spatial light modulator (DLP R Discovery 4100, .7" XGA, 1024×768 pixels, pixel size 13.6 µm). The DMD generates a total of N img random patterns (30% of DMD pixels turned on). These random illumination patterns are projected onto the object (with demagnification of 60×) through a 4f system composed of a 200 mm convex lens and a 60× objective lens with NA= 0.8 (Nikon CFI). The resulting fluorescent light is then collected with another 4f system formed by the same 60× objective and a 400 mm convex lens (magnification 120×). A dichroic mirror blocks the reflected illumination light (as in a typical epi-illumination setup). The images are taken with an sCMOS camera (PCO.edge 5.5, 2560×2160 pixels, pixel size 6.5 µm). Patterns are shifted on a 20×20 grid in the x and y directions with a step size of 0.6 FWHM of the PSF, while collecting images at each step. Our test object is carboxylate-modified red fluorescent beads (Excitation wavelength: 580 nm/Emission wavelength: 605 nm) having mean diameter of 210 nm (F8810, Life Technologies).
Reconstruction results are shown in Fig. 5, demonstrating improved resolution using our PE-SIMS algorithm, as compared to standard widefield or deconvolved widefield images.
To quantitatively analyze the experimental results, we measure the resolved feature size of the reconstructed image and compare it to our theory. As shown in the cutline in Fig. 5, two fluorescent beads separated by 328 nm can clearly be resolved using our method, which are otherwise unresolvable in either widefield or deconvolved widefield images. The contrast of this two-Gaussian shape shows these two Gaussian are separated by 1.16× FWHM, so the FWHM of the reconstructed beads is around 283 nm. Assuming the bead can be modeled as a Gaussian function with FWHM of 140 nm (210 nm in diameter for the beads), we can then deconvolve the bead shape out of the reconstruction and get the FWHM of the PSF for this case equal to 240 nm, which is below the diffraction limit λ/2N A = 371 nm.
Our algorithm can be used on other types of SIM datasets, as long as the patternpattern covariance gives a point-like function at the center. As an example, we tested our algorithm on a dataset from a previous method, Multispot SIM (MSIM) [10]. In MSIM, the patterns are a shifting grid of diffraction-limited spots. Since the previous MSIM implementation assumes known patterns, a calibration step captured an extra dataset with a uniform fluorescence sample in order to measure the patterns directly. Our algorithm ignores this calibration data, yet accurately reconstructs both the object and patterns (see Fig. 6). The MSIM result using the calibration data is shown for comparison. The sample is microtubules stained with Alexa Fluor 488 in a fixed cell observed under a TIRF 60× objective with N A = 1.45. Our PE-SIMS-PR reconstruction gives a similar result to the known-pattern MSIM reconstruction.

Conclusion
We have proposed a robust algorithm that can give 2× resolution improvement compared to widefield fluoresence imaging using a SIM dataset without knowing the imposed patterns. Our algorithm first estimates each illumination pattern from a low-resolution approximate object and measured intensities by solving a constrained convex optimization problem. We then synthesize a high-resolution image by calculating the covariance between the estimated patterns and the measured intensity images, followed by a deconvolution and shading correction to get to the final reconstruction. We quantified the limits on resolution of our algorithm by the reconstructed contrast of a simulated Siemens star target. In simulations, we showed that our algorithm gives better resolution compared to previously proposed blind algorithms [16,20]. Experimentally, we demonstrated this improvement experimentally on both random speckle pattern illumination and multi-spot scanned illumination. Figure 6: Comparison of our algorithm on dataset from Multispot SIM (MSIM) which uses withN img = 224 scanned multi-spot patterns from [10]. We show the deconvolved widefield image and the reconstructions using MSIM with known patterns, as well as our blind PE-SIMS algorithm with and without pixel reassignment.
In this paper, we used 400 random speckle illumination patterns to reconstruct the image, far more than the 9-image requirement of conventional SIM [4]. This large number of images was required for high-quality reconstructions because the average and variance of the illumination patterns must be sufficiently flat in order to avoid shading variations. Recall that we want α t (r) ≈ γ t ∆p 2 (r) in Eq. (13) to be close to a constant, which suggests that the variance of the random patterns is constant. When the number of images N img goes down, this statistical assumption is not true any more. We use a shading correction algorithm (Sec. 2.3) to fix this problem by estimating the nonuniform variance, but it is still only an estimate. Hence, when the degree of variance nonuniformity increases (as the number of images decresases), the shading correction algorithm incurs errors. Figure 7 shows simulations demonstrating the effect of reducing the number of images. We use the same random pattern as in Sec. 2.3 and shift by step sizes of 0.6 FWHM of the PSF. As we decrease the number of images from 400 to 36, the reconstruction becomes worse, due to shading errors. The shading map, α t (r)o(r), is shown in the bottom row of Fig. 7. We can see the artifacts happen at the region where the α t (r) is dim and changing. Without knowing the patterns a priori it is not possible to fully correct these shading effects. Figure 7: Results with simulated and experimental (fluorescent beads) datasets comparing random speckle and multi-spot illumination patterns. (middle row) Shading maps overlaid on the object. Decreasing the number of random patterns results in shading artifacts in the reconstruction. The random patterns are scanned in 20 × 20, 10 × 10, and 6 × 6 steps with the same step size of 0.6 FWHM of the PSF, while the multi-spot pattern is scanned with 6 × 6 steps.
Since we know that the artifacts that appear with too few images are due to a non-uniform α t (r), we can attempt to design patterns that will be uniform with a minimal number of images. We would like ∆p 2 (r) to give a uniform map. Consider the contribution from a single pattern; ∆p 2 (r) is similar to the original pattern but with sharper bright spots. The ensemble average over sums up all these bright spots after shifting the pattern around. For a shifted random pattern, we must capture many images in order for the summation of the bright spots to give a uniform map. One efficient way to get a sum of bright spots to become a uniform map is to use a periodic multi-spot pattern (see Fig. 7) [10,11,13]. The period of this multi-spot pattern is designed to be 6 shifting step sizes. Thus, we can use 6 × 6 scanning steps to give a uniform shading map α t (r). The reconstruction is also shown in Fig. 7 to be almost as good as the one illuminated with 400 shifted random patterns.
Experimentally, we see similar trends in image reconstruction quality for different illumination strategies (see the bottom row of Fig. 7). Results from random pattern illumination of fluorescent beads with N img = 400 and multi-spot illumination with N img = 36 give very similar results, and shading artifacts become prominent as the number of patterns is reduced. Note that we use the same algorithm for both the random and multi-spot illuminated datasets because the PSFs of the pattern-intensity covariance images I cov (r) for both cases are the same.
To show that the PSF for the pattern-intensity covariance image with random and multi-spot illumination are the same, we must derive the pattern-pattern covariance ∆p (r)∆p (r ) as we did in Part 2b in Sec. 2.3. To calculate the patternpattern covariance, we need to calculate the covariance of the patterns on the DMD ∆t (r)∆t (r ) and plug it into Eq. (9) to get pattern-pattern covariance ∆p (r)∆p (r ) . For the multi-spot case, we can express the pattern on the DMD and its zero-mean pattern as where r mn = (mΛ, nΛ), m and n are integers, and Λ is the period of the pattern. Then, we can calculate the covariance of the pattern on the DMD as where η is a constant that enforces unit consistency. Plugging this into Eq. (9), we can then calculate the pattern-pattern covariance as ∆p (r)∆p (r ) = (h illu h illu )(r − r ) ⊗ Λ 4 η m,n δ(r − r − r mn ).
Although the pattern-pattern covariance is only a replica of the (h illu h illu )(r), the PSF of the covariance image, I cov (r), only depends on the multiplication of h det (r) and (h illu h illu )(r) ⊗ Λ 4 η m,n δ(r − r mn ) as Eq. (7) derived. If the period of the multi-spot pattern is large compared to (h illu h illu )(r), we can still have our PSF as [(h illu h illu ) · h det ](r), which is the same as the case of random pattern illumination.  In the real space, the PSF after doing pixel reassignment looks fatter than the one without pixel reassignment. However, the OTF of the one with pixel reassignment has larger value in the high-frequency region, where the noise severely degrade the image resolution. Thus, we get better SNR by summing up all the information we have and have a OTF that better deals with noise at high-frequency region. Since we know the PSF, [(h illu h illu ) ⊗ h det ](2r), and the shading map, α t (r), of this pixel-reassigned image I PR (r), we can again apply the deconvolution and the shading correction operation described in Sec. 2.3 to get a PE-SIMS-PR reconstruction. Figure 8(b) compares the reconstruction result of fluorescent beads using 6 × 6 multi-spot illumination with and without applying pixel reassignment algorithm. Pixel reassignment results in sharper contrast when two beads are close to each other and helps clean up some background deconvolution errors. A cut-line plot of the fluorescent beads in Fig. 8(b) shows that the FWHM of the reconstructed bead from SIMS (300.3 nm) is larger than for SIMS-PR with pixel reassignment (254 nm), giving better resolution.