# PS2F: Polarized Spiral Point Spread Function for Single-Shot 3D Sensing

Bhargav Ghanekar, Vishwanath Saragadam, Dushyant Mehra, Anna-Karin Gustavsson,
Aswin C. Sankaranarayanan, Ashok Veeraraghavan B.G., V.S., A.V. are with the
Department of Electrical and Computer Engineering, Rice University, Houston,
TX.
E-mail<EMAIL_ADDRESS>D.M., A.-K.G. are with the Department of
Chemistry, Rice University, Houston TX. A.C.S. is with the Department of
Electrical and Computer Engineering, Carnegie Mellon University, Pittsburgh
PA.

###### Abstract

We propose a compact snapshot monocular depth estimation technique that relies
on an engineered point spread function (PSF). Traditional approaches used in
microscopic super-resolution imaging such as the Double-Helix PSF (DHPSF) are
ill-suited for scenes that are more complex than a sparse set of point light
sources. We show, using the Cramér-Rao lower bound, that separating the two
lobes of the DHPSF and thereby capturing two separate images leads to a
dramatic increase in depth accuracy. A special property of the phase mask used
for generating the DHPSF is that a separation of the phase mask into two
halves leads to a spatial separation of the two lobes. We leverage this
property to build a compact polarization-based optical setup, where we place
two orthogonal linear polarizers on each half of the DHPSF phase mask and then
capture the resulting image with a polarization-sensitive camera. Results from
simulations and a lab prototype demonstrate that our technique achieves up to
$50\%$ lower depth error compared to state-of-the-art designs including the
DHPSF and the Tetrapod PSF, with little to no loss in spatial resolution.

###### Index Terms:

Computational Photography, 3D Sensing, Microscopy, Phase Mask Design,
Polarization-encoded PSFs

##  1 Introduction

3D scanning is crucial to a wide range of applications, including microscopy
[1], autonomous driving [2], and robot-assisted surgeries [3]. Among the
multitude of approaches to measure depth, perhaps the hardest are those that
involve passive, monocular, and snapshot measurements—a scenario that
prioritizes compactness and time resolution. Estimating 3D information in this
context relies on cues such as shading [4] or defocus [5, 6]; however, the
underlying inverse problem is challenging and ill-posed.

A key property of real cameras is that the defocus blur changes with the depth
of scene points, which has been leveraged by prior work to obtain depth.
However, the blur produced by conventional pupils is not conducive to robust
depth estimation; to mitigate this, there has been significant interest in the
design of engineered pupil plane masks in the form of an amplitude [5], or
phase mask [7, 8, 9, 10, 11, 12, 13]. In particular, the seminal work of
Pavani et al. [9] has demonstrated that the so-called Double-Helix PSF
(DHPSF), consisting of two lobes rotating about a center, can provide high
spatial and depth resolutions, at least in the context of super-resolution
localization microscopy.

The majority of prior work on engineered PSFs focus on isolated point light
sources; this is a consequence of their intended application—namely super-
resolved localization of fluorescent particles—and the simplification that
this assumption provides. However, real-world scenes often consist of more
complex geometric primitives such as lines, edges, and curves [14, 15]; a PSF
optimized for point sources is inadequate at recovering such complex
geometries due to ambiguities in depth estimation. This is all the more
important in applications involving linear structures including blood
vasculature, neuronal networks (of the biological kind), and microtubules.

We propose and evaluate polarized spiral PSF (PS2F), a novel engineered PSF
that is well-suited for linear structures. Our key enabling observation is
that asymmetric PSFs achieve higher depth accuracy for scenes comprising of
linear structures. We achieve this asymmetry by observing that the DHPSF is
comprised of two lobes that rotate about a common center; by capturing images
with the individual lobes, we can reliably estimate depth. We propose a
compact realization of our approach with a novel polarization-based imaging
system. PS2F is generated by combining a DHPSF mask with orthogonal linear
polarizers on the two halves of the mask. We then use a polarization sensor,
capable of measuring images along four polarization angles in a snapshot by
using a Bayer-like tiling. Finally, the 3D scene is estimated by solving a
depth-dependent deconvolution problem. An overview of depth estimation with
PS2F is illustrated in Fig. 1.

Contributions.

We propose a new engineered PSF for imaging linear structures and make the
following contributions.

Figure 1: Polarized Spiral PSF (PS2F) for single-shot 3D scene estimation. We
propose a novel monocular depth imager based on an engineered PSF. Numerous
real-world scenes can be modeled as linear structures. We observe that
asymmetric PSFs, such as a single lobe of the DHPSF, enable us to accurately
estimate the 3D geometry of such linear structures. We leverage this
observation and propose a compact polarization-based setup that produces
accurate depth maps in a snapshot manner.

  * •

We theoretically demonstrate that the DHPSF is ill-suited for 3D imaging of
linear structures by analysing its Cramér-Rao lower bound (CRLB).

  * •

We propose a compact polarization-based setup, using a novel polarizer-phase
mask encoding that is capable of better 3D estimation in a snapshot manner.

  * •

We demonstrate the advantages of PS2F through several experiments with our lab
prototype over a wide range of scene geometries.

Limitations.

At its core, PS2F assumes that the scene intensity is unpolarized; if the
scene is polarized, the resultant depth estimate may not be accurate. Our
technique, like most techniques that rely on defocus blur, performs poorly for
regions with little or no textures or structures.

##  2 Related Work

PS2F draws motivation from the prior work on engineered PSFs. We discuss the
relevant ones in this section.

###  2.1 Mask-based PSF encoding

Encoding PSFs by designing masks has been used in several works, including
extended depth-of-field imaging [16, 17, 18], 3D sensing in the macroscopic
[6] and microscopic [7, 10] regimes, and in a lensless imaging setting [19].
The designed masks can be placed in front of the imaging lens or in the pupil
plane of the imaging system. Earlier works include designing amplitude masks
or coded apertures [5, 6] – such amplitude masks however lead to poor light
efficiency and very coarse depth information. Several follow up works utilized
phase-only modulation masks which include the rotating PSFs [7, 20, 21], and
PSFs obtained by optimizing the Fisher information [10]. Phase masks have also
been designed in an end-to-end optimization setting, both for macroscopic [22,
23] and microscopy [12, 13] tasks.

###  2.2 3D super-resolution microscopy using phase masks

Coupled with STORM [24] and PALM [25, 26], masks that result in rotating PSFs
with depth have been able to perform super-resolution microscopy [7, 20] of
individual fluorophore emitters to nanometer-scale accuracy. The most
successful of such designs is the DHPSF [7], a popular choice for 3D
localization of point sources [24, 27, 11]. 3D localization performance have
been further improved by designing Tetrapod PSFs [10], which were generated by
directly optimizing the Fisher information of the phase mask for $(x,y,z)$
localization of point sources. Recently, deep learning techniques have enabled
an end-to-end design framework that simultaneously learn phase masks as well
as neural networks for obtaining 3D localizations from raw captures [12, 13].

Of particular interest is the work by Nehme et al. [13], where a pair of masks
were designed in an end-to-end learning process. The two masks are placed in
two parallel $4f$ imaging channel. The mask pair, coupled with a learnt deep
network, allows for high accuracy 3D localization of single fluorophores over
a large depth range and for a high density of point-like emitters. Their work
focuses on super-resolution microscopy and on 3D localization of point-like
emitters, while our result tackles extended linear structures.

###  2.3 3D sensing of linear and extended structures

Human/animal bodies have a rich, dense 3D network of blood vessels running
through them. Capturing the 3D structure of vasculature is usually done using
variants of light-sheet microscopy (LSM) [28, 29, 30], confocal microscope
(CM) [31, 32], optical projection tomography [33], and optical coherence
tomographic angiography [34]. However, confocal and LSM methods involves
scanning individual points, lines, or planes which makes the systems bulky and
decreases the achievable imaging rate, making them non-ideal for fast dynamic
3D imaging.

###  2.4 Depth sensing using rotating PSFs

There is a rich body of work devoted to PSFs with lobes that rotate with
defocus. These can be clustered into two groups based on their use of Gauss-
Laguerre (GL) modes [9] or Fresnel rings [21].

GL-based rotating PSFs [7, 20] have been used for super-resolving and
localizing single molecules in 3D [7, 35, 36, 37]. Both GL-based [38] and
Fresnel-ring-based [39, 40] rotating PSFs have also been used for single-shot
depth estimation of scenes. However, the first work employs a two phase-mask
solution. The latter two employ a patch-based technique involving processing
in the cepstrum domain. They only recover very simple scenes at low spatial
and depth resolutions. Such recovery approaches do not extend to more complex
geometry such as linear structures.

###  2.5 Using polarization channels for PSF encoding

There have been a few works that leverage polarization to help in 3D super-
resolution imaging [41, 13, 42]. These works typically use a beam-splitter to
split the incoming light into two separate $4f$ optical system channels which
enabled individual channels to be encoded with a different mask. The images
were then captured over two non-overlapping regions of a single sensor, or on
two separate sensors. Such solutions however call for bulkier optics, and more
importantly, a need for sub-pixel alignment to achieve high spatial
resolution. We instead leverage snapshot polarimetric cameras that are
naturally well-aligned, and hence easier to work with, and more compact.

##  3 Polarized Spiral Point Spread Function

We begin with a brief background on the DHPSF, an analysis of its properties,
and its performance for linear/extended structures. Subsequently, we introduce
our proposed PSF and its associated reconstruction technique.

###  3.1 Background: DHPSF and mask design

This subsection outlines the background on the generation of the DHPSF [7] and
other GL-based rotating PSFs [20].

####  3.1.1 Phase mask for GL-based rotating PSFs

Any propagating paraxial wave/beam can be expressed in terms of the orthogonal
basis of GL modes, which are characterized by two integers $(n,m)$. Different
combinations of the GL modes lead to varying properties of paraxial beams. In
[43], it was shown that beams that constituted of GL modes that are in
arithmetic progression continuously rotated and scaled upon propagation. Using
a beam corresponding to GL modes of (1,1), (5,3), (9,5), (13,7), it was shown
in [8] that the amplitude and phase profile of the beam at $z=0$ can be used
as a mask in front of a lens (or in the Fourier plane of a $4f$ lens system)
to generate a PSF that rotates with defocus. Subsequently, this mask was
implemented in a phase-only manner using a variant of Gerchberg-Saxton
optimization procedure with constraints included for the GL modal domain [9].
The resulting phase function creates a PSF that rotates with changing point
source depth, and is popularly known as the DHPSF due to its double helical
structure. This PSF has been successfully used in 3D localization of
fluorescent particles [7][35].

####  3.1.2 Properties of GL-based rotating PSFs

Piestun et al. [43] present a theory on paraxial rotating beams using GL
modes, and characterizing the resulting rotation range, rotation rates, and
beam scaling rates. Specifically, given a $4f$ system with lenses of focal
length $f$, a GL-based mask designed with beam width $w_{0}$ for wavelength
$\lambda$, and having GL modes lying on a single line with slope $V_{1}$, the
following can be established:

  * •

The total rotation in one direction of defocus is $(V_{1}\pi/2)$

  * •

The angle of the rotating PSF

$\phi(z)=\phi_{0}+V_{1}\ \text{tan}^{-1}\left(\frac{z}{\lambda f^{2}/\pi
w_{0}^{2}}\right)$

  * •

The PSF will rotate by $(V_{1}\pi/4)$ over a depth of $\frac{\lambda
f^{2}}{\pi w_{0}^{2}}$ .

Even though phase-only masks for rotating PSFs are generated after an
optimization procedure [9], the above rates and ranges are still a good
approximation to understand the properties of any GL-based rotating PSF.

###  3.2 Challenges in imaging linear/extended structures

Phase masks such as the DHPSF [7], or the Tetrapod PSF [10], were designed for
3D estimation of point sources. To image scenes with extreme fluorophore
density, the concept of a scene made up of point-like emitters is not
accurate. With increasing fluorophore/light source density, scenes consist of
edge-like and linear structures. Using a DHPSF for recovering depth of a line
then leads to an ambiguity between depth and orientation of the line. This is
illustrated in Fig. 2. Given a line image formed using the DHPSF, there will
always be a global depth ambiguity—two different depths can give rise to the
same line image. This makes the DHPSF ill-suited for imaging edge/linear
structures.

Figure 2: Ambiguities caused by DHPSF for lines. Estimation of depth and
orientation of a line using the DHPSF is an ill-conditioned problem. There is
an inherent global ambiguity in depth estimation, since for every line image,
two possible candidate depths are possible.

###  3.3 Designing a PSF for 3D lines

A line with constant intensity blurred with an arbitrary PSF produces another
line – thereby leading to loss of information along the line direction. In
case of the DHPSF, the loss of information implies that two lines at two
different depths ($-\Delta z$ and $+\Delta z$) produce the same measurement.
This ambiguity can however be resolved if we measure two separate images with
each lobe of the PSF. Imaging with the individual DHPSF lobes separately will
remove the global degeneracy issues (up to a $2\pi$ lobe rotation range).
Next, we provide a theoretical analysis of the advantages to be obtained by
such a separation. For estimation from data (image(s) in our case), we can use
theoretical tools to understand what ultimate precision is possible. The
Fisher information and Cramér-Rao Lower bound (CRLB) give the best possible
precision that can exist by any estimator. In the case of an image of an
isolated point source, the Fisher information matrix for estimating parameters
$\theta=(x,y,z)$ is given by [10]:

$[FI(\theta)]_{ij}=\sum_{k=1}^{N_{p}}\frac{1}{\mu_{\theta}(k)+\beta}\left(\frac{\partial\mu_{\theta}}{\partial\theta}\right)_{i}\left(\frac{\partial\mu_{\theta}}{\partial\theta}\right)^{T}_{j},$
(1)

where $N_{p}$ is the number of pixels, $k$ is a variable for summing over all
pixels, $T$ corresponds to the transpose operator, $\beta$ is the background
Poisson noise level, and $\mu_{\theta}$ is the scaled PSF image for a point
source at $\theta=(x,y,z)$ (scaled by the number of photons $N$). The diagonal
elements of the inverse of the Fisher information matrix provides the CRLB for
estimation of each of the parameters $\theta=(x,y,z)$.

$[\textrm{CRLB}(\theta)]_{i}\equiv[FI(\theta)]^{-1}_{ii}$

We next extend the Cramér-Rao analysis to imaging of lines. Given a line image
(or a one-sided edge image), the parameters to be estimated are the associated
depth $z$ of the line/edge (at center of the patch), the associated
orientation $\phi$, i.e. $\theta=(z,\phi)$. Thus, given a line image $\psi$,
the CRLB estimates for depth $z$ is given by

$\displaystyle[FI(\theta)]_{ij}=\sum_{k=1}^{N_{p}}\frac{1}{\psi_{\theta}(k)+\beta}\left(\frac{\partial\psi_{\theta}}{\partial\theta}\right)_{i}\left(\frac{\partial\psi_{\theta}}{\partial\theta}\right)^{T}_{j}$
(2) $\displaystyle[\textrm{CRLB}(\theta)]_{z}\equiv[FI(\theta)]^{-1}_{zz}$ (3)

The estimation of line orientation from a line path is easier than depth
estimation, as orientation can be computed from the global structure of the
line image. Hence, for simplicity of exposition, we focus on the analysis of
the CRLBz values here. Analysis regarding the CRLBϕ values is shown in the
Supplementary. Using Eqns 2, 3 we can calculate the $\sqrt{\textrm{CRLB}_{z}}$
values for a given line image for all $(z,\phi)$. This calculation can be
readily extended to line images captured using two separate PSFs. In such a
case, the final Fisher information matrix is the sum of the individual Fisher
matrices. Assuming a $4f$ system with 50 mm focal length lenses, and
$N=100,000$ photons, $\beta=5$ photons/pixel, we compute the
$\sqrt{\textrm{CRLB}_{z}}$ plots (as a function of $z$, $\phi$) for the DHPSF
[9], the Tetrapod PSF [10], and the individual DHPSF lobes with $25\%$ photons
in each lobe. Note that $\sqrt{\textrm{CRLB}_{z}}$ obtained with a mask pair
should be compared with a single mask having 2x SNR (2x total photons), as
shown in [13]. However, our polarizer-phase mask setup causes a $50\%$ light
loss (see further in Section 3). Hence, for an appropriate comparison, we
compare the individual DHPSF lobe pair, with other PSFs having $4\times$ the
SNR. There are possible ways to remove this loss, with a setup demonstrated in
[13], which will further lower the estimated $\sqrt{\textrm{CRLB}_{z}}$
values.

Fig. 3 illustrates the $\sqrt{\textrm{CRLB}_{z}}$ values (log-scale) for a
given line at depth $z$ and having orientation $\phi$. The mean and standard
deviation of the $\sqrt{\textrm{CRLB}_{z}}$ values are also shown in the
individual insets. The DHPSF plot in Fig. 3 shows several peaks with large
magnitude. This occurs at specific $(z,\phi)$ where the line connecting the
two lobes of the PSF and the scene line are either perpendicular or parallel
to each other. The ability to discern lines at slightly differing depths is
greatly reduced at these points due to symmetry. For the Tetrapod PSF [10] and
the PhaseCam3D PSF [22], the plots in Fig. 3 also show more peaks and ridges,
but with lower magnitude compared to the DHPSF. For the PSF pair of the
individual DHPSF lobes, we see that most of the peaks are removed, except for
a central peak in the plot. The peak corresponds to the case when the PSF is
in focus and the the line orientation is perpendicular to the PSF lobes. This
central peak does not affect the effectiveness of separating lobes, because it
is a narrow peak over a only small area in the parameter space of $(z,\phi)$.
The average $\sqrt{\textrm{CRLB}_{z}}$ values for the individual DHPSF lobe
pair is similar to the DHPSF with 4x SNR. Thus, with 4x lower SNR, the PS2F is
able to achieve a comparable $\sqrt{\textrm{CRLB}_{z}}$ value with other PSFs,
and does not show consistent peaks and ridges as seen in the plots. Note that
CRLB calculations only provide insight about local ambiguities or precision
levels. However, we also observe that a PSF pair created out of separating the
DHPSF lobes removes global ambiguities as well. Thus overall, a PSF pair
created from the individual DHPSF lobes is better for estimating depths of
line/edge patches.

Figure 3: CRLB comparison for various phase masks. The figure shows
$\sqrt{CRLB_{z}}$ log-intensity (log-mm) plots for line images, as a function
of line depth $(z)$ and line orientation $(\theta)$. PS2F achieves a
comparable CRLB over the whole range of depth and orientations with fewer
peaks, even with 4x lower SNR. In contrast, competing approaches have curves
with large CRLB values, which introduce ambiguities in the depth estimate.

###  3.4 Using polarizers to separate out the DHPSF lobes in a single $4f$
optical system

Imaging with a PSF pair is possible by having two separate parallel $4f$
imaging systems with separate sensors, as demonstrated in [13]. However, the
DHPSF mask (and the Fresnel ring-based two-lobe masks) enjoy the property of
PSF lobe separability. Partitioning such a mask into two halves along a
particular axis, and allowing light through only one half leads to the
creation of only a single rotating PSF lobe, as depicted in Fig. 4. The
primary reason is that the light falling on one half is modulated to form one
lobe, while the light falling on the other half is modulated to form the other
lobe. Such a partitioning leads to partitioning of the DHPSF into two
separate, distinct lobes.

Figure 4: Effect of separating the phase mask. The phase mask for generating
the DHPSF has a special property that separating the phase profile along the
highlighted red line yields a separation of the two lobes. We leverage this to
implement a compact realization of Polarized Spiral PSF (PS2F) using a
polarization camera sensor and two polarizers in the pupil plane.

The separability of DHPSF lobes can be leveraged to build a compact imaging
system. In PS2F, we add two linear polarizers on each mask half, one oriented
along $s$-polarization and the other oriented along $p$-polarization. This
ensures that the $s$-polarization component is modulated using one mask half,
and consequently the $p$-polarized component is modulated using the other mask
half. The PSF created out of such a polarizer-phase mask encoding contains two
lobes, but with each lobe being in a different polarization (either $s$ or
$p$). This addition of polarizers causes a $50\%$ light loss, which was taken
into consideration while estimating the $\sqrt{\textrm{CRLB}_{z}}$ values.
Such a polarized spiral PSF can be captured efficiently using a polarized
camera sensor. A polarized camera sensor contains a 2D array of pixels with a
Bayer pattern that has 90, 45, 135, 0 degree polarizers. After addition of the
polarizers in the mask-plane, the orthogonally polarized DHPSF lobes can be
imaged in their respective channels of a polarization sensor. This allows for
a compact $4f$ system, allowing the imaging of a PSF pair using just a single
$4f$ system channel, making the proposed method a single-shot, single-sensor
method.

###  3.5 Imaging model and reconstruction procedure

We assume that the light from the scene is unpolarized in our imaging model.
The imaging of a 3D scene with a depth-dependent PSF can be approximated as a
sum of 2D convolutions between the depth-dependent PSF and the per-plane scene
intensity:

$I_{c}(x,y)=\sum_{z=z_{start}}^{z_{end}}h_{c}(x,y;z)*s(x,y,z)$ (4)

where $I_{c}(x,y)$ is the image intensity at $(x,y)$ in polarization channel
$c$, $h_{c}(x,y;z)$ is the 2D PSF corresponding to depth $z$ and polarization
channel $c$, and $s(x,y,z)$ is the scene intensity at point $(x,y,z)$.

The goal of 3D reconstruction is to estimate a 3D matrix $s(x,y,z)$ from the
noisy measurements of $I(x,y)$. There are numerous ways to solve this problem
including regularized least-squares [19, 44, 45], data learning-based [22,
23], and more recently, ones based on deep network-based regularizer [46]. To
keep the exposition simple, we demonstrate recovery by modelling the 3D
reconstruction problem as a regularized least squares optimization problem.

Specifically, we formulate the 3D scene estimation problem as:

$\underset{\textbf{x}}{\text{argmin}}\left\lVert\begin{bmatrix}I_{0}\\\
I_{90}\end{bmatrix}-S\begin{bmatrix}H_{0}\\\
H_{90}\end{bmatrix}\textbf{x}\right\rVert^{2}_{2}+\lambda_{TV}\lVert\Psi(\textbf{x})\rVert_{1}+\lambda_{L1}\lVert\textbf{x}\rVert_{1}$
(5)

where x is a 3D matrix of scene intensities, $H_{0},H_{90}$ are the depth-
specific PSF operators corresponding to the individual two lobes that are
polarized to $0^{\circ}$ and $90^{\circ}$ states respectively, and
$I_{0},I_{90}$ are captured images in the $0^{\circ}$ and $90^{\circ}$
polarization channels respectively. $S$ is the summing operator that sums
across the depth channels of $H_{0}\textbf{x}$ and $H_{90}\textbf{x}$
separately. We employ TV and L1 regularizers as scene priors, with
$\lambda_{TV},\lambda_{L1}$ as hyperparameters to control their regularization
effects respectively. We solve the optimization problem using autograd
functionality in PyTorch and the Adam optimizer [47] to leverage the speed of
graphical processing units (GPUs).

A key requirement for our technique is that the scene emits unpolarized light.
This is required so as to ensure that both the lobes have produce PSFs with
similar intensity levels that we can calibrate a priori. When the incident
light is polarized, the blur kernels in Eq. 4 will have an unknown scaling,
which leads to a model mismatch; in an extreme scenario, we lose all the
information in one of the lobes if the polarization angle of the incident
light is orthogonal to the corresponding polarizer in the pupil plane. While
this is a limitation of our technique, we largely encountered unpolarized
light in our intended application of fluorescence microscopy.

##  4 Simulations

Figure 5: Simulations on vascusynth dataset with varying complexity. The
figure shows reconstructions with various PSF designs on scenes with medium
complexity (top row) and high complexity (bottom row). PS2F is superior in
depth estimation, especially in regions where there is higher vasculature
(scene) complexity (see insets). Scale bars indicate $0.25$ mm.

To evaluate the performance of our proposed PS2F, we perform extensive
simulations on the VascuSynth 2013 dataset [48, 49], which provides 10
simulated 3D volumes of vascular trees for each of 12 levels of complexity.
Scene complexity is set by the number of bifurcations in the vascular tree
(ranging from 1 to 56 in steps of 5).

Figure 6: Accuracy vs. scene complexity. Depth estimation performance, for the
proposed PS2F and other comparing PSFs (averaged over 10 scenes for each
complexity level, shaded regions show variance across the 10 scenes). This
demonstrates the suitability of PS2F, and that the separation of the DHPSF
lobes helps in 3D reconstruction.

Each 3D vasculature volume in the dataset consists of $101\times 101\times
101$ voxels. These volumes were converted into 3D scenes of size $256\times
256\times 256$ through interpolation methods. This new volume was assumed to
cover a volume of $1.76\times 1.76\times 5.00$ mm3 with each voxel occupying a
$6.9\times 6.9\times 19.5~{}\mu$m3.

We further assume a $4f$ system with both lenses having focal lengths $f=50$
mm. We render images of the vascular structures using several different masks
of diameter $3$ mm in the Fourier plane. We assume that the light from the
vasculature is monochromatic of wavelength $\lambda=532$ nm. We render images
assuming occlusion of vasculature and that light is received only from the
surface of the vasculature. We assume that the center of the 3D scene is
$f=50$ mm away from the first lens, thus, the scene spans a defocus range of
[$-2.5$ mm, $+2.5$ mm]. The masks (and their corresponding PSFs) used in
rendering the simulated images were:

  * •

DHPSF: We optimize for a GL-based DHPSF mask (following procedure in [9]),
based on the rotating paraxial beam with $w_{0}=0.4$ mm and having GL modes
corresponding to $(1,1),(3,5),(5,9),(7,13)$. Full 180 degree lobe rotation was
achievable over $2\frac{\lambda f^{2}}{\pi w_{0}^{2}}=5.3$ mm.

  * •

PS2F: Using the above mask, we construct an equivalent mask pair corresponding
to PS2F by partitioning the mask into two halves appropriately.

  * •

Tetrapod PSF: We obtain a Tetrapod PSF mask designed for $550$ nm wavelength
from [36] and repurpose it for imaging over our specified $5$ mm depth range.

We then added Poisson and Gaussian read-out noise to the final measurements.
For an appropriate comparison, we perform 3D reconstructions with individual
single lobes (PS2F), two-shot capture with DHPSF, and a two-shot capture using
the Tetrapod PSF. To accurately account for the $50\%$ light loss in the PS2F
case, we also render out images with half the signal level as compared to the
signal levels in DHPSF and Tetrapod PSF case. After obtaining the 3D
volumetric estimate, we filter out points whose sum across the z-stack is
lower than a fixed threshold. We then estimate depth using the index
corresponding to a maximum-intensity-projection (MIP) of the 3D estimate.

###  4.1 Comparison with DHPSF and Tetrapod PSF

Fig. 6 shows the depth estimation performance for the proposed PS2F, DHPSF,
and the Tetrapod PSF. We compare the PSFs depth estimation using the Mean
Absolute Error (MAE), Root Mean Square Error (RMSE), and the multi-scale
Structural Similarity Index Measure [50] (MS-SSIM) metrics. The entire set of
results can be seen in the Supplementary. The proposed PS2F method performs
$\sim 2\times$ better in terms of the RMS error for depth estimation. The
depth estimation performance improvement is even higher ($\sim 3\times$) for
vasculature with more bifurcations, i.e. greater complexity. Depth estimation
results for two example vasculature scenes can be seen in Fig. 5. In Fig. 5
(top row), the vasculature scene is not very complex. The depth map estimation
results are fairly similar for all the three PSFs (PS2F, DHPSF, Tetrapod PSF),
but the zoom-in insets highlight the errors seen in the DHPSF and Tetrapod PSF
depth estimation. In Fig. 5 (bottom row), the scene has greater complexity,
and the PS2F produces much better depth estimation than the other two PSFs.
The zoom-in insets especially highlight the same.

Figure 7: Simulated reconstruction results of a USAF target at a skew angle.
(A) Shows the depth map estimates for PS2F, DHPSF, and Tetrapod PSF
respectively. The XZ maximum intensity projections are shown in (B) in their
respective columns, with ground truth shown in the inset. (C) shows zoom-in
insets of depth maps (after thresholding to remove small artifacts) of
highlighted regions in (A). PS2F has a better reconstruction in lateral and
axial dimensions as compared to the DHPSF and the Tetrapod PSF. Scale bars for
(A), (B), (C) indicate $0.5$ mm, $0.5$ mm, and $0.1$ mm respectively.

###  4.2 Resolution performance

Furthermore, we analyse the reconstruction resolution performance by
reconstructing a simulated skew USAF target. We assume a USAF target places at
the focus plane, but tilted with respect to the optical axis so that its depth
changing horizontally across the target; we render using the same three masks
in the same $4f$ optical system as above. The results are shown in Fig. 7. The
proposed PS2F generates a better reconstruction than the DHPSF and Tetrapod
PSF reconstructions, both in lateral and axial dimensions. As seen from the XZ
MIP plots in Fig. 7, the PS2F and Tetrapod PSF reconstructions are able to
obtain the skew angle of the target, with the former being much better. The
DHPSF fails to get the skew angle, due to the global depth ambiguity issues
explained in Section 3. In the lateral dimensions, the Tetrapod PSF is unable
to reconstruct Group 4 and 5. This could be attributed to the fact that the
PSF has a larger support. On the other hand, the PS2F and DHPSF have
concentrated lobes, and thus are able to reconstruct elements of Group (5,3)
(12.40 $\mu$m linewidth) and (5,2) (13.92 $\mu$m linewidth) respectively.

The obtained $xy$-resolution can be theoretically justified. The separation of
lobes in the PS2F method allows for resolution of finer scene features due to
the compactness of the PSF. The resolution will be roughly determined by the
width of the individual PSF lobes. Fitting a 2D Gaussian to each individual
lobe, the average $2\sigma$-value for the Gaussian fit is $\sim
1.74~{}\text{pixels}=12~{}\mu$m in the object space. This readily matches with
the fact that the PS2F can resolve elements of Group (5,3), which has line
elements that are spaced at $12.40~{}\mu$m distance. The diffraction-limited
resolution of an Airy disk PSF for the same optical system configuration is
given by $1.22\lambda f/D=10.81~{}\mu$m, indicating that the $xy$-resolution
obtained by the PS2F is close to the diffraction-limited resolution
($1.15\times$) as well.

For PS2F, we also analysed the spread of reconstruction signal across the
z-channel. For pixels with a significant signal level, we fit a 1D Gaussian
curve to estimate this spread. The median $2\sigma$-value across z-channel
obtained was $2\times 304=608~{}\mu$m, which is an estimate of the axial
resolution performance of the proposed PS2F. Note that the diffraction-limited
axial performance of an Airy disk PSF for the same optical system
configuration is given by $4\lambda(f/D)^{2}=591~{}\mu$m. We are able to
achieve $\sim 1.03\times$ the diffraction-limited axial resolution with the
proposed PS2F.

##  5 Experimental Results

Figure 8: Lab prototype. (A) Our lab prototype is a $4f$ system consists of a
polarization sensor (FLIR BFS-U3-51S5P-C) equipped with a $50$ mm F-mount
lens. The phase mask and polarizers are on an XY stage. An achromat
(AC254-050-A) forms the other half of the $4f$ system. We used an additional
z-axis motorized translation stage to capture the PSF stack. (B) Two
orthogonal polarizer pieces. (C) Images of the polarizers and phase mask with
a camera focused on the pupil-plane.

We perform real-world imaging with PS2F to verify our simulations. In this
section we detail out the experiments and compare the PS2F with the
corresponding DHPSF.

###  5.1 Imaging setup description

We built a $4f$ imaging setup with the first lens being an achromatic lens
with focal length $50$ mm, and the second lens being a Canon camera lens with
effective focal length $=50$ mm. The polarizer-phase mask is added to the
pupil plane of the $4f$ system. Fig. 8 shows the prototype.

####  5.1.1 Polarizer-Mask design

We fabricated a $3$ mm diameter DHPSF mask using a 3-D printer employing two-
photon lithography. See Supplementary for more details on fabrication. The
fabricated mask was the same mask used in the Section 4 to perform simulations
of the VascuSynth dataset. To create a single-mask, single-sensor design for
the PS2F as outlined in Section 3, two rectangular pieces of thin-film
polarizers were laser cut and put side-by-side in a 3D printed holder. Both
the polarizer pieces were laser cut such that their polarization orientation
was $0\deg$ and $90\deg$ respectively. The fabricated DHPSF mask was attached
to the back side of the 3D printed holder, carefully aligning the two mask
halves w.r.t. to the polarizer edges under a bright-field 4x microscope. This
polarizer-mask aperture is illustrated in Fig. 8(B),(C). Further details can
be found in the Supplementary figures.

####  5.1.2 PSF capture and calibration

To reconstruct 3D scenes, we use an experimentally captured PSF stack for the
PS2F. Using the polarization sensor, we image a $5~{}\mu$m pinhole over a
range of [$-2.5$ mm, $+2.5$ mm] on either side of the focus plane. The
experimentally-captured PSFs can be seen in Fig. 9.

###  5.2 Imaging experiments

With the fabricated mask and optical setup, we perform a series of experiments
with 3D skew planar targets (using transmissive illumination) and with 3D
fluorescence targets. We image using a polarization camera sensor (FLIR
BFS-U3-51S5P-C) to obtain the individual PS2F lobes in two polarization
channels. By averaging across the Bayer pattern, we can readily obtain the
image of the scene as imaged by a DHPSF, which we use for comparison purposes.

Figure 9: PSF stacks captured with a lab prototype. We achieved high quality
separation between the two lobes in 0 and 90 degree images with our setup. We
noticed no significant effect of any misalignment between the phase mask and
the two polarizer halves. Scale bars indicate $100~{}\mu$m. Figure 10: Skew
USAF target reconstruction. Showing reconstructed depth map, YZ maximum
intensity projection, zoom-in inset of highlighted ROI for (A) PS2F, and for
(B) DHPSF. Column (C) shows slice plots for marked yellow lines in the zoom-in
insets, corresponding to Group (5,3) elements for PS2F (top plot) and DHPSF
(bottom plot). We were able to resolve up to Group (5,3) elements (12.40
$\mu$m linewidth) using the PS2F, whereas only up to Group (4,5) elements
(19.69 $\mu$m linewidth) using the DHPSF. Text below YZ MIP plots indicate the
depth RMS error obtained for reconstruction when compared to ground truth,
showing that the PS2F performs $\sim 2\times$ better. All scale bars indicate
$1$ mm. Figure 11: 3D reconstructed fluorescent bead volumes. The scatter
plots show 3D reconstructed views with PS2F and the DHPSF for bead densities
approximately corresponding to (A) 3 beads/mm3, (B) 5 beads/mm3, and (C) 8
beads/mm3. We captured ground truth information with a semi-aligned confocal
ground truth (hence an approximate ground truth). Highlighted regions indicate
various errors: incorrect estimation by DHPSF (solid green), missing
estimation by DHPSF (dashed green), incorrect estimation by PS2F (solid blue),
and incorrect estimation by both PSFs (solid yellow). Axes values are in units
of $mm$. Figure 12: Experimental results on fluorescent cotton strands. (A)
PS2F and (B) DHPSF. First column in (A), (B) show the captured scene images
(polarization channels in the PS2F image are shown in red and green colors).
Second column in (A), (B) shows the reconstructions. Highlight ROIs (in
yellow) are shown in the last column of (A), (B) – depicting the inaccuracy of
DHPSF, as well as the ability of PS2F to reconstruct scenes with complex
linear structures. In the last row (C), columns (from left-to-right) show
color-coded depth maps obtained from confocal imaging (treated as ground
truth), PS2F reconstructions, and DHPSF reconstructions respectively. Note the
confocal GT depth map is an approximate depth map, due to imperfect
registration between GT and the imaged volumes.

####  5.2.1 3D skew planar targets

We imaged a USAF planar target placed at a skew angle, and then performed
single-shot 3D reconstruction of the target. Furthermore, we obtain a focal
stack of the target so as to obtain ground truth depth map. The ground truth
depth map was obtained from the focal stack using the Shape-from-focus method
[51]. For the FOV being imaged, the relative intensity of each of the PS2F
lobes changed across the FOV, causing a spatial variance in the PSF. This was
readily accounted for by adding a per-pixel, per-channel weight term to the 3D
scene variable x in Eq. 5, and jointly optimizing for the 3D scene x and the
weights. The results are depicted in Fig. 10. After modelling of the relative
strengths of the PS2F lobes as variables to be estimated during the
optimization, we achieved an accurate reconstruction result matching the
simulation results in terms of resolving Group (5,3) elements for the PS2F.
The YZ maximum intensity projections demonstrate that the skew plane of the
USAF has been correctly estimated using the PS2F, as compared to using the
DHPSF. When judged against the ground truth, we see a $\sim 50\%$ lower RMS
error achieved with the PS2F as compared to the DHPSF.

####  5.2.2 3D fluorescence imaging

We, next, imaged a set of prepared fluorescent samples, consisting of
$10~{}\mu$m fluorescent beads in a PDMS substrate (for more information, see
Supplementary). The bead concentration was varied from approximately $3$ to
$8$ beads/mm3. The resulting samples were $\sim 2$ mm in thickness. We imaged
these samples using the same mask as before, but in a $4f$ optical system
using lenses of focal length $30$ mm and $50$ mm resulting in a $1.6\times$
magnification. We reconstructed using the same optimization algorithm as
specified in Eq. 5, and compared against ground truth data obtained with a
confocal microscope. The results are shown in Fig. 11. Bead localization
performance using the DHPSF seems to get relatively worse as the bead
concentration increases. In contrast, PS2F enabled accurate results with
increasing bead concentration. It is important to note that in this
experiment, we use an optimization objective to perform bead reconstruction.
More sophisticated and accurate methods (such as likelihood-based methods)
exist for reconstruction of beads or point-like objects. Moreover, the ground
truth (GT) data was captured using confocal imaging method. However, due to
the complexity of matching/registering imaging volumes, we only obtain an
approximated alignment or approximate registration between the GT and the
imaging volume captured.

Using the same optical setup as above, we also image fluorescent cotton
strands suspended in a PDMS sample (see Supplementary for more information
about sample preparation). Fig. 12 shows the results for the same. PS2F method
is better able to capture the 3D structure of the fluorescent strands -
especially in parts with highly complex geometry. Imaging with PS2F results in
sharper reconstructions, which helps in reconstructing the fine, complex
structure of the strands. The estimated depth map from the PS2F is also more
accurate as compared to the DHPSF case.

##  6 Conclusion

We presented PS2F for single shot, monocular depth estimation of extended
(linear) structures where we separated the two lobes of the DHPSF into the two
orthogonal states of polarization. This separation enabled us to remove
ambiguities when estimating depth of line segments by breaking the inherent
symmetry of DHPSF with respect to depth. The snapshot capabilities of PS2F
will enable faster microscopic imaging, including high resolution light sheet
microscopy at real-time rates with fewer captures. We designed and
demonstrated a compact physical realization of PS2F with a single,
polarization-sensitive camera, and showed that PS2F results in $2\times$ or
higher accuracy compared to current state-of-the-art phase masks such as DHPSF
and the Tetrapod PSF. We believe that our approach of leveraging polarization
multiplexed phase masks combined with polarization-sensitive cameras opens new
avenues for high resolution snapshot depth imaging at micro and macro scales.

PS2F is inherently designed for extended/linear structures, and imaging such
structures will involve high photon counts than when imaging single molecules
or point-like emitters. Hence, we do not specifically deal with a low-photon
count scenario in this work. Nevertheless, PS2F can be adapted to low-photon
count scenarios by using advanced recovery techniques including deep network-
based approaches. Improvements such as modeling spatially varying PSF, as well
as accounting for polarized input light will further improve PS2F results.

## Acknowledgments

The authors would like to thank Dong Yan for his help in preparing the
fluorescent bead and fluorescent strands samples. This work was supported by
NSF awards IIS-1730574, IIS-1730147, CCF-1652569, IIS-1652633, and
EEC-1648451. A.-K.G. acknowledges partial financial support from the National
Institute of General Medical Sciences of the National Institutes of Health
(Grant No. R00GM134187), the Welch Foundation (Grant No. C-2064-20210327), and
startup funds from the Cancer Prevention and Research Institute of Texas
(Grant No. RR200025).

## References

  * Fischer et al. [2011] R. S. Fischer, Y. Wu, P. Kanchanawong, H. Shroff, and C. M. Waterman, “Microscopy in 3d: a biologist’s toolbox,” _Trends in Cell Biology_ , vol. 21, no. 12, pp. 682–691, 2011.
  * Arnold et al. [2019] E. Arnold, O. Y. Al-Jarrah, M. Dianati, S. Fallah, D. Oxtoby, and A. Mouzakitis, “A survey on 3d object detection methods for autonomous driving applications,” _IEEE Trans. Intelligent Transportation Systems_ , vol. 20, no. 10, pp. 3782–3795, 2019.
  * Reiter et al. [2014] A. Reiter, A. Sigaras, D. Fowler, and P. K. Allen, “Surgical structured light for 3d minimally invasive surgical imaging,” in _IEEE Intl. Conf. Intelligent Robots and Systems_ , 2014.
  * Ping-Sing and Shah [1994] T. Ping-Sing and M. Shah, “Shape from shading using linear approximation,” _Image and Vision Computing_ , vol. 12, no. 8, pp. 487–498, 1994.
  * Levin et al. [2007] A. Levin, R. Fergus, F. Durand, and W. T. Freeman, “Image and depth from a conventional camera with a coded aperture,” _ACM Trans. Graphics_ , vol. 26, no. 3, pp. 70–es, 2007.
  * Zhou et al. [2009] C. Zhou, S. Lin, and S. Nayar, “Coded aperture pairs for depth from defocus,” in _IEEE Intl. Conf. Computer Vision (ICCV)_ , 2009.
  * Pavani et al. [2009] S. R. P. Pavani, A. Greengard, and R. Piestun, “Three-dimensional localization with nanometer accuracy using a detector-limited double-helix point spread function system,” _Appl. Physics Letters_ , vol. 95, no. 2, p. 021103, 2009\.
  * Greengard et al. [2006] A. Greengard, Y. Y. Schechner, and R. Piestun, “Depth from diffracted rotation,” _Optics Letters_ , vol. 31, no. 2, pp. 181–183, 2006.
  * Pavani and Piestun [2008] S. R. P. Pavani and R. Piestun, “High-efficiency rotating point spread functions,” _Optics Express_ , vol. 16, no. 5, pp. 3484–3489, 2008.
  * Shechtman et al. [2014] Y. Shechtman, S. J. Sahl, A. S. Backer, and W. E. Moerner, “Optimal point spread function design for 3d imaging,” _Physical Review Letters_ , vol. 113, no. 13, p. 133902, 2014.
  * Badieirostami et al. [2010] M. Badieirostami, M. D. Lew, M. A. Thompson, and W. Moerner, “Three-dimensional localization precision of the double-helix point spread function versus astigmatism and biplane,” _Appl. Physics Letters_ , vol. 97, no. 16, p. 161103, 2010.
  * Nehme et al. [2020a] E. Nehme, D. Freedman, R. Gordon, B. Ferdman, L. E. Weiss, O. Alalouf, T. Naor, R. Orange, T. Michaeli, and Y. Shechtman, “Deepstorm3d: dense 3d localization microscopy and psf design by deep learning,” _Nature Methods_ , vol. 17, no. 7, pp. 734–740, 2020.
  * Nehme et al. [2020b] E. Nehme, B. Ferdman, L. E. Weiss, T. Naor, D. Freedman, T. Michaeli, and Y. Shechtman, “Learning an optimal psf-pair for ultra-dense 3d localization microscopy,” _arXiv preprint arXiv:2009.14303_ , 2020.
  * Cho et al. [2011] T. S. Cho, S. Paris, B. K. Horn, and W. T. Freeman, “Blur kernel estimation using the radon transform,” in _IEEE Comp. Vision and Pattern Recognition (CVPR)_ , 2011.
  * Joshi et al. [2008] N. Joshi, R. Szeliski, and D. J. Kriegman, “Psf estimation using sharp edge prediction,” in _IEEE Comp. Vision and Pattern Recognition (CVPR)_ , 2008\.
  * Castro and Ojeda-Castañeda [2004] A. Castro and J. Ojeda-Castañeda, “Asymmetric phase masks for extended depth of field,” _Appl. Optics_ , vol. 43, no. 17, pp. 3474–3479, 2004.
  * Yang et al. [2007] Q. Yang, L. Liu, and J. Sun, “Optimized phase pupil masks for extended depth of field,” _Optics Communications_ , vol. 272, no. 1, pp. 56–66, 2007.
  * Zammit et al. [2014] P. Zammit, A. R. Harvey, and G. Carles, “Extended depth-of-field imaging and ranging in a snapshot,” _Optica_ , vol. 1, no. 4, pp. 209–216, 2014.
  * Boominathan et al. [2020] V. Boominathan, J. K. Adams, J. T. Robinson, and A. Veeraraghavan, “Phlatcam: Designed phase-mask based thin lensless camera,” _IEEE Trans. Pattern Analysis and Machine Intelligence_ , vol. 42, no. 7, pp. 1618–1629, 2020.
  * Lew et al. [2011] M. D. Lew, S. F. Lee, M. Badieirostami, and W. Moerner, “Corkscrew point spread function for far-field three-dimensional nanoscale localization of pointlike objects,” _Optics Letters_ , vol. 36, no. 2, pp. 202–204, 2011\.
  * Prasad [2013] S. Prasad, “Rotating point spread function via pupil-phase engineering,” _Optics Letters_ , vol. 38, no. 4, pp. 585–587, 2013.
  * Wu et al. [2019] Y. Wu, V. Boominathan, H. Chen, A. Sankaranarayanan, and A. Veeraraghavan, “Phasecam3d—Learning phase masks for passive single view depth estimation,” in _IEEE Intl. Conf. Computational Photography (ICCP)_ , 2019\.
  * Chang and Wetzstein [2019] J. Chang and G. Wetzstein, “Deep optics for monocular depth estimation and 3d object detection,” in _IEEE Intl. Conf. Computer Vision (ICCV)_ , 2019.
  * Huang et al. [2008] B. Huang, W. Wang, M. Bates, and X. Zhuang, “Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy,” _Science_ , vol. 319, no. 5864, pp. 810–813, 2008.
  * Betzig et al. [2006] E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” _Science_ , vol. 313, no. 5793, pp. 1642–1645, 2006.
  * Hess et al. [2006] S. T. Hess, T. P. Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation localization microscopy,” _Biophysical J._ , vol. 91, no. 11, pp. 4258–4272, 2006.
  * Juette et al. [2008] M. F. Juette, T. J. Gould, M. D. Lessard, M. J. Mlodzianoski, B. S. Nagpure, B. T. Bennett, S. T. Hess, and J. Bewersdorf, “Three-dimensional sub–100 nm resolution fluorescence microscopy of thick samples,” _Nature Methods_ , vol. 5, no. 6, pp. 527–529, 2008.
  * Girkin and Carvalho [2018] J. M. Girkin and M. T. Carvalho, “The light-sheet microscopy revolution,” _J. Optics_ , vol. 20, no. 5, p. 053002, 2018.
  * Di Giovanna et al. [2018] A. P. Di Giovanna, A. Tibo, L. Silvestri, M. C. Müllenbroich, I. Costantini, A. L. A. Mascaro, L. Sacconi, P. Frasconi, and F. S. Pavone, “Whole-brain vasculature reconstruction at the single capillary level,” _Scientific Reports_ , vol. 8, no. 1, pp. 1–11, 2018.
  * Lugo-Hernandez et al. [2017] E. Lugo-Hernandez, A. Squire, N. Hagemann, A. Brenzel, M. Sardari, J. Schlechter, E. H. Sanchez-Mendoza, M. Gunzer, A. Faissner, and D. M. Hermann, “3d visualization and quantification of microvessels in the whole ischemic mouse brain using solvent-based clearing and light sheet microscopy,” _J. Cerebral Blood Flow & Metabolism_, vol. 37, no. 10, pp. 3355–3367, 2017.
  * St. Croix et al. [2005] C. M. St. Croix, S. H. Shand, and S. C. Watkins, “Confocal microscopy: comparisons, applications, and problems,” _Biotechniques_ , vol. 39, no. 6, pp. S2–S5, 2005.
  * Kelch et al. [2015] I. D. Kelch, G. Bogle, G. B. Sands, A. R. Phillips, I. J. LeGrice, and P. R. Dunbar, “Organ-wide 3d-imaging and topological analysis of the continuous microvascular network in a murine lymph node,” _Scientific Reports_ , vol. 5, no. 1, pp. 1–19, 2015.
  * Bassi et al. [2011] A. Bassi, L. Fieramonti, C. D’Andrea, G. Valentini, and M. Mione, “In vivo label-free three-dimensional imaging of zebrafish vasculature with optical projection tomography,” _J. Biomedical Optics_ , vol. 16, no. 10, p. 100502, 2011.
  * Makita et al. [2006] S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” _Opt. Express_ , vol. 14, no. 17, pp. 7821–7840, 2006.
  * Grover et al. [2011] G. Grover, S. Quirin, C. Fiedler, and R. Piestun, “Photon efficient double-helix psf microscopy with application to 3d photo-activation localization imaging,” _Biomedical Optics Express_ , vol. 2, no. 11, pp. 3010–3020, 2011.
  * Gustavsson et al. [2018] A.-K. Gustavsson, P. N. Petrov, M. Y. Lee, Y. Shechtman, and W. Moerner, “3d single-molecule super-resolution microscopy with a tilted light sheet,” _Nature Communications_ , vol. 9, no. 1, pp. 1–8, 2018.
  * Bennett et al. [2020] H. W. Bennett, A.-K. Gustavsson, C. A. Bayas, P. N. Petrov, N. Mooney, W. Moerner, and P. K. Jackson, “Novel fibrillar structure in the inversin compartment of primary cilia revealed by 3d single-molecule superresolution microscopy,” _Molecular Biology of the Cell_ , vol. 31, no. 7, pp. 619–639, 2020.
  * Quirin and Piestun [2013] S. Quirin and R. Piestun, “Depth estimation and image recovery using broadband, incoherent illumination with engineered point spread functions,” _Appl. Optics_ , vol. 52, no. 1, pp. A367–A376, 2013.
  * Berlich et al. [2016] R. Berlich, A. Bräuer, and S. Stallinga, “Single shot three-dimensional imaging using an engineered point spread function,” _Optics Express_ , vol. 24, no. 6, pp. 5946–5960, 2016.
  * Wang et al. [2017] Z. Wang, Y. Cai, Y. Liang, X. Zhou, S. Yan, D. Dan, P. R. Bianco, M. Lei, and B. Yao, “Single shot, three-dimensional fluorescence microscopy with a spatially rotating point spread function,” _Biomedical Optics Express_ , vol. 8, no. 12, pp. 5493–5506, 2017.
  * Roider et al. [2014] C. Roider, A. Jesacher, S. Bernet, and M. Ritsch-Marte, “Axial super-localisation using rotating point spread functions shaped by polarisation-dependent phase modulation,” _Optics Express_ , vol. 22, no. 4, pp. 4029–4037, 2014.
  * Ikoma et al. [2021] H. Ikoma, T. Kudo, Y. Peng, M. Broxton, and G. Wetzstein, “Deep learning multi-shot 3d localization microscopy using hybrid optical–electronic computing,” _Optics Letters_ , vol. 46, no. 24, pp. 6023–6026, 2021.
  * Piestun et al. [2000] R. Piestun, Y. Y. Schechner, and J. Shamir, “Propagation-invariant wave fields with finite energy,” _JOSA A_ , vol. 17, no. 2, pp. 294–303, 2000.
  * Yanny et al. [2020] K. Yanny, N. Antipa, W. Liberti, S. Dehaeck, K. Monakhova, F. L. Liu, K. Shen, R. Ng, and L. Waller, “Miniscope3d: optimized single-shot miniature 3d fluorescence microscopy,” _Light: Science & Applications_, vol. 9, no. 1, pp. 1–13, 2020.
  * Xue et al. [2020] Y. Xue, I. G. Davison, D. A. Boas, and L. Tian, “Single-shot 3d wide-field fluorescence imaging with a computational miniature mesoscope,” _Science Advances_ , vol. 6, no. 43, p. eabb7508, 2020.
  * Zhang et al. [2020] K. Zhang, M. R. Kellman, E. Bostan, and L. Waller, “3d fluorescence deconvolution with deep priors (conference presentation),” in _Three-Dimensional and Multidimensional Microscopy: Image Acquisition and Processing XXVII_ , 2020.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * Hamarneh and Jassi [2010] G. Hamarneh and P. Jassi, “Vascusynth: Simulating vascular trees for generating volumetric image data with ground-truth segmentation and tree analysis,” _Computerized Medical Imaging and Graphics_ , vol. 34, no. 8, pp. 605–616, 2010.
  * Jassi and Hamarneh [2011] P. Jassi and G. Hamarneh, “Vascusynth: Vascular tree synthesis software,” _Insight J._ , 2011.
  * Wang et al. [2003] Z. Wang, E. P. Simoncelli, and A. C. Bovik, “Multiscale structural similarity for image quality assessment,” in _Asilomar Conf. Signals, Systems & Computers_, 2003.
  * Peruz [2022] S. Peruz, “Shape from focus,” https://www.mathworks.com/matlabcentral/fileexchange/55103-shape-from-focus, 2022\.
  * Goodman [2005] J. W. Goodman, “Introduction to fourier optics, roberts & co,” _Publishers, Englewood, Colorado_ , 2005.
  * Bando et al. [2008] Y. Bando, B.-Y. Chen, and T. Nishita, “Extracting depth and matte using a color-filtered aperture,” in _ACM SIGGRAPH Asia 2008 papers_ , 2008, pp. 1–9.
  * [54] Nanoscribe gmbh. [Online]. Available: https://www.nanoscribe.de/

##  7 Supplementary: Methods

###  7.1 Properties of GL-based rotating PSFs

[43] presented a theory on paraxial rotating beams using GL modes, along with
a detailed estimation of rotation range, rotation rates, and beam scaling
rates. These theoretical values could be extended to give some theoretical
description about GL-based rotating PSFs. However, to the best of our
knowledge, this has not been shown. Here we present a calculation for the
rotation range and rates for GL-based rotating PSFs.

Consider a monochromatic paraxial beam with width $w_{0}$ and wavelength
$\lambda$, which is generated from GL modes that lie on a single line in the
GL modal plane with a slope $V_{1}$, expressed as:

$n_{j}=V_{1}m_{j}+V_{2}\ \ \ \ j=1,2,...$ (6)

where $V_{1},V_{2}$ are integer constants, and $\\{m_{j}\\}$ are non-negative
numbers in arithmetic progression. [43] showed that such a paraxial beam
rotates at a rate

$\frac{d\phi}{dz}=\frac{V_{1}}{1+(z/z_{R})^{2}}$ (7)

where $z_{R}=\frac{\pi w_{0}^{2}}{\lambda}$ is the Rayleigh length. A point
with coordinates $(\rho_{0},\phi_{0})$ at $z=0$, follows the trajectory [43]:

$\displaystyle\rho$ $\displaystyle=\rho_{0}\sqrt{1+(z/z_{R})^{2}}$ (8)
$\displaystyle\phi$ $\displaystyle=\phi_{0}+V_{1}(arctan(z/z_{R}))$ (9)

which implies that the maximum rotation possible (in one direction) is
$V_{1}(\pi/2)$.

Thus, for the rotating beam case, over a distance of $z_{R}$, the beam rotates
$V_{1}(\pi/4)$. The way to find out rotation and scaling rates for a GL-based
rotating PSF is to find the equivalent Rayleigh length $z_{R}$ for it in a
$4f$ optical system. Note that the GL modes arise from modal solutions of the
Fresnel diffraction operator[43]. The Fresnel diffraction integral is given as
follows [52]:

$U(u,v)=\frac{e^{jkz}}{j\lambda
z}e^{j\frac{k}{2z}(u^{2}+v^{2})}\mathcal{F}\\{U(x,y)e^{j\frac{k}{2z}(x^{2}+y^{2})}\\}_{\frac{u}{\lambda
z},\frac{v}{\lambda z}}$ (10)

where $j^{2}=-1$ and $k$ is the wave number. Eqn 10 determines the wave
profile $U(u,v)$ at $z$ distance away from an input wave profile $U(x,y)$.
Thus, the phase term that determines the _defocus_ or _propagation_ by a
distance $z$ is

$\frac{k}{2z}(u^{2}+v^{2})$ (11)

Similarly, upon solving the propagation integrals for a $4f$ system with
defocus, the phase term that affects the defocus is given by

$\text{exp}\left(j\frac{k}{2f}(\frac{\Delta z}{f})(u^{2}+v^{2})\right)$ (12)

where $\Delta z$ is the distance from the in-focus plane. Comparing Eqns 11
and 12, we can obtain the equivalent Rayleigh length in the $4f$ system (say
$z^{\prime}_{R}$)

$\displaystyle\frac{1}{z_{R}}=\frac{z^{\prime}_{R}}{f^{2}}\implies
z^{\prime}_{R}=f^{2}/z_{R}=\frac{\lambda f^{2}}{\pi w_{0}^{2}}$ (13)

Thus, given a GL-based rotating PSF mask designed with beam width $w_{0}$, for
wavelength $\lambda$ and $V_{1}$ slope in the GL modal plane:

  * •

The total rotation amount in one direction of defocus is $V_{1}(\pi/2)$

  * •

The angle of the rotating PSF
$\phi(z)=\phi_{0}+V_{1}(\text{arctan}(\frac{z}{\lambda f^{2}/\pi w_{0}^{2}}))$

  * •

And over a depth of $\frac{\lambda f^{2}}{\pi w_{0}^{2}}$ the PSF will rotate
by $V_{1}(\pi/4)$ radians.

###  7.2 Analysis of CRLBϕ

Figure 13: CRLBϕ comparison for various phase masks. The figure shows $\sqrt{CRLB_{\phi}}$ log-scale plots for line images, as a function of line depth $(z)$ and line orientation $(\theta)$. Individual insets depict the mean and standard deviation of the $\sqrt{\textrm{CRLB}_{\phi}}$ values, which are all very low. Figure 14: Simulation methodology to evaluate PS2F. Given a 3D volume, we used a simulated stack of depth-dependent PSFs to obtain the two polarization images. We then computed the depth at each spatial point with the two polarization images as inputs, resulting in the point cloud shown on the right. TABLE I: Depth Estimation performance, averaged over 10 scenes for each scene complexity (noise level$=0.02$). $B$ refers to the number of bifurcations (B) in the 3D vascular structure. Showing mean absolute error (MAE) and RMS error (RMSE) values (all units in mm), and multi-scale SSIM [50] (MS-SSIM) values. For the first two metrics, lower is better (indicated by $\downarrow$), and for MS-SSIM, higher is better (indicated by $\uparrow$). | PS2F(proposed) | DHPSF(x2) [7] | TetrapodPSF(x2) [10]
---|---|---|---
| MAE $\downarrow$ | RMSE $\downarrow$ | MS-SSIM $\uparrow$ | MAE $\downarrow$ | RMSE $\downarrow$ | MS-SSIM $\uparrow$ | MAE $\downarrow$ | RMSE $\downarrow$ | MS-SSIM $\uparrow$
B=6 | 0.198 | 0.292 | 0.990 | 0.284 | 0.46 | 0.990 | 0.352 | 0.541 | 0.990
B=16 | 0.225 | 0.334 | 0.980 | 0.494 | 0.781 | 0.979 | 0.532 | 0.811 | 0.979
B=26 | 0.211 | 0.321 | 0.972 | 0.507 | 0.839 | 0.971 | 0.574 | 0.872 | 0.970
B=36 | 0.224 | 0.344 | 0.969 | 0.548 | 0.904 | 0.967 | 0.585 | 0.906 | 0.966
B=46 | 0.261 | 0.413 | 0.966 | 0.564 | 0.933 | 0.964 | 0.769 | 1.129 | 0.963

In the CRLB analysis performed in Sec 3.3, we consider an analysis for CRLBz
only. Intuitively, the estimation of line orientation from the image of a line
patch is an easier task as compared to estimating line depth. This is because
line orientation can be readily estimated from the global structure present in
the line image. To verify this, we perform CRLBϕ calculations with the same
parameters ($N=100,000$ photon, $\beta=5$ photons/pixel). The
$\sqrt{\textrm{CRLB}_{\phi}}$ plots are shown in Supplementary Fig. 13. The
$\sqrt{\textrm{CRLB}_{\phi}}$ values are very low, indicating that our
intuition about line orientation estimation being an easier problem is
correct.

##  8 Supplementary: Simulations

###  8.1 VascuSynth dataset simulations

The VascuSynth dataset simulations pipeline is best illustrated by Fig. 14.

Table I shows the comparative results for the proposed PS2F, DHPSF, and the
Tetrapod PSF over three metrics - Mean Absolute Depth Error, RMS Depth Error,
and multi-scale SSIM [50].

We also perform simulations under varying levels of noise, to judge the
robustness of the proposed PS2F. We simulate vasculature scene renderings with
Poisson noise, and also Gaussian read-out noise of levels $0.02$, $0.05$, and
$0.1$ \- which correspond to PSNRs of 34dB, 26dB, and 20dB. Fig. 15 shows that
with increasing noise, the proposed PS2F reconstruction only worsens slightly,
showing the robustness of the proposed PSF and method to noise.

Figure 15: Depth Estimation performance, for different noise levels (averaged
over 5 scenes for each complexity level, shaded regions show variance across
the 5 scenes). This shows the robustness of the proposed PSF and technique.

###  8.2 Ablation study: Comparison with polarized stereo pair

The proposed PS2F has been created using a novel polarizer-phase mask design
and the usage of a polarization-based camera sensor. We perform an ablation
study, comparing our proposed PS2F (polarizer-phase mask design) with a PSF
created using just polarizer halves in the Fourier (pupil) plane. This creates
images in the $s$, $p$-polarized channels that correspond to different half
apertures (taking inspiration from [53]). With such a polarized stereo-vision
effect, where each polarization channel sees a slightly different perspective,
we attempt reconstruction of a 3D USAF skew target in simulation. The
reconstructed intensity plots (depth color-coded) for the PS2F and the
Polarized-Stereo PSF are shown in Fig. 16. A polarizer-phase mask design gives
better spatial resolution, and a $\sim 2\times$ lower depth RMS error in
reconstruction as compared to the polarizer-only design. This can possibly be
attributed to the fact that the synthetic polarized-stereo PSF pair has a
narrow baseline length, and also that the PSFs are no longer compact, leading
to worse spatial resolution.

Figure 16: Comparing polarizer-phase mask design vs. polarizer-only design:
Simulated 3D USAF target reconstructions shown for the PS2F case and for the
polarized-stereo (polarizers-only) case. Showing depth color-coded
reconstructed intensity plots. A polarizer-phase mask design gives a $\sim
2\times$ lower depth RMS error in reconstruction as compared to the polarizer-
only design, and poorer spatial resolution. Scale bars indicate $0.1$ mm.

##  9 Supplementary: Experiments

###  9.1 Pupil plane encoding using polarizers and phase mask

Figure 17: Polarizer-mask setup in the pupil plane. (A) shows the 3D printed
holder, in which two polarizer films (B) aligned to orthogonal polarizations
(0 and 90 degree) will be placed, along with the fabricated DHPSF mask (C).
(D) shows the entire setup kept in front of an LCD screen with polarized
illumination. This demonstrates how such a setup allows for light from only
one half of the mask to pass through. (E) shows an image of the polarizer-
phase mask after keeping it in the pupil plane (imaged using a polarization
camera).

Figure 17 illustrates how the polarizer-phase mask encoding was achieved.

###  9.2 Phase mask fabrication

We used a two-photon photolithography 3D printer (Photonic Professional GT,
Nanoscribe GmbH [54]) to print the phase mask structures using IP-Dip
photoresist on a $700~{}\mu$m thick fused silica substrate. A $3$ mm diameter
phase mask corresponding to the DHPSF phase mask profile was fabricated, with
$2~{}\mu$m discretization in the $xy$-dimensions. For reliable fabrication,
the phase profile was also quantized to 5 different levels.

###  9.3 Fluorescent sample preparation

The fluorescent bead sample was prepared by adding 10$\mu$m fluorescent beads
solutions with varying concentrations in polydimethylsiloxane (PDMS) (Sylgard,
Dow Corning; 10:1 elastomer:cross-linker weight ratio). The fluorescent cotton
strand 3D sample was prepared by using a fluorescent highlighter on a tiny
piece of a cotton ball, which consisted of strands with thicknesses ranging
from $10$-$30~{}\mu$m. The fluorescent-highlighted sample was then placed in a
PDMS solution. All samples were cured at room temperature for a minimum of 24
hours.

###  9.4 Hyperparameter selection

In reconstruction, the optimization objective (Eq 5 in the main paper) has two
hyperparameters - $\lambda_{TV}$ and $\lambda_{L1}$ \- these control total
variation (TV) and L1 regularization. The three broad types of experimental
reconstructions shown in the submission (USAF target, fluorescent beads,
fluorescent cotton strands) had different scene settings and properties, such
as transmission vs. fluorescence illumination (different background signal
statistics), extended vs. point-like scenes (different sparsity) and thus the
hyperparameters were set accordingly to -

  * •

USAF (planar) target: $(\lambda_{L1},\lambda_{TV})=(0.02,0.005)$

  * •

Fluorescent beads: $(\lambda_{L1},\lambda_{TV})=(0.05,0.00)$

  * •

Fluorescent strands: $(\lambda_{L1},\lambda_{TV})=(0.002,0.002)$

For the same type of scene, high sensitivity to hyperparameter weights was not
observed in our reconstructions.