## Abstract

The performance of structured illumination microscopy (SIM) systems
depends on the computational method used to process the raw data. In
this paper, we present a regularized three-dimensional (3D)
model-based (MB) restoration method with positivity constraint (PC)
for 3D processing of data from 3D-SIM (or 3-beam interference SIM), in
which the structured illumination pattern varies laterally and
axially. The proposed 3D-MBPC method introduces positivity in the
solution through the reconstruction of an auxiliary function using a
conjugate-gradient method that minimizes the mean squared error
between the data and the 3D imaging model. The 3D-MBPC method provides
*axial super resolution*, which is not the same as
improved optical sectioning demonstrated with model-based approaches
based on the 2D-SIM (or 2-beam interference SIM) imaging model, for
either 2D or 3D processing of a single plane from a 3D-SIM dataset.
Results obtained with our 3D-MBPC method show improved 3D resolution
over what is achieved by the standard generalized Wiener filter
method, the first known method that performs 3D processing of 3D-SIM
data. Noisy simulation results quantify the achieved 3D resolution,
which is shown to match theoretical predictions. Experimental
verification of the 3D-MBPC method with biological data demonstrates
successful application to data volumes of different sizes.

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

## 1. Introduction

Three-dimensional (3D) structured illumination microscopy (3D-SIM) [1,2], in which the structured illumination (SI) pattern varies
laterally and axially, has become one of the most effective optical
imaging modalities used in biological investigations because of its
optical sectioning and super-resolution capabilities (see review [3] and references therein). 3D-SIM has
been shown to achieve both lateral and axial super resolution by
*3D processing* of the 3D-SIM data set (which consists of a
large number of images acquired by axially-scanning the 3D sample while
changing the phase and orientation of the SI pattern) using a multistep
approach proposed by Gustafsson et al. [1], that employs a non-iterative 3D generalized Wiener filter
(3D-GWF) for the 3D deconvolution step. Two other methods known to provide
3D processing of 3D-SIM data follow a similar multistep approach as in
[1], but instead of the 3D-GWF, they
use an iterative approach for the 3D deconvolution step [4,5]. For 3D-SIM based on three-wave interference [1], fifteen raw images per focal plane are
acquired using phase-diversity (5 different phases) in three orientations
of the SI pattern.

Because SIM is based on a computational sensing paradigm, computational
methods are an integral part of the imaging system modification and have a
direct impact on system performance. Through the inclusion of constraints
in the optimization framework, model-based iterative approaches can be
developed to provide regularized solutions to the 3D-SIM inverse problem
while avoiding unrealistic negative values that are not consistent with
the desired fluorescence intensity of the underlying sample (a problem of
the results obtained with the 3D-GWF method [6]). In our prior work presented in conference
publications [7,8], we developed the first 3D model-based
(3D-MB) iterative approach to solve the 3D-SIM inverse imaging problem.
Our 3D-MB approach is based on a forward imaging model for 3D-SIM that
takes into account axial scanning of the sample, consistent with 3D
data-acquisition in commercial 3D-SIM systems. This method was shown to
provide more accurate 3D restorations than the multistep-based approach
that uses the 3D-GWF (which we will refer to in what follows as the 3D-GWF
method) in both noisy simulation [7]
and with experimental data [8], at
the expense of longer computation time. In this contribution we extend our
prior work by developing a *regularized* 3D-MB method that
*includes a positivity constraint*.

Development of model-based iterative restoration methods has been restricted to two-dimensional (2D) structured illumination microscopy (2D-SIM) [9] instead, in which the SI pattern does not vary axially, and the inverse problem can be solved easily by 2D processing of the data [3], which is faster and has modest computer memory requirements compared to 3D processing. Although in 2D-SIM axial scanning of the sample is still performed, the acquired data set (9 images per focal plane of the sample due to phase diversity introduced in the SI pattern at 3 different pattern orientations) is processed with 2D restoration methods to provide the best solution consistent with a 2D forward imaging model. Several model-based restoration methods for 2D-SIM have been developed based on different optimization approaches including constraints for regularization and positivity (see [10–17] and references therein), providing in some cases joint estimation of the SI patterns and corresponding parameters (such as, rotation angle and phase) [11–14,18] and data-acquisition reduction by taking advantage of redundancies in the forward images [10,13,19–21]. Some model-based optimization algorithms used in 2D-SIM achieved positivity of the solution by either including directly a positivity constraint [15] or through the reconstruction of an auxiliary function [18], as we do here.

Direct extension of the 2D imaging model to 3D has enabled 3D processing of 2D-SIM data [11,15]. However, for 3D processing of 3D-SIM data, the correct 3D imaging model must account that the axial illumination term modules the system’s point-spread function as was shown in [1], in order to avoid 3D restoration artifacts [7,8]. This model difference has prevented an easy extension of model-based restoration methods developed for 2D-SIM to methods that are suitable for 3D restoration from 3D-SIM data. A 2D-processing method based on the 2D-SIM imaging model and data preprocessing to remove out-of-focus light has been proposed to restore single slices obtained with a modified 3D-SIM system from a thin sample [22]. 3D processing of a single plane from a 3D-SIM dataset, with model-based approaches based on the 2D-SIM imaging model was shown to reduce out-of-focus light in the 2D image [11,15]. We have shown that 3D processing of a 3D experimental volume (obtained from a commercial SIM microscope) with the wrong 3D model (i.e., direct extension of the 2D SIM model to 3D) results in restoration artifacts, that are more evident axially [8]. The use of the correct 3D imaging model in solving the 3D inverse problem for 3D-SIM is crucial.

As mentioned above, all previously developed model-based algorithms applied to 3D-SIM data rely on either 2D or 3D processing of a single plane from the dataset, but they do not use the correct 3D imaging model as in our method [7]. Recently, a new iterative 3D restoration method for 3D-SIM [23], developed based on the same 3D-MB framework as the one we presented in prior work [7,8], proposed a different optimization method (a non-linear gradient decent) to minimize the cost function without a regularization constraint. Positivity of the solution was achieved by setting negative values to zero at the end of each iteration, while the performance of the method was evaluated only in simulation.

In this paper, we develop based on the correct 3D imaging model for 3D-SIM,
a regularized 3D restoration method, which enforces positivity of the 3D
solution through the reconstruction of an auxiliary function using a
conjugate-gradient method that minimizes the accumulative mean squared
error between the 3D-SIM data volume and the 3D imaging model. Preliminary
results from this work presented at a conference [24], demonstrated in simulation that the inclusion of the
positivity *a priori* information improves the accuracy in
the solution of the 3D-SIM inverse imaging problem. Here results from both
simulated and experimental 3D-SIM data verify the ability of our 3D-model
based algorithm and quantify the achieved 3D resolution, which is shown to
match theoretical predictions computed from the actual imaging conditions
used. To our knowledge, the method presented in this contribution is the
first regularized 3D model-based algorithm to directly solve the 3D-SIM
inverse problem using positivity *a priori* information and
experimental data from a commercial microscope.

## 2. Model and method

#### 2.1 Forward imaging model of 3D-SIM

In a 3D-SIM system, as previously described in [1,7], the intensity in the 3D image recorded using axial scanning of a 3D sample labeled with fluorescence, can be modeled as:

*z*are the transverse and axial coordinates respectively; $f({\boldsymbol{x},z} )$ is the density distribution of fluorophores within the sample; ${i_k}(z )$ and ${j_k}(\boldsymbol{x} )$ are the separated axial and lateral functions of the SI pattern, respectively; $h({\boldsymbol{x},z} )$ is the point spread function (PSF) of the widefield microscope and

*A*is the 3D convolution operator defined in Eq. (1). For the three-wave (3W) interference 3D-SIM system (3W-SIM) [1], the SI pattern is given by

The SI pattern is a separable function and it can be rewritten in the following form [1]:

#### 2.2 Restoration method

A 3D-MB approach can be used to restore the 3D image of the biological sample, $f({\boldsymbol{x},z} )$, from the acquired 3D data. Since the desired fluorescence distribution of the underlying sample is non-negative, we propose in this paper, a 3D model-based with positivity constraint (3D-MBPC) approach using the conjugate-gradient descent algorithm to reconstruct an estimate of $f({\boldsymbol{x},z} )$ through the reconstruction of an auxiliary function $\zeta $ such that:

Solving the inverse imaging problem in terms of the auxiliary function first, guarantees that the desired final solution computed using Eq. (2) is always non-negative. Positivity of the solution can be achieved in this manner without the need of applying it in the form of a penalty term as in [25]. Because the inverse problem is ill-posed, we introduce a regularization term to the cost function using a weight, as it is commonly done, that allows us to control the tradeoff between resolution and smoothness of the result.There are numerous solvers for non-negative optimization of a convex
cost function. In this paper, we choose the standard method of
optimizing the L^{2}-norm cost function with Tikhonov
regularization, because the convergence of this method has been
studied extensively [26] in
solving the Fredholm equation of the first kind (Eq. (1)). The reconstruction is
obtained by minimizing the functional:

^{2}-norm; $ g_l^{mes}\; $ and ${g_l}\; \; $ are, respectively, the ${l^{th}}$ 3D recorded 3D-SIM image and model prediction from Eq. (1); while $\alpha $ and ${{\cal R}}(f )$ are the regularization parameter and regularization function, respectively.

*L*is the total number of raw 3D-SIM images acquired for different pattern orientations and phases (in general

*L =*15 for three-wave interference 3D-SIM [1]). In this work, the regularization function is taken to be the L

^{2}-norm, i.e., ${{\cal R}}(f )\; = \; |{|f |} |$

^{2}.

The gradient of the cost function (Eq. (3)) with respect to $\zeta ({x,z} )$ is computed as:

*A*, defined by

*n*= 1, …,

*N*: with the step size ${\beta _n}$ determined at each iteration

*n*by minimizing the cost function with respect to ${\beta _n}$, i.e., $F({{{\hat{\zeta }}_n} + {\beta_n}{d_n}} ),$ and the updating direction, ${d_n}$, computed by the well-known Polak-Ribière conjugate gradient method [27]:

*n*

^{th}iteration, and $\langle{\cdot} | \cdot \rangle$ is the inner product. For the cost function in Eq. (3), minimizing $F({{{\hat{\zeta }}_n} + {\beta_n}{d_n}} )$ gives a third-order polynomial with respect to ${\beta _n}$, which can be solved analytically for ${\beta _n}.$ The initial guess, ${\hat{\zeta }_0}$, is set equal to the square root of the widefield (WF) image, which is obtained from the sum of the 15 raw 3D-SIM images. We choose the WF image as the initial guess as it can be easily computed, and also because SIM is expected to improve the resolution of the conventional microscope, which produces the WF image. The complete algorithm for the 3D-MBPC method and information about the code availability are provided in the Supplement 1 (Algorithm S1).

## 3. Simulated results

#### 3.1 Simulated data

The simulation is for the 3W-SIM [1] using a 63x/1.4 NA oil-immersion (refractive index
n = 1.515) lens with an emission wavelength $\lambda $ = 515 nm. For
these parameters, the lateral limit (based on the Rayleigh resolution
limit) and axial resolution limit for the widefield (WF) microscope
are:
*dx *= 0.61*λ*/NA = 224
nm and *dz *= $2\lambda
/\textrm{N}{\textrm{A}^2}$ = 525 nm,
respectively. The predicted SIM lateral and axial resolution limits
are
*d*_{x-SIM }= *dx*/(1+(*u*_{m}/*u*_{c})) = 124
nm (with ${u_m}$ = 0.8${u_c}$) and
*d*_{z-SIM }= *dz*/(1+(*w*_{m}/*w*_{c})) = 292
nm, respectively, with $\frac{{{w_m}}}{{{u_m}}}$ = $\frac{{{w_c}}}{{{u_c}}}$ = $\frac{{\textrm{NA}}}{4}
= $ 0.35, where ${u_m}$ and $ {w_m}\; $ are the lateral and axial modulations
frequencies of the 3D structured illumination pattern generated using
3W-SIM, and ${u_c}$ and ${w_c}\; $ are the lateral and axial cut-off
frequencies of the native WF system, respectively. Using these
parameters and a 3D point-spread function computed using the Gibson
and Lanni model [28], simulated
data was computed using Eq. (1) at three orientation angles ($\theta $ = 0°,
60°, 120°) of the SI pattern in order to achieve
isotropic resolution, with 5 SI pattern phases $\varphi $, shifted by a $({2\pi /5} )$-rads step starting with $\varphi = 0$ rads, at each orientation.

To investigate the performance of the 3D-MBPC approach and compare it to the performance of other approaches, we use a synthetic object as in Fig. 1. The numerical object is simulated on a 512×512×512 grid of a cube with a side equal to 10.22 $\mathrm{\mu }\textrm{m}$. The object is a symmetric 3D star-like object with equally spaced spokes, with the length of each spoke being approximately 3 $\mathrm{\mu }\textrm{m}$. In the XY plane of the object at the middle Z slice, there are 24 spokes. The lateral resolution limit is defined here as the center-to-center distance between two neighboring spokes along a circle centered at the center of the XY-section image. In Fig. 1(b), the red circle is for the lateral resolution limit of the conventional microscope while the blue circle is for the predicted lateral resolution limit of the simulated 3W-SIM. On the other hand, the axial resolution limit is defined as the center-to-center distance between two neighboring spokes along the arc shown in Fig. 1(c); the yellow arc is for widefield microscopy while the green arc is for 3W-SIM. We apply the forward imaging model Eq. (1) to the object to obtain the raw 3W-SIM data (Fig. 1(d) & (e)) and then down-sample it to a 256×256×256 grid to simulate the effect of data recording with a CCD camera (which has a finite pixel and grid size). In addition, Poisson noise [29] is applied to the simulated data at the SNR level of 15 dB (a level found to be reasonable in a previous simulation study [7]).

#### 3.2 3D-MBPC method performance

Using the simulated data from the 3D star-like object described in Section 3.1 we investigate the performance of the 3D-MBPC method as a function of the number iterations, starting with its result at iteration zero, which is equal to the initial guess. As iterations increase to 300, the predicted theoretical SIM resolution limits can be achieved as evident in the images shown in the supplementary Fig. S1. Additionally, the MSE and SSIM values (Table 2), computed between the 3D-MBPC result at selected iterations and the true object, quantify the improvement in the restoration as the number of iterations increases.

#### 3.3 Comparison of different restoration methods

We compare the performance of our 3D-MBPC method to the performance of
the traditional 3D-GWF method, the first method known to perform 3D
processing of 3D-SIM data. In order to evaluate the benefit of 3D
processing over 2D processing, we also process each slice of the 3D
volume using the 2D-GWF implemented in the open source software
fairSIM [30,31] (referred to here as 2D-FairSIM)
and stack the restored 2D images to generate a 3D restoration result.
To compare the results of the different approaches, the restored 3D
images of the star-like object (Fig. 2) are first normalized using the
*l*^{2}-normalization and negative values in
the solution of the 2D-FairSIM and 3D-GWF results are then set to 0.
The restoration parameters used in the fairSIM and 3D-GWF methods, as
well as the regularization parameter and number of iterations used in
the 3D-MBPC method (reported in Fig. 2) have been manually optimized to show the best
results in each case. For example, inspection of results obtained at
different number of iterations of the 3D-MBPC method (supplementary
Fig. S1) enables selection of the number of iterations required for
the final result. The mean square error (MSE) [32] and the structural similarity index measure
(SSIM) [33] computed between
the restoration and the true object are used to compare the
performance of the different computational approaches quantitatively.
The metrics in Table 3
show that the 3D-MBPC result provides the best MSE and SSIM values,
compared to those obtained from the 2D-FairSIM and 3D-GWF results. In
addition to this comparison in the space domain, frequency domain
analysis confirms that the 3D-MBPC results show better frequency
contrast compared to the results from the 2D-FairSIM and 3D-GWF
methods, which are more prone to noise amplification (supplementary
Fig. S4).

The restored results in Fig. 2 verify the theoretical prediction for the lateral and axial resolution limits denoted by the red and yellow arcs, respectively, for widefield microscopy (Fig. 2(b)) and by the blue and green arcs, respectively, for the 3W-SIM limits (Fig. 2(c), (d), & (e)). Qualitative inspection of the XY-section images (Fig. 2, top row) shows that the lateral theoretical resolution limit is achieved in the 2D-FairSIM, 3D-GWF, and 3D-MBPC results. However, intensity profiles taken along the blue arc show that the spokes are not resolved at the theoretical lateral limit in the 2D-FairSIM result (Fig. 2(f)). Furthermore, as expected, only the 3D restoration methods (3D-GWF and 3D-MBPC) can achieve the theoretically predicted axial resolution limit for 3W-SIM (see Section 3.1). In addition, the XZ-section image from the 2D-FairSIM result shows artifacts (black stripes) along the Z direction. Lateral (Fig. 2(f)) and axial (Fig. 2(g)) intensity profiles from the 3D-GWF and 3D-MBPC results show better contrast than those from the 2D-FairSIM result. Finally, the axial intensity profiles (Fig. 2(g)) show that the spokes are better resolved in the 3D-MBPC result (Fig. 2(e)) compared to the 3D-GWF result (Fig. 2(d)).

## 4. Experimental results

#### 4.1 Experimental settings

Two experimental datasets are used in this study, taken from the publicly available fairSIM project data repository [30]. Both datasets were captured using a Delta Vision OMX v4 3D-SIM system (Applied Precision, GE Healthcare) with a 60x/1.42 NA oil-immersion lens. Table 4 summarizes parameters used for reconstruction provided at the fairSIM project data repository. The lateral modulation frequency, ${u_m}$, of the structured illumination pattern is determined from the Fourier transform of the raw data by measuring the distance between the outmost side peak and the central peak. Similarly, the axial modulation frequency, ${w_m}$, is computed from the experimental optical transfer function (OTF) data provided by the fairSIM project developer [34].

The experimental OTF data contains both the conventional OTF, $H({u,v,w} )$, which is the Fourier transform (FT)
of the widefield PSF, $h({\boldsymbol{x},z}
)$, and the modulated OTF, which is the
FT of the axially-modulated PSF, i.e. ${H_m}({u,v,w} )= \;
$ FT $\{ h({\boldsymbol{x},z}
)$ $cos ({2\pi {w_m}z}
)$}, where $({u,v,w} )$ are the 3D spatial frequencies. The
experimental OTF data are given in the frequency domain on a grid size
of 129 × 129 × 65 voxels, with the lateral and axial
pixel sizes of 4.88e^{-2} cycles/$\mu $m and 1.23e^{-1} cycles/$\mu $m, respectively. Since the 3D-SIM
datasets are on a different grid size than the experimental OTF data,
we used spline interpolation to interpolate the experimental OTF data
to the corresponding dataset’s grid size. One can compute the
experimental PSF, $h({\boldsymbol{x},z}
)$, by taking the inverse Fourier
transform of the conventional experimental OTF. From the experimental
OTF data, we computed the ratio $\textrm{IFT}\{
{H_m}({u,v,w} )\} /\textrm{IFT}\{{H({u,v,w} )} \}$ (on a discrete set of non-zero values
of the denominator, where IFT{·} is the inverse
Fourier transform operator, to get the $cos({2\pi {w_m}z}
)$ term, from which we determined the
axial modulation frequency, ${w_m}$, by measuring the distance between
the outmost side peak and the central peak in the FT of the cosine
term. From these imaging parameters, the theoretical lateral and axial
resolutions for widefield microscopy are
*dx *= 0.61λ/NA = 226 nm
and *dz *= $2\lambda
/\textrm{N}{\textrm{A}^2}$ = 521 nm,
respectively. Similarly, the predicted SIM lateral and axial
resolution limits are
*d*_{x-SIM }= *dx*/(1+*u*_{m}/*u*_{c}) = 125 nm
(with ${u_m}$ = 0.8${u_c}$) and
*d*_{z-SIM }= *dz*/(1+*w*_{m}/*w*_{c}) = 330 nm,
respectively, where ${w_m} = \;
0.58$ *w*_{c} and
*w*_{c} = 1/*dz*.

#### 4.2 Results from a small data volume

In this section, we compare restoration results obtained with the 3
different methods described in Section 3.3, from the “OMX LSEC Actin
525 nm” dataset, which has only 7 axial slices
(Table 4). The Wiener
parameter used in the fairSIM and 3D-GWF methods and the
regularization parameter and number of iterations used for the 3D-MBPC
method have been manually optimized to show the best results in each
case. For a quantitative comparison of the results, the restored 3D
images are first normalized using
*l*^{2}-normalization and then negative values
in the 2D-FairSIM and 3D-GWF results are set to 0.

Figure 3 shows that although the results from the 2D processing method using the fairSIM software (Fig. 3(b)) are consistent with the results from the 3D restoration methods (Fig. 3(c)&(d)), the contrast in the latter is better and therefore fine features are better resolved. Furthermore, the results from the 3D restoration methods show better optical sectioning than the results from 2D processing, as evidenced by the features pointed by the white and yellow arrows in Fig. 3(e), (f) & (g), which look different at the three different depth slices in the 3D-GWF and 3D-MBPC results, but they are almost the same in the corresponding 3 depth slices in the 2D-FairSIM result. Additionally, the yellow arrows in the bottom row of Fig. 3(e), (f) & (g) point to fine features that are better resolved in the 3D-MBPC result compared to the results from the other two methods. Frequency domain analysis of these results (supplementary Fig. S6) shows that a better frequency contrast is evident in the Fourier transform of the 3D-MBPC result compared to the Fourier transform of the results obtained from the other methods.

#### 4.3 Results from a large data volume

To assess axial resolution achieved with 3D restoration, and to further demonstrate the advantages of 3D processing over 2D processing of 3D-SIM data, we compare restoration results obtained from the “OMX U2OS Actin 525 nm” dataset (Table 4), which has 53 axial planes, using the 3 methods discussed in the previous sections. The size of each 3D image in the raw data is 40.96 µm × 40.96 µm × 6.625 µm, on a 512×512×53 grid. The entire raw dataset consists of 15 3D images of this size due to the 3 orientations and 5 phases of the SI pattern used to capture the data at each focal plane within the sample. A comparison between the results obtained from this dataset using the 3D-MBPC method and the 3D-MB method [7,8], which does not include the positivity and regularization constraints, shows the benefits of using the two constraints in solving the 3D-SIM inverse imaging problem (supplementary Fig. S3).

Figure 4 compares results
obtained from 2D restoration (Fig. 4(b)) and from the two 3D restoration methods
(Fig. 4(c) & (d).
Restoration parameters (e.g., Wiener parameter used the 2D-FairSIM and
3D-GWF methods or regularization parameter and number of iterations
used with the 3D-MBPC method) were manually optimized to show the best
results obtained with each method. As in the other two studies
presented earlier, the restored 3D images are first normalized using
the *l*^{2}-normalization and then negative
values in the 2D-FairSIM and 3D-GWF results are set to 0. For example,
visual inspection of results obtained at different number of
iterations of the 3D-MBPC method (supplementary Fig. S2) enables
selection of the number of iterations for the final result. In this
comparison study, the difference between the results obtained from the
different methods is more evident, especially in the XZ-section images
(bottom row in Fig. 4(a-d)). The white and yellow arrows in Fig. 4 (e), (f) and (g), which show
XY-section images at different Z depths from the 3D restored volumes,
point to the improved lateral and axial resolution achieved by the 3D
restoration methods over the 2D-FairSIM method, and moreover, by our
3D-MBPC method over the traditional 3D-GWF method. The improved axial
resolution achieved by 3D processing of the data is quantified in
Fig. 5. The axial
intensity profiles in Fig. 5 show that the predicted 3W-SIM axial super-resolution of
330 nm (see Section 4.1) is only achieved by the 3D-MBPC method. These restoration
results demonstrate that our 3D-MBPC approach can provide improved
resolution and optical sectioning over the standard 3D-GWF method.
This result is also confirmed by a frequency domain analysis
(supplementary Fig. S5) in which better frequency contrast is evident
in the 3D-MBPC result compared to the results from the other methods,
which are more prone to noise amplification, as evident by larger
uniform values in high frequencies. Additionally, more high frequency
bands are evident in the Fourier transform of the 3D MBPC result,
where the highest frequency attained is closest to the predicted
frequency limit of 3D-SIM, than in the Fourier transform of the
results obtained with the other two methods.

## 5. Discussion and conclusion

The 3D model-based with positive constraint (3D-MBPC) algorithm demonstrated here provides 3D processing of an entire data volume acquired with 3D-SIM (based on a 3D SI pattern that varies both laterally and axially). To our knowledge this is the first time this 3D inverse problem has been solved using an optimization method with a positivity constraint, which does not require any preprocessing demodulation steps as in the case of the traditional 3D generalized Wiener filter (3D-GWF) method. The most important part of our contribution is that we demonstrate quantitatively the ability to achieve the 3D-SIM potential of axial super resolution as predicted by theory, in both simulation and experiment, using the 3D-MBPC method. Moreover, this work shows that the 3D model-based optimization approach provides more accurate results than the ones obtained with the standard 3D-GWF method, which yields negative values that are not consistent with fluorescence intensities and need to be mitigated.

Our choice of optimization method and constraints is based on approaches we have used before, and we know that they work well in solving ill-posed inverse problems in 3D microscopy. Other optimization methods and constraints can be used to solve the 3D inverse problem using the correct 3D imaging model for the 3D-SIM system (Eq. (1)). The accuracy of the 3D restoration results can be compromised if a 2D restoration method is used to process 3D-SIM data. Here, we have shown this problem by results obtained from 3D-SIM data using the 2D-GWF approach, implemented in fairSIM, and comparing them to the results obtained with the 3D-GWF method.

Our optimization method computes the gradient and step-size analytically, and although we have optimized its implementation and we use multiple high-performance computing (HPC) clusters available to us to process the data, the performance of the 3D-MBPC method currently implemented in MATLAB (MathWorks 2018a, Natick, MA, USA) is slow. For example, it took 12 hours to obtain the results of the small data volume (Section 4.2) and 20 hours for the large data volume (Section 4.3) using 12 Intel Skylake Gold 6148 Processors with 84 GB DDR4 RAM, and EDR InfiniBand in our HPC clusters. There are different ways to improve the processing speed of our 3D-MBPC method, for example, one can estimate the gradient and step-size using the LBFGS optimization method [35].

In this contribution we implemented and applied the 3D-MBPC method to 3D-SIM data from a system in which the illumination is created via 3-beam interference [1]. However, the 3D-MBPC method is suitable for use with any 3D-SIM system in which the 3D SI pattern is separable to a lateral and axial function (Eq. (1)), such as the prism-based 3D-SIM system we have been developing [2], as we have recently shown in simulation [36].

In this work, the 3D-MBPC method uses the correct phases and orientations of the SI pattern. Our work can be extended to estimate simultaneously the object intensity and the illumination parameters, using the alternating optimization approach implemented in solving the 2D-SIM inverse problem [11]. Another possible method is using automatic differentiation [37] to compute the gradient (of both the object intensity and the illumination parameters) and a line search for the step-size computation. Finally, a deep learning approach recently applied to SIM [38] does not require any information of the model parameters but it provides only 2D processing of 3D-SIM data. Extension of this effort to 3D processing, as well as integration of the 3D imaging model with a neural network to develop a physics-guided deep learning method for 3D-SIM is expected to provide more accurate and interpretable results [39]. We will report on this effort in a future publication.

## Funding

Directorate for Biological Sciences (DBI-1353904).

## Acknowledgements

The authors thank Dr. Hasti Shabani (former PhD student at the Computational Imaging Research Laboratory, the University of Memphis) for prior work related to this project and Dr. Anne Sentenac (Institut Fresnel, Marseille, France) for her suggestion to use this restoration framework in 3D-SIM and her collaboration on prior related work. The authors also thank Drs. Michael Unser and Emmanuel Soubies (Biomedical Imaging Group at the Swiss Federal Institute of Technology in Lausanne -EPFL) for the MATLAB code that generates the star-like object. Last but not least, the authors thank Dr. Marcel Müller for the experimental OTF and the necessary parameters used in the restoration of the experimental data, as well as for helpful discussions.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## Data availability

Data underlying the results presented in this paper are available online (see [30]).

## Supplemental document

See Supplement 1 for supporting content.

## References

**1. **M. G. L. Gustafsson, L. Shao, P. M. Carlton, C. J. R. Wang, I. N. Golubovskaya, W. Z. Cande, D. A. Agard, and J. W. Sedat, “Three-dimensional
resolution doubling in wide-field fluorescence microscopy by
structured illumination,” Biophys.
J. **94**(12),
4957–4970
(2008). [CrossRef]

**2. **A. Doblas, H. Shabani, G. Saavedra, and C. Preza, “Tunable-frequency
three-dimensional structured illumination microscopy with reduced
data-acquisition,” Opt.
Express **26**(23),
30476 (2018). [CrossRef]

**3. **R. Heintzmann and T. Huser, “Super-resolution
structured illumination microscopy,”
Chem. Rev. **117**(13),
13890–13908
(2017). [CrossRef]

**4. **J. Atkins, J. Yin, K. Chu, P. Goodwin, P. J. McMillan, S. Wachsmann-Hogiu, S. Lane, and Z. J. Smith, “Image reconstruction
for structured-illumination microscopy with low signal
level,” Opt. Express **22**(7),
8687–8702
(2014). [CrossRef]

**5. **H. Wang, J. Liao, S. Lang, Y. Gong, and Y. Zhang, “Super-resolution
algorithm based on Richardson–Lucy deconvolution for
three-dimensional structured illumination
microscopy,” J. Opt. Soc. Am.
A **36**(2),
173–178
(2019). [CrossRef]

**6. **C. Karras, M. Smedh, R. Förster, H. Deschout, J. Fernandez-Rodriguez, and R. Heintzmann, “Successful optimization
of reconstruction parameters in structured illumination microscopy
– a practical guide,” Opt.
Commun. **436**,
69–75 (2019). [CrossRef]

**7. **H. Shabani, S. Labouesse, A. Sentenac, and C. Preza, “Three-dimensional
deconvolution based on axial-scanning model for structured
illumination microscopy,” in
Proceedings of the International Symposium on Biomedical
Imaging (ISBI2019), pp.
552–555.

**8. **C. T. S. Van, H. Shabani, and C. Preza, “Experimental
verification of 3D model-based restoration for 3D-SIM with data
reduction,” inImaging and Applied Optics Congress, 2020 OSA Technical
Digest Series (Optical Society of
America, 2020), paper
CF4C.3.

**9. **M. G. L. Gustafsson, “Surpassing the lateral
resolution limit by a factor of two using structured illumination
microscopy,” J. Microsc. **198**(2),
82–87 (2000). [CrossRef]

**10. **F. Orieux, E. Sepulveda, V. Loriette, B. Dubertret, and J. C. Olivo-Marin, “Bayesian estimation for
optimized structured illumination microscopy,”
IEEE Trans. Image Process. **21**(2),
601–614
(2012). [CrossRef]

**11. **A. Jost, E. Tolstik, P. Feldmann, K. Wicker, A. Sentenac, and R. Heintzmann, “Optical sectioning and
high resolution in single-slice structured illumination microscopy by
thick slice blind-SIM reconstruction,”
PLoS One **10**(7),
e0132174 (2015). [CrossRef]

**12. **S. Labouesse, A. Negash, J. Idier, S. Bourguignon, T. Mangeat, P. Liu, A. Sentenac, and M. Allain, “Joint reconstruction
strategy for structured illumination microscopy with unknown
illuminations,” IEEE Trans. on Image
Process **26**(5),
2480–2493
(2017). [CrossRef]

**13. **F. Orieux, V. Loriette, J. C. Olivo-Marin, E. Sepulveda, and A. Fragola, “Fast myopic 2D-SIM
super resolution microscopy with joint modulation pattern
estimation,” Inverse Probl. **33**(12), 125005
(2017). [CrossRef]

**14. **L.-H. Yeh, L. Tian, and L. Waller, “Structured illumination
microscopy with unknown patterns and a statistical
prior,” Biomed. Opt. Express **8**(2), 695
(2017). [CrossRef]

**15. **E. Soubies and M. Unser, “Computational
super-sectioning for single-slice structured-illumination
microscopy,” IEEE Trans. Comput.
Imaging **5**(2),
240–250
(2019). [CrossRef]

**16. **J. Luo, C. Li, Q. Liu, J. Wu, H. Li, C. Kuang, X. Hao, and X. Liu, “Super-resolution
structured illumination microscopy reconstruction using a
least-squares solver,” Front.
Phys. **8**, 118
(2020). [CrossRef]

**17. **W. Yu, Y. Li, S. Jooken, O. Deschaume, F. Liu, S. Wang, and C. Bartic, “Second-order optimized
regularized structured illumination microscopy (sorSIM) for
high-quality and rapid super resolution image reconstruction with low
signal level,” Opt. Express **28**(11), 16708
(2020). [CrossRef]

**18. **R. Ayuk, H. Giovannini, A. Jost, E. Mudry, J. Girard, T. Mangeat, N. Sandeau, R. Heintzmann, K. Wicker, K. Belkebir, and A. Sentenac, “Structured illumination
fluorescence microscopy with distorted excitations using a filtered
blind-SIM algorithm,” Opt.
Lett. **38**(22),
4723 (2013). [CrossRef]

**19. **J. Boulanger, N. Pustelnik, L. Condat, L. Sengmanivong, and T. Piolot, “Nonsmooth convex
optimization for structured illumination microscopy image
reconstruction,” Inverse
Problems **34**(9),
095004 (2018). [CrossRef]

**20. **S. Dong, J. Liao, K. Guo, L. Bian, J. Suo, and G. Zheng, “Resolution doubling
with a reduced number of image acquisitions,”
Biomed. Opt. Express **6**(8), 2946
(2015). [CrossRef]

**21. **F. Ströhl and C. F. Kaminski, “Speed limits of
structured illumination microscopy,”
Opt. Lett. **42**(13),
2511 (2017). [CrossRef]

**22. **P. Vermeulen, H. Zhan, F. Orieux, J. -C. Olivo-Marin, Z. Lenkei, V. Loriette, and A. Fragola, “Out-of-focus background
subtraction for fast structured illumination super-resolution
microscopy of optically thick samples,”
J. Microsc. **259**(3),
257–268
(2015). [CrossRef]

**23. **C. Kuang, H. Zhu, J. Han, L. Yin, M. Cai, Q. Liu, X. Hao, X. Liu, Y. Sun, C. Kuang, C. Kuang, C. Kuang, X. Hao, X. Liu, and X. Liu, “3D super-resolution
microscopy based on nonlinear gradient descent structured
illumination,” Opt. Express **29**(14),
21428–21443
(2021). [CrossRef]

**24. **C. T. S. Van and C. Preza, “Three-dimensional (3D)
model-based restoration for structured illumination microscopy based
on a 3D illumination pattern,” in
Imaging and Applied Optics Congress, 2021 OSA Technical
Digest Series (Optical Society of
America, 2021), paper 3Th4D.
1.

**25. **W. A. Carrington, R. M. Lynch, E. D. W. Moore, G. Isenberg, K. E. Fogarty, and F. S. Fay, “Superresolution
three-dimensional images of fluorescence in cells with minimal light
exposure,” Science **268**(5216),
1483–1487
(1995). [CrossRef]

**26. **C. Groetsch, * The Theory of Tikhonov
Regularization for Fredholm Equations of the First
Kind* (Pitman
Publishing, 1984).

**27. **E. Polak and G. Ribiere, “Note sur la convergence
de méthodes de directions
conjuguées,” ESAIM: Math. Model.
Numer. Anal. - Modélisation Mathématique Anal.
Numérique **3**(R1),
35–43
(1969).

**28. **S. F. Gibson and F. Lanni, “Experimental test of an
analytical model of aberration in an oil-immersion objective lens used
in three-dimensional light microscopy,”
J. Opt. Soc. Am. A **8**(10), 1601
(1991). [CrossRef]

**29. **D. L. Snyder, R. L. White, and A. M. Hammoud, “Image recovery from
data acquired with a charge-coupled-device
camera,” J. Opt. Soc. Am. A **10**(5), 1014
(1993). [CrossRef]

**30. **fairSIM,
The fairSIM project, Github,
2021, https://www.fairsim.org/.

**31. **M. Müller, V. Mönkemöller, S. Hennig, W. Hübner, and T. Huser, “Open-source image
reconstruction of super-resolution structured illumination microscopy
data in ImageJ,” Nat. Commun. **7**(1), 10980
(2016). [CrossRef]

**32. **Z. Wang and A. C. Bovik, “Mean squared error:
Love it or leave it? A new look at signal fidelity
measures,” IEEE Signal Process.
Mag. **26**(1),
98–117
(2009). [CrossRef]

**33. **Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality
assessment: from error visibility to structural
similarity,” IEEE Trans. on Image
Process. **13**(4),
600–612
(2004). [CrossRef]

**34. **M. Müller, Biophotonics Group, Bielefeld University:
Bielefeld, Nordrhein-Westfalen, Germany (personal communication,
2020).

**35. **J. Nocedal and S. J. Wright, * Numerical
Optimization*
(Springer,
1999).

**36. **J. Mohammed, C. T. S. Van, and C. Preza, “Evaluation of different
illumination designs for tunable 3D structured illumination microscopy
through model-based restoration,” Proc.
SPIE **11649**, 116490M
(2021).

**37. **L. B. Rall, * Automatic Differentiation:
Techniques and Applications*
(Springer-Verlag,
1981).

**38. **L. Jin, B. Liu, F. Zhao, S. Hahn, B. Dong, R. Song, T. C. Elston, Y. Xu, and K. M. Hahn, “Deep learning enables
structured illumination microscopy with low light levels and enhanced
speed,” Nat. Commun. **11**(1), 1934
(2020). [CrossRef]

**39. **V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling:
interpretable, efficient deep learning for signal and image
processing,” IEEE Signal Process.
Mag. **38**(2),
18–44 (2021). [CrossRef]