# Quantitative Jones matrix imaging using vectorial Fourier ptychography

Xiang Dai  Department of Biomedical Engineering, Duke University, Durham, NC
27708  Shiqi Xu  Department of Biomedical Engineering, Duke University,
Durham, NC 27708  Xi Yang  Department of Biomedical Engineering, Duke
University, Durham, NC 27708  Kevin C. Zhou  Department of Biomedical
Engineering, Duke University, Durham, NC 27708  Carolyn Glass  Department of
Pathology, Duke University, Durham, NC 27708  Pvan Chandra Konda  Department
of Biomedical Engineering, Duke University, Durham, NC 27708  Roarke
Horstmeyer  Department of Biomedical Engineering, Duke University, Durham, NC
27708  Corresponding author<EMAIL_ADDRESS>
###### Abstract

This paper presents a microscopic imaging technique that uses variable-angle
illumination to recover the complex polarimetric properties of a specimen at
high resolution and over a large field-of-view. The approach extends Fourier
ptychography, which is a synthetic aperture-based imaging approach to improve
resolution with phaseless measurements, to additionally account for the
vectorial nature of light. After images are acquired using a standard
microscope outfitted with an LED illumination array and two polarizers, our
vectorial Fourier Ptychography (vFP) algorithm solves for the complex 2x2
Jones matrix of the anisotropic specimen of interest at each resolved spatial
location. We introduce a new sequential Gauss-Newton-based solver that
additionally jointly estimates and removes polarization-dependent imaging
system aberrations. We demonstrate effective vFP performance by generating
large-area (29 mm2), high-resolution (1.24 $\mu$m full-pitch) reconstructions
of sample absorption, phase, orientation, diattenuation, and retardance for a
variety of calibration samples and biological specimens.

## Introduction

Polarization imaging methods provide a useful means to access information
about the molecular arrangement of anisotropic samples in a label-free manner
[1], with applications spanning pathology [2, 3, 4, 5, 6, 7, 8, 9],
developmental biology [10], and mineralogy [11]. As such, polarization is an
important and exciting intrinsic contrast mechanism in microscopy, especially
in combination with phase contrast and fluorescence imaging approaches [12].
However, polarimetric imaging over wide areas at high resolution remains an
outstanding challenge. Current microscopes that image at micrometer-level
detail can observe a field-of-view (FOV) that covers several square
millimeters [13], which makes polarimetric imaging of large specimens
difficult. To overcome this challenge, most existing approaches scan in a
step-and-repeat fashion to create large-area image composites [14]. As a fully
quantitative polarimetric recording typically requires multiple snapshots per
measurement, this scanning process can be tedious and time-consuming. What’s
more, it is also challenging to obtain phase-sensitive polarimetric
measurements (e.g., retardance) with relatively simple microscope hardware, to
step-and-scan while maintaining accurate phase sensitivity, and to minimize
the impact of polarization-dependent properties of the imaging system on the
acquired measurements.

Figure 1: Overview of vectorial Fourier ptychography (vFP). (a) vFP optical
setup consists of a standard polarized light microscope modified with an LED
array source to sequentially illuminate specimens with light polarized by
generator. Light scattered in polarization-dependent manner passes through an
analyzer before being recorded by the image sensor. (b) Images are acquired
for each LED under four different generator-analyzer configurations. (c)
4-channel Jones Matrix vFP reconstruction of the cross-section of a broad bean
root (color shows phase). (d) Various polarization-dependent sample properties
extracted from the Jones matrix using eigenvector analysis. (e) vFP
simultaneously measures a polarization-dependent pupil function as it removes
associated per-channel aberrations.

In this work, we present a new approach for phase-sensitive polarimetric
imaging that can address the above challenges and which produces high-
resolution, wide area polarimetric measurements using relatively simple
hardware. Our approach effectively extends an existing computational imaging
method, termed Fourier ptychography [15], to account for the vectorial nature
of light during high resolution reconstruction. This extension, termed
vectorial Fourier ptychography (vFP), first acquires a sequence of intensity
images under variable-angle illumination provided by an LED array. It then
utilizes a phase retrieval algorithm to recover complex polarimetric image
reconstructions, summarized by a complex-valued Jones Matrix at each resolved
pixel, from which one can directly measure quantities such as specimen
retardance, orientation and diattenuation. The resolution of these
reconstructions surpasses the standard resolution limit defined by the
microscope objective lens: here, we use a 4X 0.1 NA objective lens and
$f=100mm$ tube lens to provide 2X overall system magnification (29 $mm^{2}$
FOV), with an expected two-point resolution limit of 6.32$\mu$m. We
reconstruct complex-valued images with a final two-point resolution of
1.24$\mu$m, typically associated with 10X-20X objective lens performance. Our
computational strategy also jointly estimates and removes the imaging system’s
polarization-dependent aberrations, which points the way towards a new class
of simple and relatively inexpensive microscopes that can produce high-quality
and quantitative complex polarimetric measurements over large areas without
step-and-repeat scanning.

In the remainder of this paper, we will first review existing approaches for
polarimetric imaging, then present the forward and inverse mathematical models
for vectorial Fourier ptychography, before experimentally demonstrating its
effective quantitative performance.

##  1 Background

There are three general categories of polarization-sensitive microscopic
imaging. The first category captures qualitative polarimetric images,
including the use of differential interference contrast (DIC), and placing a
specimen between two polarizing films to visually assess variations of
transmitted light [16]. To obtain quantitative polarimetric measurements, a
second category of approach relies upon incoherent illumination to measure a
specimen’s Muller matrix [17], which describes the specimen’s optical response
via a 4x4 real-valued matrix at each spatial location of interest. In this
class of incoherent approaches, it is typical to describe light via a 4-entry
Stokes vector, and a complete specimen reconstruction typically requires at
least sixteen measurements per imaged pixel [18]. There have been a variety of
optical geometries developed for full Mueller matrix-based polarization
measurement [19, 20, 1, 18, 9]. By adopting circularly polarized specimen
illumination and applying certain specimen approximations, various setups can
obtain fewer image measurements to obtain specimen properties like
birefringence and orientation[9].

A third category of quantitative polarimetric microscopy relies upon coherent
light and can be described using the Jones Matrix formalism [21]. Within this
category, the optical field is commonly modeled as two-dimensional complex
vector, and associated optical transforms (including the interaction of light
with thin specimens) as $2\times 2$ complex matrix transforms. Typically,
phase-sensitive measurements are needed to infer the $2\times 2$ specimen
Jones Matrix at each resolved location in the specimen. Prior work has
considered this type of approach while detecting light with interferometry
[22, 23, 24, 25, 26, 27, 28, 29], optical coherence tomography [30, 20], and
laser scanning confocal microscopy [31], for example. Unlike incoherent
techniques, phase-sensitive methods enable direct measurement of specimen
phase retardance for topological information, as well as specimen
diattenuation, which is the complex differential transmittance within
independent polarization channels.

To alleviate the complexity introduced by interferometric detection methods, a
number of strategies have been recently developed that jointly utilize a post-
processing algorithm to assist with recovering a specimen’s complex
polarization-dependent properties [32, 33, 34, 35, 36, 37, 38]. These
computational imaging-based techniques offer a variety of unique benefits,
such as improved imaging speed and reduced system complexity. In this work, we
achieve complex-valued Jones matrix imaging at high resolution over a large
field-of-view by using principles from Fourier ptychography (FP) [15]. Recent
work has considered how FP can improve polarimetric imaging [36], but did not
demonstrate the ability to recover the specimen’s complex-valued Jones matrix,
like many of the interferometric approaches above. Complex Jones matrix
reconstruction offers a complete picture of the polarimetric properties of
thin specimens and has been demonstrated in scanning-based ptychography [39,
37], which uses a completely different measurement arrangement but is
mathematically similar to FP [40, 41]. In this work, we build upon these prior
demonstrations in ptychography to jointly reconstruct quantitative, complex-
valued Jones matrix measurements at high-resolution over a large area, along
with the additional benefit of computationally accounting for polarization-
dependent imaging system aberrations.

##  2 Methods

###  2.1 vFP forward model

We begin by developing a mathematical expression for Fourier ptychographic
imaging using Jones calculus, which will allow us to describe our data
collection strategy using vectorial optical fields. A diagram of the vFP
optical setup is in Fig. 1(a). Jones calculus represents an optical field as a
$2\times 1$ complex vector containing two orthogonal components, which we
define here along the $x$ and $y$ axes. The light’s polarization state is
defined by the amplitude and phase of these two components at each spatial
location of interest. General optical transformations (e.g., passage through
an optical component) are summarized by $2\times 2$ complex matrices. We note
that the Jones calculus employed here does not directly consider depolarizing
effects, or 3D effects, for example, which can easily be added in at the
expense of a more complex model.

Our setup consists of a programmable LED array as illumination source, a
generator polarizer placed between the sample and the LED array, standard
microscope optics, a CMOS camera, and an analyzer polarizer placed between the
optics and camera. The $n$th tilted plane wave emitted by the $n$th LED within
the array is described by $\text{exp}(i\mathbf{k}_{n}\cdot\mathbf{r})$, which
is polarized by the generator to form a $2\times 1$ vector field,
$\mathbf{G_{n}^{m}}=\left[g_{x}^{m}\text{exp}(i\mathbf{k}_{n}\cdot\mathbf{r}),g_{y}^{m}\text{exp}(i\mathbf{k}_{n}\cdot\mathbf{r})\right]^{T}$,
which is obliquely-incident on the sample of interest. Here,
$m\in\\{0^{\circ},45^{\circ},...\\}$ specifies a particular polarization
state, and we use $\mathbf{g_{m}}=\left[g_{x}^{m},g_{y}^{m}\right]$ as
complex-valued scalars to define the appropriate generator-dependent weighting
for each vector component. The polarized plane wave $\mathbf{G_{n}^{m}}$ then
interacts with the thin specimen of interest, which we model as a product with
a $2\times 2$ complex matrix $\boldsymbol{\bar{O}}(\mathbf{r})$.

The resulting field then propagates to the imaging system’s back focal plane,
which we model (under certain approximations) with a $2\times 2$ matrix
$\mathbf{F}$ that contains the two dimensional Fourier transform operator
$\mathcal{F}$ along its diagonal. Scalar FP models imaging system aberrations
with a scalar pupil function $P(k)$ [42]. In this work, we extend this model
to also allow for polarization-dependent aberrations with a $2\times 2$
complex pupil function matrix $\mathbf{P}(\mathbf{k})$, which accounts for
both standard aberrations (e.g., defocus, astigmatism) and more complex
birefringent effects, for example as introduced by plastic lenses.

After passage through the pupil function (i.e., multiplication with
$\mathbf{P}(\mathbf{k})$), the vectorized field then propagates to the image
plane, which we model here with a $2\times 2$ (inverse) Fourier transform
matrix, although alternative transforms can be directly inserted. The vector
field finally passes through a second “analyzer" polarizer, $\mathbf{A_{l}}$,
$(l\in\\{0^{\circ},45^{\circ},\mathrm{right\ circular,}...\\})$, to produce
the $n$th vector field $\mathbf{E_{n}^{(l,m)}}$ at the image plane,

$\mathbf{E_{n}^{(l,m)}}=\mathbf{A}_{l}\mathbf{F}^{-1}\mathbf{P}(\mathbf{k})\mathbf{F}\mathbf{\bar{O}}(\mathbf{r})\mathbf{G_{n}^{m}}$
(1)

where

$\bar{\mathbf{O}}(\mathbf{r})=\begin{bmatrix}&\mkern-18.0mu\bar{O}_{xx}(\mathbf{r})&\bar{O}_{xy}(\mathbf{r})\\\
&\mkern-18.0mu\bar{O}_{yx}(\mathbf{r})&\bar{O}_{yy}(\mathbf{r})\end{bmatrix},\mathbf{P}(\mathbf{k})=\begin{bmatrix}&\mkern-18.0muP_{xx}(\mathbf{k})&P_{xy}(\mathbf{k})\\\
&\mkern-18.0muP_{yx}(\mathbf{k})&P_{yy}(\mathbf{k})\end{bmatrix}$ (2)

are the Jones matrices of the spectrum of the sample and pupil. We can take
advantage of the Fourier shift theorem to remove the first Fourier transform
matrix, and instead describe a shifted sample spectrum across the pupil plane,
caused by each plane wave traveling at angle $\mathbf{k}_{n}$, as
$\mathbf{O}(\mathbf{k}-\mathbf{k}_{n})$, where $\mathbf{O}$ is the Fourier
transform of $\mathbf{\bar{O}}$. This leads to the following compact form of
the vFP forward model for each detected intensity image $i$ under illumination
from the $n$th LED angle with generator and analyzer configurations $(l,m$):

$i_{n}^{l,m}\left(\mathbf{O},\mathbf{P}\right)=\lvert{\mathbf{a}_{l}}^{T}\mathbf{F}^{-1}\mathbf{P}(\mathbf{k})\mathbf{O}(\mathbf{k}-\mathbf{k}_{n})\mathbf{g}_{m}\rvert^{2},$
(3)

where $\mathbf{a}_{l}$ and $\mathbf{g}_{m}$ is the Jones vector of the
$l^{th}$ analyzer and $g^{th}$ generator, respectively. In our first
demonstration we use linear polarizers for both generator and analyzer, as
will be discussed in the next subsection, we have
$\mathbf{a}_{l}=\left[\cos{\alpha_{l}},\sin{\alpha_{l}}\right]^{T}$, and
$\mathbf{g}_{m}=\left[\cos{\gamma_{m}},\sin{\gamma_{m}}\right]^{T}$, where
$\alpha_{l}$ and $\gamma_{m}$ are the angles of the linear polarizers of the
$l^{th}$ analyzer and $g^{th}$ generator, respectively. The goal of vFP is to
jointly recover the complex $2\times 2$ matrices describing the specimen and
pupil function, $\boldsymbol{O}$ and $\mathbf{P}$, from an acquired image
dataset $i_{n}^{l,m}$ for multiple illumination angles $n=1$ to $N$ and
multiple configurations of the generator $l$ and the analyzer $m$. There are
various possible selections for these measurement variables, and we next
describe one effective measurement strategy that yields a well-conditioned
image dataset for inverse problem estimation.

###  2.2 Measurement strategy

Obtaining four unique measurements per illumination angle is one effective
strategy to solve for the four unknowns summarized by the sample matrix
$\mathbf{O}(\mathbf{k})$. There are, of course, many possible configurations
of the generator and analyzer to utilize for such measurements[39]. Specific
selections significantly reduce the complexity of Eq. (3). For example,
selecting the generator($G$) and the analyzer($A$) to produce and admit zero-
degree linearly polarized light will result in $\mathbf{g}_{m}=[1,0]^{T}$ and
$\mathbf{a}_{l}=[1,0]^{T}$. The resulting image intensity will be described by
$i_{n}^{l=0\degree,m=0\degree}=\lvert\mathcal{F}^{-1}[P_{xx}O_{xx}(\mathbf{k}-\mathbf{k}_{n})+P_{xy}O_{yx}(\mathbf{k}-\mathbf{k}_{n})]\rvert^{2}$
for plane wave illumination at an angle $\mathbf{k}_{n}$. Setting the
generator and the analyzer at 90 degrees will likewise yield additional
complementary information about $\mathbf{O}$ and $\mathbf{P}$. The above two
configurations provide the helpful benefit of high optical transmission for
most specimens of interest, which facilitates the use of lower exposure times
and reduced noise. To balance the experimental benefits of high optical
transmission with the computational benefits of measurement redundancy, we
next set the generator ($G$) at $45\degree$ for our remaining two
measurements. Our four polarization configurations were thus,
$(l,m)\in\\{G\mkern-6.0mu:0\degree,A\mkern-6.0mu:0\degree;G\mkern-6.0mu:90\degree,A\mkern-6.0mu:90\degree;G\mkern-6.0mu:45\degree,A\mkern-6.0mu:0\degree;G\mkern-6.0mu:45\degree,A\mkern-6.0mu:90\degree\\}$,
as shown in Fig. 1. We note that alternative generator/analyzer configurations
are certainly possible and compatible with our forward and inverse solver, and
will be the subject of future study.

###  2.3 Inverse problem

Using the measurement strategy proposed in Section 22.2, the four amplitude
measurements for the $n^{th}$ LED are

$\displaystyle i_{n}^{0\degree,0\degree}(\mathbf{r})=\lvert\mathcal{F}^{-1}$
$\displaystyle\\{P_{xx}(\mathbf{k})\cdot
O_{xx}(\mathbf{k}-\mathbf{k}_{n})+P_{xy}(\mathbf{k})\cdot
O_{yx}(\mathbf{k}-\mathbf{k}_{n})\\}\rvert^{2}$ (4a) $\displaystyle
i_{n}^{90\degree,90\degree}(\mathbf{r})=\lvert\mathcal{F}^{-1}$
$\displaystyle\\{P_{yx}(\mathbf{k})\cdot
O_{xy}(\mathbf{k}-\mathbf{k}_{n})+P_{yy}(\mathbf{k})\cdot
O_{yy}(\mathbf{k}-\mathbf{k}_{n})\\}\rvert^{2}$ (4b) $\displaystyle
i_{n}^{45\degree,90\degree}(\mathbf{r})=\lvert\mathcal{F}^{-1}$
$\displaystyle\\{P_{xx}(\mathbf{k})\cdot
O_{xx}(\mathbf{k}-\mathbf{k}_{n})+P_{xx}(\mathbf{k})\cdot
O_{xy}(\mathbf{k}-\mathbf{k}_{n})$ $\displaystyle+$ $\displaystyle
P_{xy}(\mathbf{k})\cdot
O_{yx}(\mathbf{k}-\mathbf{k}_{n})+P_{xy}(\mathbf{k})\cdot
O_{yy}(\mathbf{k}-\mathbf{k}_{n})\\}\rvert^{2}$ (4c) $\displaystyle
i_{n}^{45\degree,90\degree}(\mathbf{r})=\lvert\mathcal{F}^{-1}$
$\displaystyle\\{P_{yx}(\mathbf{k})\cdot
O_{xx}(\mathbf{k}-\mathbf{k}_{n})+P_{yx}(\mathbf{k})\cdot
O_{xy}(\mathbf{k}-\mathbf{k}_{n})$ $\displaystyle+$ $\displaystyle
P_{yy}(\mathbf{k})\cdot
O_{yx}(\mathbf{k}-\mathbf{k}_{n})+P_{yy}(\mathbf{k})\cdot
O_{yy}(\mathbf{k}-\mathbf{k}_{n})\\}\rvert^{2}.$ (4d)

We have developed an iterative algorithm to estimate the complex-valued object
and pupil matrices (extending scalar phase recovery [43, 44] to the vector
case) by minimizing the Euclidean distance between predicted and measured
amplitude,

$\mathcal{D}(\mathbf{O},\mathbf{P})=\frac{1}{2N}\sum_{n}^{N}\sum_{(l,m)}\lVert\sqrt{i_{n}^{l,m}\left(\mathbf{O},\mathbf{P}\right)}-\sqrt{\hat{i}(\mathbf{r})_{n}^{l,m}}\rVert^{2}_{2}.$
(5)

Here, $N$ is the total number of illumination angles, while
${i}_{n}^{l,m}\left(\mathbf{O},\mathbf{P}\right)$ and $\hat{i}_{n}^{l,m}$ are
the vFP forward model prediction (eq.4) and the actual experimental image
measurements for the $n$th LED illumination using the $(l,m)$th generator-
analyzer configuration. To minimize the loss function
$\mathcal{D}(\mathbf{O},\mathbf{P})$, we use a Gauss-Newton-based sequential
solver, which alternatively updates the object and pupil matrix estimations
based upon the 4 polarization measurements from each LED. Our method is
inspired by prior ptychographic phase retrieval solvers like the sequential
Gauss-Newton(GN) method [43] and regularized ptychographic iterative engine
(rPIE) [44], which use CR-calculus [45] to estimate the gradient of the loss
with respect to both sample and pupil. In our approach, we stack 4 image
measurements into a per-LED measurement vector and also compute a complex-
valued loss function gradient. Unlike prior work, minimizing the loss function
with a GN method requires inversion of a non-diagonal system matrix due mixing
of polarization states. This adds several complexities as compared to matching
computations with scalar fields, which can effectively use a direct division
for inversion [43, 44]. We refer the interested reader to a complete
presentation and derivation of our algorithm in the Supplementary Material.

Figure 2: Verification of vFP diattenuation estimation accuracy, as well as
resolution gain, using a resolution target and linear polarizer. (a)
Diattenuation (brightness) and orientation (hue) of vFP reconstruction with
USAF target half covered (in the upper left) by a linear polarizer. (b)
Incoherent sum of squared channel amplitudes (brightness) and orientation
(hue) of USAF target fully covered by the linear polarizer, rotated by 48
degrees as compared to placement in (a). (c) Zoom-in of the raw image captured
under illumination from center LED. (d) Zoom-in of the vFP reconstruction in
(b) verifies resolution gain.

Once we reconstruct the $2\times 2$ Jones matrix for the specimen and pupil,
we use a simple eigendecomposition method [1] to derive the homogeneous
anisotropic properties of the specimen, including retardance, diattenuation
and orientation. Specifically, for the 2x2 Jones matrix
$\bar{\mathbf{O}}(\mathbf{r})$ per reconstructed pixel, we first compute
eigenvalues $\xi_{q}$, $\xi_{r}$ and eigenvectors $\mathbf{E}_{q}$,
$\mathbf{E}_{r}$. Under the assumption of homogeneous anisotropic specimen,
the calculated eigenvectors indicate two perpendicular specimen axes with
associated eigenvalues that can be used to compute a per-pixel retardance,

$B=\frac{\lambda}{2\pi}\lvert\angle\xi_{q}-\angle\xi_{r}\rvert,$ (6)

where $\angle$ denotes complex angle. Similarly, the two eigenvalues may be
used to compute an amplitude-dependent diattenuation measurement:

$D=\Bigg{\lvert}\frac{\lvert\xi_{q}\rvert^{2}-\lvert\xi_{r}\rvert^{2}}{\lvert\xi_{q}\rvert^{2}+\lvert\xi_{r}\rvert^{2}}\Bigg{\rvert}.$
(7)

In this work, we adopt a definition for orientation as the angle between the
reference axis given by the linear polarizer within the setup and the slow
axis within the specimen plane, where the slow axis is defined by the
direction of the eigenvector whose corresponding eigenvalue phase is the
larger of $\xi_{q}$ and $\xi_{r}$. Mathematically, orientation to the slow
axis $\omega$ can be computed as the angle of the eigenvector,

$\omega=\tan^{-1}(\nicefrac{{{E}_{i,x}}}{{{E}_{i,y}}}),\quad
i=\operatorname*{arg\,max}_{i\in\\{q,r\\}}\angle\xi_{i},$ (8)

where $E_{\cdot,x}$ and $E_{\cdot,y}$ are the real $x$ and $y$ components of
the eigenvector, respectively.

Figure 3: Orientation measurement validation. (a) vFP orientation measurement
of monosodium urate (MSU) crystals($y$-axis) exhibit expected linear
relationship with crystal rotation angle ($x$-axis). Each point reports
average and standard deviation (error bar) of per-crystal vFP orientation
measurement averaged over pixels within each crystal (examples in 3 insets,
color is orientation and intensity is retardance). (b) Raw image and vFP
reconstruction of human connective tissue specimen. (b1) Example raw image and
(b2) zoom-in. (b3) vFP reconstruction shown as incoherent summation of matrix
elements. (c) vFP retardance and orientation of the same area as in (b2-3),
where (c1) shows retardance in grayscale intensity and orientation false-
colored and (c2, c3) highlight how specimen orientation follows fiber growth
direction and exhibits spatial correlations with retardance.

##  3 Results

We verified our vFP model and inverse solver in a series of experiments using
a standard microscope outfitted with an LED array. The vFP setup uses an
objective lens (4$\times$, $NA_{o}$=0.1) and digital camera containing 2.4
$\mu m$ pixels (Basler acA5472-17um, 5496 × 3672 pixel count). The LED array
contains 15 × 15 surface-mounted elements (center wavelengths: 632 nm, 525 nm
and 470 nm per LED color with 4 mm LED pitch). We positioned the LED array 68
mm beneath the sample to create a maximum illumination NA of $NA_{i}=$ 0.41.
This leads to an effective synthetic NA for our Fourier ptychographic
reconstructions of $NA=NA_{i}+NA_{o}=0.51$, or a 5x improvement in spatial
resolution in each lateral dimension. Linear wire grid polarizers
(polarization wavelength range between 420 - 700 nm) were inserted directly
before the LED array (generator) and before the camera (analyzer). Both linear
polarizers were installed in rotation mounts with resonant piezoelectric
motors (bidirectional accuracy, 0.002∘).

###  3.1 Validation of vFP resolution and diattenuation

To first validate that the polarization-dependent metrics of vFP exhibit the
expected resolution enhancement commonly observed with scalar FP, we imaged a
sample consisting of a plastic linear polarizer (Edmund Optics) placed above a
high-resolution absorptive USAF target (Ready Optics). The linear polarizer
provides a simple means to create an expected constant orientation and a
uniform diattenuation across the regions it covers, while the USAF target
provides a standard resolution measure.

The expected resolution gain of a vFP reconstruction follows the same
principles as with scalar FP [46], wherein the smallest resolvable full-pitch
resolution (i.e., distance between the center of two similar bars) can be
approximated by, $\frac{\lambda}{NA_{i}+NA_{o}}$. Using red LED illumination
(center wavelength 632 nm) in this experiment, we accordingly expect a 1.24
$\mu$m full-pitch resolution cutoff after executing vFP, which corresponds to
a cutoff slightly above USAF target Group 9 Element 5. The vFP orientation
reconstruction achieves this resolution as shown in Fig. 2(d), which provides
an appreciable gain over the raw captured image intensity in (c).

This linear polarizer experiment additionally provides a means to assess
quantitative diattenuation and orientation measurement accuracy. As only the
upper-left of the USAF target is covered by the polarizer in Fig. 2(a), we
expect a diattenuation value closer to 1 within this region, and close to 0 in
the lower right, which is uncovered. The mean diattenuation and variance,
averaged across a 5002 pixel area within these two respective regions, are
$(\bar{D}=0.9822,\sigma^{2}_{D}=0.001)$ and
$(\bar{D}=0.0693,\sigma^{2}_{D}=0.011)$, respectively. We hypothesize that
this experimentally measured diattenuation of less than 1 is due to the use of
a plastic linear polarizer, whose transmission at 632 nm is close to 45%,
rather than 50%. In addition, the transition boundary between the image area
covered and uncovered by linear polarizer exhibits inaccurate values due to
angle-dependent shadowing effects introduced by its finite thickness (i.e.,
the polarizer does not obey the required thin specimen approximation, leading
to such artifacts). The diattenuation estimates at locations of absorptive
USAF bars may also be inaccurate, as the bars are optically opaque.

Figure 4: vFP retardance and phase validation experiments. (a) Retardance (top
two rows) and phase (bottom two rows) reconstructions of 10-µm monodisperse
polystyrene microspheres immersed in oils with refractive indices (RIs) of
1.56, 1.57 and 1.58. Profiles through select beads show RI-dependent
retardance shift at bead-oil boundary measured with red-dotted line, while
phase profiles closely match ground truth shift. (b) Reconstructed maximum
retardance shift (red dashed line in (a), average of 20 beads per point)
follows expected linear trend with RI difference between microspheres and oil.
(c) Reconstructed max. phase shift also follows ground truth (average of 20
beads per point).

###  3.2 Quantitative orientation measurement

Orientation is a useful polarimetric quantity used within material and
biological specimen analysis. To further investigate the quantitative accuracy
of vFP for specimen orientation measurement, we next imaged a collection of
randomly-oriented monosodium urate (MSU) crystals, which are are rod-like
microstructures whose optical orientation properties are expected to exhibit a
linear relationship with crystal rotation angle within the x-y plane. A
summary of this experiment’s results is presented in Fig. 3(a). Images of 3
example reconstructed MSU crystals are included as insets (color denotes
reconstructed orientation and intensity denotes reconstructed retardance).
Plotting the vFP-measured per-crystal orientation versus the expected
orientation, as measured from crystal rotation angle within the x-y plane,
reveals a clear linear relationship. Measured orientation points were averaged
over all pixels within each crystal, where the error bar shows the resulting
standard deviation. Ground-truth orientations were computed by fitting a line
to each MSU crystal captured within a single FOV.

In a final experiment to assess vFP’s ability to measure orientation, we
imaged a thin section of human connective tissue containing both muscle and
collagen (Fig. 3(b-c)). In this type of tissue specimen, it is common to
observe polarization-dependent effects due in part to the response of collagen
fiber[47, 48] To improve visualization, we overlay orientation as an
undirected vector field on top of the retardance shown in grayscale in Fig.
3(c) to reveal a spatial correlation between retardance and orientation.
Reconstructed orientation vector alignment is consistent with fiber growth
direction within the connective tissue. We note the FOV in (c1) matches that
of the example raw image is shown in (b2) and the zoom of the vFP
reconstruction in (b3).

Figure 5: Large-FOV vFP reconstruction of thin cardiac tissue section stained
with Congo red dye. (a) Phase reconstruction. (b) Retardance reconstruction.
(c,d) Zoom-ins of (a). (e,f) Zoom-ins of (b). (g) Difference of the difference
between the real $xx$ and $yy$ components of the Jones matrix reconstruction
between red and green channels, corresponding to the red boxes in (a) and (b).
(h) Zoom-in of (g). (i) Validation via bright-field microscope imaging with
crossed polarizers. (j) Spatially-varying pupils, plotted as the magnitude and
phase of the Jones matrix components. White dots denote the center of the
pupil to improve visualization of shifts.

###  3.3 Retardance and phase validation

To validate vFP’s ability to measure specimen retardance, we followed the
findings of a recent study [12] that demonstrated a linear relationship
between retardance and refractive index (RI) difference at the boundary of
microspheres and their embedding medium. Though microspheres do not
intrinsically exhibit birefringence, a sharp RI transition at their boundary
will lead to a measurable retardance shift that is linearly proportional to
the RI difference at the transition boundary [49]. To explore this connection,
we prepared 3 specimens of oil-immersed polystyrene microspheres with oils of
3 different RI to generate 3 unique RI transition differences. Images of vFP-
reconstructed retardance are in Fig. 4(a) top, from which we observe a clear
peak of retardance that consistently exists at the bead profile edge (i.e.,
forms a ring around each bead). Profiles of randomly selected beads below
confirm such peak formation, from which we measure the maximum per-bead
retardance shift (red dashed line). Plotting the average maximum retardance
shift (20 beads per point) as a function of RI difference in Fig. 4(b)
confirms the previously observed linear trend. Through this experimental
analysis, we were additionally able to confirm the quantitative accuracy of
vFP reconstructed phase values, as shown at the bottom of Fig.4(a) and
summarized in plot of average maximum phase shift versus RI difference in Fig.
4(c).

###  3.4 Large-area amyloid plaque detection

We next applied vFP to reconstruct the vectorial response of a specimen that
is commonly viewed using polarization optics within the pathology lab: thin
sections of fixed cardiac tissue stained with Congo red dye. Congo Red is
known to stain specimen areas that contain amyloid fibrils in a manner that
highlights birefringence as anomalous colors (typically, as an “apple green”
shading) when viewed between two crossed polarizers using white light
illumination in a transmission geometry [50]. Visual identification of such
apple-green areas within the interstium and vessel walls of endomyocardial
tissue is often an important diagnostic indicator of the presence of cardiac
amyloidosis. In Fig. 5(a,b), we present a large-FOV vFP reconstruction (8.27
$mm^{2}$ shown) of the phase and retardance of Congo red-stained cardiac
tissue, wherein the arterial wall is clearly visible. Interestingly, in Fig.
5(f), the red blood cells (small, dark circular features) located within the
vessels each exhibit a higher retardance along their border, which is
consistent with our observations with polystyrene beads (Fig. 4). Examining
the real component of the specimen Jones matrix entries $\bar{O}_{xx}$ and
$\bar{O}_{yy}$ within this area shows minimal visual disparity under both
green and red illumination wavelengths. However, by computing the difference
($\bar{O}_{xx}-\bar{O}_{yy}$) for each color channel, and subtracting the red
channel’s difference from the green, it is possible to produce a high-contrast
map of Congo-red stained areas, as presented as a green overlay in Fig. 5(f).
Comparing the resulting channel difference map to what is typically observed
with polarization optics in a standard microscope (white-light illumination,
NA=0.25) demonstrates a close qualitative mapping. We hypothesize that vFP can
thus provide a useful tool for examining the polarimetric response of such
Congo-red-stained cardiac tissue over large areas at high resolution.

Reconstructing a large FOV with vFP additionally allows us to examine how the
recovered Jones matrix of the imaging system pupil function differs between
on-axis and off-axis sample positions (Fig. 5(j)). Here, we divided the
reconstruction FOV into $3\times 3$ regions and averaged all pupils within
each region. The Jones matrix pupil phase patterns experience a shift from the
center for off-axis locations [51]. Interestingly, the $J_{xy}$ and $J_{yx}$
pupil amplitudes are spatially varying, unlike $J_{xx}$ and $J_{yy}$. Higher-
NA objectives or cheaper plastic lenses may exhibit stronger polarization
effects, which will be the subject of future investigation.

##  4 Discussion and Conclusion

This paper proposes a new computational microscopy method, termed vectorial
Fourier ptychography (vFP), to measure the complex amplitude and phase
transmission of a thin sample along with its polarization properties contained
within its complex Jones matrix. As a synthetic aperture imaging technique,
vFP provides a means to obtain high-resolution Jones matrix measurements over
a wide FOV without any complex moving parts. Through a simple modification of
a conventional polarized light microscope (i.e., replacing the source with an
LED array) and a new computational reconstruction algorithm, we demonstrated
that our method can jointly recover a vectorial image of complex-valued 2x2
Jones matrices in a variety of calibrated phantoms and biological samples,
along with estimation of the 2x2 Jones matrix summarizing the imaging system’s
vectorial response (i.e., its polarization-dependent aberrations). Our method
thus exhibits capabilities that neither conventional FP nor polarized light
microscopy can achieve alone.

As this initial demonstration was a proof of concept, there are several
avenues for future direction. To improve acquisition speed, vFP can certainly
borrow LED multiplexing [43], multispectral, [52] or multi-aperture [53]
strategies previously demonstrated in conventional FP. Furthermore, the
approach could also utilize a polarization-sensitive digital image sensor [35,
38], which could reduce the amount of required captured data or otherwise
improve vFP reconstruction fidelity. In addition, a screen display-based
illumination mechanism could be adopted [54] to miniaturize the imaging
system. While the present work was restricted to thin samples, additional
theoretical development of the vFP framework may allow for 3D polarization
microscopy [25, 38]in the future. Moreover, additional development may also
account for the partial coherence of the LED illumination employed by vFP
within the Stokes formalism [38]. Finally, since vFP can jointly estimate and
account for the polarization properties of the imaging system, our method
could be applied to low-cost microscopes employing plastic lenses, which are
known to exhibit birefringence effects under strain. Given the experimental
simplicity of our method and the broad range of anisotropic biological
samples, vFP could be widely adopted as a significant extension of
conventional FP and polarized light microscopy.

##  5 Funding

This work was supported in part by a Duke-Coulter Translational Partnership
Grant and a 3M Non-tenured Faculty Award.

##  6 Acknowledgments

We would like to thank Dr. Shalin Mehta, Dr. Ian Sigal, Dr. Bin Yang, Dr. Li-
Hao Yeh, and Dr. Lars Loetgering for helpful guidance and feedback.

##  7 Disclosures

RH: Ramona Optics Inc. (F,I,P,S)

## References

  * [1] Russell A Chipman, Wai Sze Tiffany Lam, and Garam Young.  Polarized light and optical systems.  CRC press, 2018.
  * [2] Moshe Wolman.  Polarized light microscopy as a tool of diagnostic pathology.  Journal of Histochemistry & Cytochemistry, 23(1):21–50, 1975.
  * [3] Lee-Way Jin, Kacey A Claborn, Miki Kurimoto, Morten A Geday, Izumi Maezawa, Faranak Sohraby, Marcus Estrada, Werner Kaminksy, and Bart Kahr.  Imaging linear birefringence and dichroism in cerebral amyloid pathologies.  Proceedings of the National Academy of Sciences, 100(26):15294–15298, 2003.
  * [4] Bhari Sharanesha Manjunatha, Aditi Agrawal, Vandana Shah, et al.  Histopathological evaluation of collagen fibers using picrosirius red stain and polarizing microscopy in oral squamous cell carcinoma.  Journal of cancer research and therapeutics, 11(2):272, 2015.
  * [5] Chao He, Honghui He, Jintao Chang, Yang Dong, Shaoxiong Liu, Nan Zeng, Yonghong He, and Hui Ma.  Characterizing microstructures of cancerous tissues using multispectral transformed mueller matrix polarization parameters.  Biomedical optics express, 6(8):2934–2945, 2015.
  * [6] Yibo Zhang, Seung Yoon Celine Lee, Yun Zhang, Daniel Furst, John Fitzgerald, and Aydogan Ozcan.  Wide-field imaging of birefringent synovial fluid crystals using lens-free polarized microscopy for gout diagnosis.  Scientific reports, 6(1):1–14, 2016.
  * [7] Ning-Jiun Jan, Jonathan L Grimm, Huong Tran, Kira L Lathrop, Gadi Wollstein, Richard A Bilonick, Hiroshi Ishikawa, Larry Kagemann, Joel S Schuman, and Ian A Sigal.  Polarization microscopy for characterizing fiber orientation of ocular tissues.  Biomedical optics express, 6(12):4705–4718, 2015.
  * [8] Ali H Badreddine, Tomas Jordan, and Irving J Bigio.  Real-time imaging of action potentials in nerves using changes in birefringence.  Biomedical optics express, 7(5):1966–1973, 2016.
  * [9] Chao He, Honghui He, Jintao Chang, Binguo Chen, Hui Ma, and Martin J Booth.  Polarisation optics for biomedical and clinical applications: a review.  Light Sci Appl 10, 194, 2021.
  * [10] Maki Koike-Tani, Tomomi Tani, Shalin B Mehta, Amitabh Verma, and Rudolf Oldenbourg.  Polarized light microscopy in reproductive and developmental biology.  Molecular reproduction and development, 82(7-8):548–562, 2015.
  * [11] Neha Panwar and Ravi Sharma.  A review on influence of mineralogy and diagenesis on spectral induced polarization measurements in carbonate rocks.  In Petro-physics and Rock Physics of Carbonate Reservoirs, pages 115–125. Springer, 2020.
  * [12] Syuan-Ming Guo, Li-Hao Yeh, Jenny Folkesson, Ivan E Ivanov, Anitha P Krishnan, Matthew G Keefe, Ezzat Hashemi, David Shin, Bryant B Chhun, Nathan H Cho, et al.  Revealing architectural order with quantitative label-free imaging and deep learning.  elife, 9:e55502, 2020.
  * [13] Guoan Zheng, Xiaoze Ou, Roarke Horstmeyer, Jaebum Chung, and Changhuei Yang.  Fourier ptychographic microscopy: A gigapixel superscope for biomedicine.  Opt. Photon. News, 25(4):26–33, Apr 2014.
  * [14] Ning-Jiun Jan, Kira Lathrop, and Ian A Sigal.  Collagen architecture of the posterior pole: high-resolution wide field of view visualization and analysis using polarized light microscopy.  Investigative ophthalmology & visual science, 58(2):735–744, 2017\.
  * [15] Guoan Zheng, Roarke Horstmeyer, and Changhuei Yang.  Wide-field, high-resolution fourier ptychographic microscopy.  Nature photonics, 7(9):739–745, 2013.
  * [16] Michael Bass, Casimer DeCusatis, Jay Enoch, Vasudevan Lakshminarayanan, Guifang Li, Carolyn Macdonald, Virendra Mahajan, and Eric Van Stryland.  Handbook of optics, volume i: Geometrical and physical optics, polarized light, components and instruments (set).  McGraw-Hill, Inc., 2009.
  * [17] J Larry Pezzaniti and Russell A Chipman.  Mueller matrix imaging polarimetry.  Optical engineering, 34(6):1558–1568, 1995.
  * [18] Dennis H Goldstein.  Polarized light.  CRC press, 2017.
  * [19] Andreas H Hielscher, Angelia A Eick, Judith R Mourant, Dan Shen, James P Freyer, and Irving J Bigio.  Diffuse backscattering mueller matrices of highly scattering media.  Optics express, 1(13):441–453, 1997.
  * [20] Gang Yao and Lihong V Wang.  Two-dimensional depth-resolved mueller matrix characterization of biological tissue by optical coherence tomography.  Optics Letters, 24(8):537–539, 1999.
  * [21] R Clark Jones.  A new calculus for the treatment of optical systemsi. description and discussion of the calculus.  Josa, 31(7):488–493, 1941.
  * [22] Zhuo Wang, Larry J Millet, Martha U Gillette, and Gabriel Popescu.  Jones phase microscopy of transparent and anisotropic samples.  Optics letters, 33(11):1270–1272, 2008.
  * [23] Sherazade Aknoun, Michel Aurrand-Lions, Benoit Wattellier, and Serge Monneret.  Quantitative retardance imaging by means of quadri-wave lateral shearing interferometry for label-free fiber imaging in tissues.  Optics Communications, 422:17–27, 2018.
  * [24] Yuheng Jiao, Mikhail E Kandel, Xiaojun Liu, Wenlong Lu, and Gabriel Popescu.  Real-time jones phase microscopy for studying transparent and birefringent specimens.  Optics Express, 28(23):34190–34200, 2020.
  * [25] Amirhossein Saba, Joowon Lim, Ahmed B Ayoub, Elizabeth E Antoine, and Demetri Psaltis.  Polarization-sensitive optical diffraction tomography.  Optica, 8(3):402–408, 2021.
  * [26] Tristan Colomb, Pia Dahlgren, Didier Beghuin, Etienne Cuche, Pierre Marquet, and Christian Depeursinge.  Polarization imaging by use of digital holography.  Applied optics, 41(1):27–37, 2002.
  * [27] Youngchan Kim, Joonwoo Jeong, Jaeduck Jang, Mahn Won Kim, and YongKeun Park.  Polarization holographic microscopy for extracting spatio-temporally resolved jones matrix.  Optics Express, 20(9):9948–9955, 2012.
  * [28] Taeseok Daniel Yang, Kwanjun Park, Yong Guk Kang, Kyoung J Lee, Beop-Min Kim, and Youngwoon Choi.  Single-shot digital holographic microscopy for quantifying a spatially-resolved jones matrix of biological specimens.  Optics express, 24(25):29302–29311, 2016.
  * [29] Baoliang Ge, Qing Zhang, Rui Zhang, Jing-Tang Lin, Po-Hang Tseng, Che-Wei Chang, Chen-Yuan Dong, Renjie Zhou, Zahid Yaqoob, Irmgard Bischofberger, et al.  Single-shot quantitative polarization imaging of complex birefringent structure dynamics.  arXiv preprint arXiv:2106.06118, 2021.
  * [30] Johannes F De Boer, Christoph K Hitzenberger, and Yoshiaki Yasuno.  Polarization sensitive optical coherence tomography–a review.  Biomedical optics express, 8(3):1838–1873, 2017.
  * [31] F Massoumian, R Juškaitis, MAA Neil, and T Wilson.  Quantitative polarized light microscopy.  Journal of Microscopy, 209(1):13–22, 2003.
  * [32] Rudolf Oldenbourg.  Polarized light microscopy: principles and practice.  Cold Spring Harbor Protocols, 2013(11):pdb–top078600, 2013.
  * [33] Shalin B Mehta, Michael Shribak, and Rudolf Oldenbourg.  Polarized light imaging of birefringence and diattenuation at high resolution and high sensitivity.  Journal of Optics, 15(9):094007, 2013.
  * [34] Seungwoo Shin, KyeoReh Lee, Zahid Yaqoob, Peter TC So, and YongKeun Park.  Reference-free polarization-sensitive quantitative phase imaging using single-point optical phase conjugation.  Optics express, 26(21):26858–26865, 2018.
  * [35] Bijie Bai, Hongda Wang, Tairan Liu, Yair Rivenson, John FitzGerald, and Aydogan Ozcan.  Pathological crystal imaging with single-shot computational polarized light microscopy.  Journal of Biophotonics, 13(1):e201960036, 2020.
  * [36] Seungri Song, Jeongsoo Kim, Sunwoong Hur, Jaewoo Song, and Chulmin Joo.  Large-area, high-resolution birefringence imaging with polarization-sensitive fourier ptychographic microscopy.  ACS Photonics, 8(1):158–165, 2021.
  * [37] Qinghua Song, Arthur Baroni, Rajath Sawant, Peinan Ni, Virginie Brandli, Sébastien Chenot, Stéphane Vézian, Benjamin Damilano, Philippe de Mierry, Samira Khadir, et al.  Ptychography retrieval of fully polarized holograms from geometric-phase metasurfaces.  Nature communications, 11(1):1–8, 2020.
  * [38] Li-Hao Yeh, Ivan E Ivanov, Syuan-Ming Guo, Bryant B Chhun, Ezzat Hashemi, May H Han, and Shalin B Mehta.  upti: uniaxial permittivity tensor imaging of intrinsic density and anisotropy.  In Novel Techniques in Microscopy, pages NM3C–4. Optical Society of America, 2021.
  * [39] Patrick Ferrand, Arthur Baroni, Marc Allain, and Virginie Chamard.  Quantitative imaging of anisotropic material properties with vectorial ptychography.  Optics letters, 43(4):763–766, 2018.
  * [40] Roarke Horstmeyer and Changhuei Yang.  A phase space model of fourier ptychographic microscopy.  Optics express, 22(1):338–358, 2014.
  * [41] Pavan Chandra Konda, Lars Loetgering, Kevin C Zhou, Shiqi Xu, Andrew R Harvey, and Roarke Horstmeyer.  Fourier ptychography: current applications and future promises.  Optics Express, 28(7):9603–9630, 2020.
  * [42] Xiaoze Ou, Guoan Zheng, and Changhuei Yang.  Embedded pupil function recovery for fourier ptychographic microscopy.  Opt. Express, 22(5):4960–4972, Mar 2014.
  * [43] Lei Tian, Xiao Li, Kannan Ramchandran, and Laura Waller.  Multiplexed coded illumination for fourier ptychography with an led array microscope.  Biomedical optics express, 5(7):2376–2389, 2014.
  * [44] Andrew Maiden, Daniel Johnson, and Peng Li.  Further improvements to the ptychographical iterative engine.  Optica, 4(7):736–745, 2017.
  * [45] Ken Kreutz-Delgado.  The complex gradient operator and the cr-calculus.  arXiv preprint arXiv:0906.4835, 2009.
  * [46] Xiaoze Ou, Roarke Horstmeyer, Guoan Zheng, and Changhuei Yang.  High numerical aperture fourier ptychography: principle, implementation and characterization.  Opt. Express, 23(3):3472–3491, Feb 2015.
  * [47] Valery V Tuchin.  Polarized light interaction with tissues.  Journal of biomedical optics, 21(7):071114, 2016.
  * [48] Bin Yang, Po-Yi Lee, Yi Hua, Bryn Brazile, Susannah Waxman, Fengting Ji, Ziyi Zhu, and Ian A Sigal.  Instant polarized light microscopy for imaging collagen microarchitecture and dynamics.  Journal of Biophotonics, 14(2):e202000326, 2021.
  * [49] Rudolf Oldenbourg.  Analysis of edge birefringence.  Biophysical journal, 60(3):629–641, 1991.
  * [50] Alexander J Howie, Douglas B Brewer, Daniel Howell, and Adrian P Jones.  Physical basis of colors seen in congo red-stained amyloid in polarized light.  Laboratory investigation, 88(3):232–242, 2008.
  * [51] James P McGuire and Russell A Chipman.  Polarization aberrations. 1. rotationally symmetric optical systems.  Applied optics, 33(22):5080–5100, 1994.
  * [52] Lars Loetgering, Xiaomeng Liu, Anne CC De Beurs, Mengqi Du, Guido Kuijper, Kjeld SE Eikema, and Stefan Witte.  Tailoring spatial entropy in extreme ultraviolet focused beams for multispectral ptychography.  Optica, 8(2):130–138, 2021.
  * [53] Pavan Chandra Konda, Jonathan M Taylor, and Andrew R Harvey.  Multi-aperture fourier ptychographic microscopy, theory and validation.  Optics and Lasers in Engineering, 138:106410, 2021.
  * [54] Kyung Chul Lee, Kyungwon Lee, Jaewoo Jung, Se Hee Lee, Donghyun Kim, and Seung Ah Lee.  A smartphone-based fourier ptychographic microscope using the display screen for illumination.  ACS Photonics, 8(5):1307–1315, 2021.
  * [55] Li-Hao Yeh, Jonathan Dong, Jingshan Zhong, Lei Tian, Michael Chen, Gongguo Tang, Mahdi Soltanolkotabi, and Laura Waller.  Experimental robustness of fourier ptychography phase retrieval algorithms.  Optics express, 23(26):33214–33240, 2015.
  * [56] Jianliang Qian, Chao Yang, A Schirotzek, F Maia, and S Marchesini.  Efficient algorithms for ptychographic phase retrieval.  Inverse Problems and Applications, Contemp. Math, 615:261–280, 2014\.
  * [57] John M Rodenburg and Helen ML Faulkner.  A phase retrieval algorithm for shifting illumination.  Applied physics letters, 85(20):4795–4797, 2004.

##  Appendix A A Sequential Gauss-Newton-based vectorial Fourier ptychography
solver

Using the configuration discussed in the Measurement strategy section, the
general vectorial Foruier ptychography(vFP) forward model can be simplified as

$\displaystyle{i_{1i}(\mathbf{r})}$ $\displaystyle=\lvert
z_{1i}(\mathbf{r})\rvert^{2}$ (9a) $\displaystyle{i_{2i}(\mathbf{r})}$
$\displaystyle=\lvert z_{2i}(\mathbf{r})\rvert^{2}$ (9b)
$\displaystyle{i_{3i}(\mathbf{r})}$ $\displaystyle=\lvert
z_{3i}(\mathbf{r})\rvert^{2}$ (9c) $\displaystyle{i_{4i}(\mathbf{r})}$
$\displaystyle=\lvert z_{4i}(\mathbf{r})\rvert^{2},$ (9d)

where

$\displaystyle z_{1i}(\mathbf{r})=\mathcal{F}^{-1}\\{P_{xx}(\mathbf{u})$
$\displaystyle\cdot O_{xx}(\mathbf{u}-\mathbf{u}_{i})+P_{xy}(\mathbf{u})\cdot
O_{yx}(\mathbf{u}-\mathbf{u}_{i})\\}$ (10a) $\displaystyle
z_{2i}(\mathbf{r})=\mathcal{F}^{-1}\\{P_{yx}(\mathbf{u})$ $\displaystyle\cdot
O_{xy}(\mathbf{u}-\mathbf{u}_{i})+P_{yy}(\mathbf{u})\cdot
O_{yy}(\mathbf{u}-\mathbf{u}_{i})\\}$ (10b) $\displaystyle
z_{3i}(\mathbf{r})=\mathcal{F}^{-1}\\{P_{xx}(\mathbf{u})$ $\displaystyle\cdot
O_{xx}(\mathbf{u}-\mathbf{u}_{i})+P_{xx}(\mathbf{u})\cdot
O_{xy}(\mathbf{u}-\mathbf{u}_{i})$ $\displaystyle+$ $\displaystyle
P_{xy}(\mathbf{u})\cdot
O_{yx}(\mathbf{u}-\mathbf{u}_{i})+P_{xy}(\mathbf{u})\cdot
O_{yy}(\mathbf{u}-\mathbf{u}_{i})\\}$ (10c) $\displaystyle
z_{4i}(\mathbf{r})=\mathcal{F}^{-1}\\{P_{yx}(\mathbf{u})$ $\displaystyle\cdot
O_{xx}(\mathbf{u}-\mathbf{u}_{i})+P_{yx}(\mathbf{u})\cdot
O_{xy}(\mathbf{u}-\mathbf{u}_{i})$ $\displaystyle+$ $\displaystyle
P_{yy}(\mathbf{u})\cdot
O_{yx}(\mathbf{u}-\mathbf{u}_{i})+P_{yy}(\mathbf{u})\cdot
O_{yy}(\mathbf{u}-\mathbf{u}_{i})\\}.$ (10d)

Assuming everything is Nyquist rate sampled and stored as discrete signal,
using vectorization notation, we can write the object as
$\mathbf{x}=[\mathbf{x}_{1};\mathbf{x}_{2};\mathbf{x}_{3};\mathbf{x}_{4}];\mathbf{x}\in\mathbb{C}^{4m^{2}}$,
where
$\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\mathbf{x}_{4}\in\mathbb{C}^{m^{2}}$
are vectorization notations for
$O_{xy}(\mathbf{u}),O_{xy}(\mathbf{u}),O_{xy}(\mathbf{u}),O_{xy}(\mathbf{u})$.
$[;]$ denotes row concatenation. Similarly, we write measurements as
$\mathbf{y}_{i}=[\mathbf{y}_{1i};\mathbf{y}_{2i};\mathbf{y}_{3i};\mathbf{y}_{4i}];\mathbf{y}_{i}\in\mathbb{C}^{n^{2}}$,
where
$\mathbf{y}_{1i},\mathbf{y}_{1i},\mathbf{y}_{1i},\mathbf{y}_{1i}\in\mathbb{C}^{n^{2}}$
are vectorization notation for measurements
$i_{1i}(\mathbf{r}),i_{2i}(\mathbf{r}),i_{3i}(\mathbf{r}),i_{4i}(\mathbf{r})$,
and
$\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}\in\mathbb{C}^{4n^{2}}$
are vectors representing pupils $P_{xx},P_{xy},P_{yx},P_{yy}$. $m^{2}$ and
$n^{2}$ are pixel size of the object and camera. We can then rewrite the
forward model as

$\mathbf{y}_{i}=\lvert\mathbf{H}\mathbf{P}\mathbf{D}_{i}\hat{\mathbf{x}}\rvert,$
(11)

where $\hat{\mathbf{x}}$ is the ground truth of $\mathbf{x}$, $\mathbf{D}_{i}$
is composed of four truncating matrix
$\mathbf{Q}_{i}\in\mathbb{B}^{n^{2}\times m^{2}},\mathbb{B}=\\{0,1\\},$ that
crops and keeps the low-passed frequency region for the $i^{th}$ illumination,
and $\mathbf{H}$ contains four 2D inverse Fourier transform matrix
$\mathbf{F}\in\mathbb{C}^{n^{2}\times n^{2}}$,

$\mathbf{D}_{i}=\begin{bmatrix}&\mathbf{Q}_{i}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{Q}_{i}&\mathbf{0}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{0}&\mathbf{Q}_{i}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{Q}_{i}\end{bmatrix},\mathbf{H}=\begin{bmatrix}&\mathbf{F}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{F}&\mathbf{0}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{0}&\mathbf{F}&\mathbf{0}\\\
&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{F}\end{bmatrix};$ (12)

$\mathbf{0}$ is a all-zero matrix. $\mathbf{P}$ is the pupil matrix

$\mathbf{P}=\begin{bmatrix}&\text{Diag}(\mathbf{p}_{1})&\mathbf{0}&\text{Diag}(\mathbf{p}_{2})&\mathbf{0}\\\
&\mathbf{0}&\text{Diag}(\mathbf{p}_{3})&\mathbf{0}&\text{Diag}(\mathbf{p}_{4})\\\
&\text{Diag}(\mathbf{p}_{1})&\text{Diag}(\mathbf{p}_{1})&\text{Diag}(\mathbf{p}_{2})&\text{Diag}(\mathbf{p}_{2})\\\
&\text{Diag}(\mathbf{p}_{3})&\text{Diag}(\mathbf{p}_{3})&\text{Diag}(\mathbf{p}_{4})&\text{Diag}(\mathbf{p}_{4})\end{bmatrix},$
(13)

where $\text{Diag}(\mathbf{a})$ represents a diagonal matrix with the vector
$\mathbf{a}$ on its diagonal. We aim to recover $\hat{\mathbf{x}}$ from a
sequence of $L$ measurements $y_{i}$ from $i^{th}$ LED, by minimizing the
Euclidean distance between predicted and measured amplitude, as an amplitude-
based loss gives empirically improved results comparing to intensity-based
loss[55],

$f(\mathbf{x})=\sum_{i=1}^{L}\lVert\lvert\mathbf{z}_{i}\rvert-\mathbf{y}_{i}\rVert_{2}^{2},$
(14)

where

$\mathbf{z}_{i}\coloneqq\mathbf{H}\mathbf{P}\mathbf{D}_{i}\mathbf{x}$ (15)

is the measurement with phase, and $L$ is the total LED numbers. We attempt to
minimize this quadratic objective function using a sequential Gauss-Newton-
based method, which minimize a part of the loss function
$f_{i}(\mathbf{x})=\frac{1}{2}\lVert\lvert\mathbf{z}_{i}\rvert-\mathbf{y}_{i}\rVert_{2}^{2}$
in sequence. Using the $\mathbb{C}\mathbb{R}$ \- calculus[45] that treats
$\mathbf{x}$ and its complex conjugate $\bar{\mathbf{x}}$ as independent
variables, we can write out the Jacobian as 111note that we drop the
$\frac{1}{2}$ in eq. 14 that is commonly used to cancel the
”$\lVert\cdot\rVert^{2}$” introduced by the quadratic term, thus, the Jacobian
in eq.16 does not have to have the $\frac{1}{2}$ due to the derivative of
”$\sqrt{\cdot}$”

$\mathbf{J}_{i}=\frac{\partial\lvert\mathbf{z}_{i}\rvert}{\partial\mathbf{x}}=\text{Diag}(\overline{\mathbf{z}_{i}}/\lvert\mathbf{z}_{i}\rvert)\mathbf{H}\mathbf{P}\mathbf{D}_{i},.$
(16)

’$/$’ is the element-wise division, and the Hessian can be approximated with

$\mathbf{J}_{i}^{\text{H}}\mathbf{J}_{i}=\mathbf{D}_{i}^{\text{H}}\mathbf{P}^{\text{H}}\mathbf{P}\mathbf{D}_{i}.$
(17)

$(\cdot)^{\text{H}}$ is the Hermitian adjoint. Starting with $\mathbf{x}^{0}$,
the iteration proceeds as

$\mathbf{x}^{t,i+1}=\mathbf{x}^{t,i}-\gamma\big{(}\mathbf{J}_{i}^{\text{H}}\mathbf{J}_{i}\big{)}^{-1}\mathbf{J}_{i}^{\text{H}}(\lvert\mathbf{z}_{i}(\mathbf{x}^{t,i})\rvert-\mathbf{y}_{i}).$
(18)

Moreover,

$\displaystyle\big{(}\mathbf{J}_{i}^{\text{H}}\mathbf{J}_{i}\big{)}^{-1}\mathbf{J}_{i}^{\text{H}}=$
$\displaystyle\mathbf{D}_{i}^{\text{H}}\big{(}\mathbf{P}^{\text{H}}\mathbf{P}\big{)}^{-1}\mathbf{D}_{i}^{\text{H}}\mathbf{D}_{i}\mathbf{P}^{\text{H}}\mathbf{H}^{\text{H}}\text{Diag}(\mathbf{z}_{i}/\lvert\mathbf{z}_{i}\rvert)$
(19a) $\displaystyle=$
$\displaystyle\mathbf{D}_{i}^{\text{H}}\big{(}\mathbf{P}^{\text{H}}\mathbf{P}\big{)}^{-1}\mathbf{P}^{\text{H}}\mathbf{H}^{\text{H}}\text{Diag}(\mathbf{z}_{i}/\lvert\mathbf{z}_{i}\rvert)$
(19b) $\displaystyle=$
$\displaystyle\mathbf{D}_{i}^{\text{H}}\mathbf{P}^{\dagger}\mathbf{H}^{\text{H}}\text{Diag}(\mathbf{z}_{i}/\lvert\mathbf{z}_{i}\rvert),$
(19c)

where

$\displaystyle\mathbf{P}^{\dagger}$
$\displaystyle\coloneqq\big{(}\mathbf{P}^{\text{H}}\mathbf{P}\big{)}^{-1}\mathbf{P}$
(20)
$\displaystyle=\begin{bmatrix}&\text{Diag}(\mathbf{p}_{4}/\mathbf{q})&\text{Diag}(\mathbf{p}_{2}/\mathbf{q})&\mathbf{0}&\text{Diag}(\mathbf{p}_{2}/\mathbf{q})\\\
&\text{Diag}(-\mathbf{p}_{4}/\mathbf{q})&\text{Diag}(\mathbf{p}_{3}/\mathbf{q})&\text{Diag}(\mathbf{p}_{4}/\mathbf{q})&\mathbf{0}\\\
&\text{Diag}(\mathbf{p}_{1}/\mathbf{q})&\text{Diag}(\mathbf{p}_{1}/\mathbf{q})&\text{Diag}(\mathbf{p}_{2}/\mathbf{q})&\text{Diag}(\mathbf{p}_{2}/\mathbf{q})\\\
&\text{Diag}(\mathbf{p}_{3}/\mathbf{q})&\text{Diag}(\mathbf{p}_{3}/\mathbf{q})&\text{Diag}(\mathbf{p}_{4}/\mathbf{q})&\text{Diag}(\mathbf{p}_{4}/\mathbf{q})\end{bmatrix},$
(21)

$\mathbf{q}=\mathbf{p}_{1}\cdot\mathbf{p}_{4}-\mathbf{p}_{2}\cdot\mathbf{p}_{3}$,
and again, ’$(/)$’ and ’$(\cdot)$’ are the element-wise division and
multiplication, respectively. It’s worth mentioning that
$\mathbf{P}^{\dagger}$ is the left inverse of the $\mathbf{P}$ for the benefit
of readers who are new to the field; hence, historically it is also widely
referred as an alternating projection method(see discussions in [56]).
Moreover, the adjoint of the cropping matrix, $\mathbf{D}_{i}^{\text{H}}$, is
a zero-padding matrix. Thus, the update in 18 can be explicitly implemented as

$\displaystyle O_{xx}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{xx}(\mathbf{u}-\mathbf{u}_{i})-Q^{*}(\mathbf{u})\Psi_{1i}(\mathbf{u})/\lvert
Q(\mathbf{u})\rvert^{2}$ (22a) $\displaystyle
O_{xy}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{xy}(\mathbf{u}-\mathbf{u}_{i})-Q^{*}(\mathbf{u})\Psi_{2i}(\mathbf{u})/\lvert
Q(\mathbf{u})\rvert^{2}$ (22b) $\displaystyle
O_{yx}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{yx}(\mathbf{u}-\mathbf{u}_{i})-Q^{*}(\mathbf{u})\Psi_{3i}(\mathbf{u})/\lvert
Q(\mathbf{u})\rvert^{2}$ (22c) $\displaystyle
O_{yy}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{yy}(\mathbf{u}-\mathbf{u}_{i})-Q^{*}(\mathbf{u})\Psi_{4i}(\mathbf{u})/\lvert
Q(\mathbf{u})\rvert^{2},$ (22d)

where

$Q(\mathbf{u})=P_{xx}(\mathbf{u})P_{yy}(\mathbf{u})-P_{xy}(\mathbf{u})P_{yx}(\mathbf{u}),$
(23)

$\displaystyle\Psi_{1}(\mathbf{u})$
$\displaystyle=P_{yy}(\mathbf{u})\Phi_{1i}(\mathbf{u})+P_{xy}(\mathbf{u})\Phi_{2i}(\mathbf{u})+P_{xy}(\mathbf{u})\Phi_{4i}(\mathbf{u})$
(24a) $\displaystyle\Psi_{2}(\mathbf{u})$
$\displaystyle=-P_{yy}(\mathbf{u})\Phi_{1i}(\mathbf{u})+P_{yx}(\mathbf{u})\Phi_{2i}(\mathbf{u})+P_{yy}(\mathbf{u})\Phi_{3i}(\mathbf{u})$
(24b) $\displaystyle\Psi_{3}(\mathbf{u})$
$\displaystyle=P_{xx}(\mathbf{u})\Phi_{1i}(\mathbf{u})+P_{xx}(\mathbf{u})\Phi_{2i}(\mathbf{u})+P_{xy}(\mathbf{u})\Phi_{3i}(\mathbf{u})+P_{xy}(\mathbf{u})\Phi_{4i}(\mathbf{u})$
(24c) $\displaystyle\Psi_{4}(\mathbf{u})$
$\displaystyle=P_{yx}(\mathbf{u})\Phi_{1i}(\mathbf{u})+P_{yx}(\mathbf{u})\Phi_{2i}(\mathbf{u})+P_{yy}(\mathbf{u})\Phi_{3i}(\mathbf{u})+P_{yy}(\mathbf{u})\Phi_{4i}(\mathbf{u}),$
(24d)

and

$\Phi_{mi}(\mathbf{u})=\mathcal{F}\Big{\\{}\Big{(}\mathbf{z}_{mi}(\mathbf{r})-\frac{\mathbf{z}_{mi}(\mathbf{r})}{\lvert\mathbf{z}_{mi}(\mathbf{r})\rvert}\sqrt{i_{mi}(\mathbf{r})}\Big{)}\Big{\\}},m=1,2,3,4$
(25)

respectively. In addition, we introduce a $\lvert Q\rvert_{\text{max}}^{2}$
regularization to the division, which was shown to be robust and stable[57,
44], and the final proposed update becomes

$\displaystyle O_{xx}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{xx}(\mathbf{u}-\mathbf{u}_{i})-\frac{Q^{*}(\mathbf{u})\Psi_{1i}(\mathbf{u})}{(1-\alpha)\lvert
Q(\mathbf{u})\rvert^{2}+\alpha\lvert Q(\mathbf{u})\rvert_{\text{max}}^{2}}$
(26a) $\displaystyle O_{xy}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{xy}(\mathbf{u}-\mathbf{u}_{i})-\frac{Q^{*}(\mathbf{u})\Psi_{2i}(\mathbf{u})}{(1-\alpha)\lvert
Q(\mathbf{u})\rvert^{2}+\alpha\lvert Q(\mathbf{u})\rvert_{\text{max}}^{2}}$
(26b) $\displaystyle O_{yx}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{yx}(\mathbf{u}-\mathbf{u}_{i})-\frac{Q^{*}(\mathbf{u})\Psi_{3i}(\mathbf{u})}{(1-\alpha)\lvert
Q(\mathbf{u})\rvert^{2}+\alpha\lvert Q(\mathbf{u})\rvert_{\text{max}}^{2}}$
(26c) $\displaystyle O_{yy}(\mathbf{u}-\mathbf{u}_{i})$
$\displaystyle=O_{yy}(\mathbf{u}-\mathbf{u}_{i})-\frac{Q^{*}(\mathbf{u})\Psi_{4i}(\mathbf{u})}{(1-\alpha)\lvert
Q(\mathbf{u})\rvert^{2}+\alpha\lvert Q(\mathbf{u})\rvert_{\text{max}}^{2}}.$
(26d)

Similarly, we can derive the update for the pupils

$\displaystyle P_{xx}(\mathbf{u})$
$\displaystyle=P_{xx}(\mathbf{u})-\frac{T_{i}^{*}(\mathbf{u})\Lambda_{1i}(\mathbf{u})}{(1-\alpha)\lvert
T_{i}(\mathbf{u})\rvert^{2}+\alpha\lvert
T_{i}(\mathbf{u})\rvert_{\text{max}}^{2}}$ (27a) $\displaystyle
P_{xy}(\mathbf{u})$
$\displaystyle=P_{xy}(\mathbf{u})-\frac{T_{i}^{*}(\mathbf{u})\Lambda_{2i}(\mathbf{u})}{(1-\alpha)\lvert
T_{i}(\mathbf{u})\rvert^{2}+\alpha\lvert
T_{i}(\mathbf{u})\rvert_{\text{max}}^{2}}$ (27b) $\displaystyle
P_{yx}(\mathbf{u})$
$\displaystyle=P_{yx}(\mathbf{u})-\frac{T_{i}^{*}(\mathbf{u})\Lambda_{3i}(\mathbf{u})}{(1-\alpha)\lvert
T_{i}(\mathbf{u})\rvert^{2}+\alpha\lvert
T_{i}(\mathbf{u})\rvert_{\text{max}}^{2}}$ (27c) $\displaystyle
P_{yy}(\mathbf{u})$
$\displaystyle=P_{yy}(\mathbf{u})-\frac{T_{i}^{*}(\mathbf{u})\Lambda_{4i}(\mathbf{u})}{(1-\alpha)\lvert
T_{i}(\mathbf{u})\rvert^{2}+\alpha\lvert
T_{i}(\mathbf{u})\rvert_{\text{max}}^{2}},$ (27d)

where

$T_{i}=O_{xx}(\mathbf{u}-\mathbf{u}_{i})O_{yy}(\mathbf{u}-\mathbf{u}_{i})-O_{xy}(\mathbf{u}-\mathbf{u}_{i})O_{yx}(\mathbf{u}-\mathbf{u}_{i}),$
(28)

and

$\displaystyle\Lambda_{1i}(\mathbf{u})$
$\displaystyle=O_{yx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{2i}(\mathbf{u})+O_{yy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{2i}(\mathbf{u})-O_{yx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{3i}(\mathbf{u})$
(29a) $\displaystyle\Lambda_{2i}(\mathbf{u})$
$\displaystyle=-O_{xx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{2i}(\mathbf{u})-O_{xy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{2i}(\mathbf{u})+O_{xx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{3i}(\mathbf{u})$
(29b) $\displaystyle\Lambda_{3i}(\mathbf{u})$
$\displaystyle=-O_{yx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{1i}(\mathbf{u})-O_{yy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{1i}(\mathbf{u})+O_{yy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{4i}(\mathbf{u})$
(29c) $\displaystyle\Lambda_{4i}(\mathbf{u})$
$\displaystyle=O_{xx}(\mathbf{u}-\mathbf{u}_{i})\Phi_{1i}(\mathbf{u})+O_{xy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{1i}(\mathbf{u})-O_{xy}(\mathbf{u}-\mathbf{u}_{i})\Phi_{4i}(\mathbf{u})$
(29d)

##  Appendix B Analyzing polarization state with a linear polarizer

Here we derive a simple way to compute the detected intensity value after the
electric field $\mathbf{e}=[E_{x},E_{y}]^{T}$ passing through the $l^{th}$
analyzer, made with a linear polarizer placed at angle $\alpha_{l}$ with
respect to the x-polarization. This handy equation is used in the eq(3) in the
main text. The Jones matrix for the analyzer is [1]

$\mathbf{M}_{\alpha_{l}}=\begin{bmatrix}\text{cos}^{2}(\alpha_{l})&\text{cos}(\alpha_{l})\text{sin}(\alpha_{l})\\\
\text{cos}(\alpha_{l})\text{sin}(\alpha_{l})&\text{sin}^{2}(\alpha_{l})\end{bmatrix}.$
(30)

Hence, the electric field exiting the polarizer is

$\displaystyle\mathbf{e}^{\prime}=\mathbf{M}_{\alpha_{l}}\mathbf{e}$ (31)
$\displaystyle\mkern-36.0mu=\begin{bmatrix}\text{cos}^{2}(\alpha_{l})E_{x}+\text{cos}(\alpha_{l})\text{sin}(\alpha_{l})E_{y}\\\
\text{cos}(\alpha_{l})\text{sin}(\alpha_{l})E_{x}+\text{sin}^{2}(\alpha_{l})E_{y}\end{bmatrix}.$
(32)

Similarly, we denote the x-y polarization components of the existing electric
field $\mathbf{e}^{\prime}=[E_{x}^{\prime},E_{y}^{\prime}]^{T}$. After a bit
tedious math, the detected intensity can be calculated as

$\displaystyle i_{\alpha_{l}}=$
$\displaystyle\lvert\mathbf{e}^{\prime}\rvert^{2}$ (33) $\displaystyle=$
$\displaystyle\lvert E_{x}^{\prime}\rvert^{2}+\lvert E_{y}^{\prime}\rvert^{2}$
(34) $\displaystyle=$
$\displaystyle\lvert\text{cos}\alpha_{l}E_{x}+\text{sin}\alpha_{l}E_{y}\rvert^{2}$
(35) $\displaystyle=$
$\displaystyle\lvert\textbf{a}_{l}^{T}\textbf{e}\rvert^{2},$ (36)

where $\mathbf{a}_{l}=[\text{cos}(\alpha_{l}),\text{sin}(\alpha_{l})]^{T}$. In
addition, an alternative intuitive graphic explanation is also provided as
fig. 6

Figure 6: Graphic explaination of eq.36. The slit is the linear polarizer.
$\alpha$ is the angle between the linear polarizer and the x-polarization.
Projecting the $E_{x}$ and $E_{y}$ components of the incoming electric field
to the polarizer direction yields a passing electric field
$E_{x}\text{cos}(\alpha_{l})+E_{y}\text{sin}(\alpha_{l})$


# Evaluation of CNN-based Automatic Music Tagging Models

###### Abstract

Recent advances in deep learning accelerated the development of content-based
automatic music tagging systems. Music information retrieval (MIR) researchers
proposed various architecture designs, mainly based on convolutional neural
networks (CNNs), that achieve state-of-the-art results in this multi-label
binary classification task. However, due to the differences in experimental
setups followed by researchers, such as using different dataset splits and
software versions for evaluation, it is difficult to compare the proposed
architectures directly with each other. To facilitate further research, in
this paper we conduct a consistent evaluation of different music tagging
models on three datasets (MagnaTagATune, Million Song Dataset, and MTG-
Jamendo) and provide reference results using common evaluation metrics (ROC-
AUC and PR-AUC). Furthermore, all the models are evaluated with perturbed
inputs to investigate the generalization capabilities concerning time stretch,
pitch shift, dynamic range compression, and addition of white noise. For
reproducibility, we provide the PyTorch implementations with the pre-trained
models.

##  1 Introduction

Automatic music tagging is a multi-label binary classification task that aims
to predict relevant tags for a given song. Typically these tags carry useful
semantic music information that can later be used for other applications such
as music recommendation or music retrieval. To tackle the problem of music
tagging, recent studies in music information retrieval (MIR) adopted deep
neural networks, mostly based on convolutional neural networks (CNNs),
motivated by their huge success in other domains (e.g., computer vision,
natural language processing). The introduction of deep learning helped to
break the previous glass ceiling in the performance of music tagging systems
and MIR researchers started actively proposing ingenious architecture designs.
As a result, the hand-crafted feature-based approaches were replaced by data-
driven feature learning approaches in most recent automatic music tagging
research.

To maximize the advantages of CNNs, a deep fully convolutional network (FCN)
was proposed for music tagging [1]. It uses a stack of $3\times 3$ rectangular
filters followed by a max-pooling layer. As an alternative, the Musicnn [2],
also a Mel-spectrogram-based CNN, tried to incorporate domain knowledge into
its filter designs so that the model can capture timbral characteristics and
temporal patterns using vertical filters and horizontal filters, respectively.
Sample-level CNN [3] pursued an assumption-free end-to-end model by applying
1D convolution directly to a raw audio waveform, and the following research
[4] improved the performance by adding squeeze-and-excitation blocks [5].
Different from images, however, music is sequential. For this reason, a
convolutional recurrent neural network (CRNN) [6] was proposed to extract
local features using CNNs and summarize them with recurrent neural networks
(RNNs). Another sequence modeling approach [7] adapted the self-attention
mechanism [8] to summarize the temporal sequence of the extracted local
features by CNNs. Finally, Harmonic CNN [9] used a harmonically stacked
trainable representation to preserve spectro-temporal locality in convolution
layers.

Unfortunately, due to different experimental setups followed by the authors of
these approaches when reporting results (e.g., dataset splits, library
versions, computing environments, and optimization methods), it is difficult
to compare the proposed architectures directly with each other. In this paper,
we address this issue and report experimental results for various state-of-
the-art music tagging models using three different datasets (MagnaTagATune,
Million Song Dataset, and MTG-Jamendo dataset) with a consistent experimental
setup. In addition, we conduct experiments to assess the robustness of these
architectures against four different types of deformations [10] and determine
their generalization abilities. For the reproducibility, we provide PyTorch
implementations for all the models considered and their pre-trained
models.111https://github.com/minzwon/sota-music-tagging-models/

The paper is organized as follows. Section 2 describes automatic music tagging
tasks with details of multiple instance problem, evaluation metrics, and
dataset descriptions. Section 3 elaborates on the model designs. We report
model performances and their generalization capabilities in Section 4 and
Section 5, respectively. Finally, Section 6 concludes the paper.

##  2 Automatic Music Tagging

###  2.1 Multiple Instance Problem

Various semantic information characterizing music, such as genres, subgenres,
moods, instruments, decades, and languages, can be expressed in the form of
music tags. Automatic content-based music tagging is a task that aims to
predict such relevant tags for a given song based on its acoustic
characteristics. However, depending on the tag, relevant characteristics are
not necessarily predominant in the entire music recording. For example, when a
song has a tag female vocal, it does not imply that the female vocal appears
in every time segment of the song. In essence, this is a multiple instance
problem [11]. A song can have many acoustic characteristics (instances) that
describe it, but often only specific characteristics (instances) are
responsible for the associated music tag. In most cases, we do not have time
precise instance-level annotations because the precise labeling can be
laborsome, therefore a music tag associated with a song is simply applied to
all music excerpts (instances) of the song during training [12].

There are two approaches to handle the multiple instance problem in music
tagging. One is to train the model on full songs and produce song-level
predictions from a song-level input. From a given song-level input, the model
has to predict relevant tags. Another one is to train the model on short audio
chunks (instances) and generate chunk-level predictions, which can be later
aggregated (e.g., majority vote, average) in the evaluation phase. In this
paper, we refer to these two approaches as song-level training and chunk-level
training, respectively.

###  2.2 Evaluation Metrics

A common evaluation metric of binary decision problems is the area under
receiver operating characteristic curve (ROC-AUC). However, the area under
precision-recall curve (PR-AUC) can be more informative for evaluation on
highly skewed datasets [13]. Hence, we use both macro ROC-AUC and PR-AUC to
evaluate all considered music tagging models with both scores being averaged
across all the tags the models operate with. For our robustness studies, we
report the variance of ROC-AUC and PR-AUC scores obtained on perturbed audio
inputs.

###  2.3 Datasets

MagnaTagATune (MTAT) [14] is one of the most commonly used datasets for
benchmarking automatic music tagging systems. It contains multi-label
annotations by genre, mood, and instrumentation for 25,877 audio segments,
each 30s long. The audio is in the MP3 format (32 Kbps bitrate and 16 kHz
sample rate). Originally the dataset is split into 16 folders, and commonly
the first 12 folders are used for training, the 13th for validation, and the
last three are used for testing. Only 50 most frequent tags are typically used
for the task. The top 50 tags include genre and instrumentation labels, as
well as decades (e.g., ’80s’ and ’90s’) and moods.

In our work, we follow the same data split and we use the top 50 tags to be
consistent with results reported in the majority of previous
studies.222https://github.com/jongpillee/music_dataset_split/tree/master/MTAT_split
However, we have noted that the performances reported in some studies are
inconsistent, as their authors discarded audio segments without any associated
tags (leading to slightly higher values of performance). Also, some of the
previous studies are re-using those reports unintentionally compare
incompatible performance values, which is one of the main motivations of our
work.

Million Song Dataset (MSD) [15] is a dataset of audio features for one million
songs, partially expanded by the MIR community with crowdsourced tags from
Last.fm as well as a mapping to 30s audio preview segments originally obtained
from 7digital.333https://www.7digital.com In total, this subset of the dataset
contains 241,904 annotated song segments and it is commonly used for
benchmarking on a larger scale. The tags cover genre, instrumentation, moods
and decades. The audio segments vary in quality, being encoded as MP3s with a
bitrate from 64 to 128 Kbps and sample rate of 22 kHz or 44 kHz.

Again, in our study we follow the dataset split commonly used by
researchers444https://github.com/jongpillee/music_dataset_split [3, 2] and use
only 50 most frequent tags for consistency with previous studies. This split
includes 201,680 songs for training, 11,774 for validation and 28,435 for
testing. Unfortunately, similarly to MTAT, some previous studies report
inconsistent comparisons of automatic tagging approaches, as there also exists
a slightly different
split,555https://github.com/keunwoochoi/MSD_split_for_tagging containing audio
segments of shorter duration and on which lower performance values were
reported. Note that the tag annotations available for this dataset are
inherently noisy as they come from a free-form social tagging application for
music enthusiasts and are used without any preprocessing intended to improve
the quality of tags [16].

MTG-Jamendo Dataset [17] contains audio for 55,701 full songs and is built
using music publicly available on the Jamendo666https://jamendo.com music
platform under Creative Commons licenses. The minimum duration of each song is
30s, and they are provided in the MP3 (320 Kbps bitrate). Thus, this dataset
contains significantly larger audio segments with higher encoding quality than
MTAT and MSD. The tracks in the dataset are annotated by 692 different tags
covering genres, instrumentation, moods and themes. All tags were originally
provided by the artists submitting music to Jamendo, but they were
preprocessed with the goal of tag cleaning by the creators of the dataset.

Multiple splits of the data are provided for training, validation and test. In
this work we use the split-0 and 50 most frequent
tags.777https://github.com/MTG/mtg-jamendo-dataset As this dataset has been
released only recently, not many studies are reporting the performance of the
models using it yet. Yet, it is a useful addition for the evaluation
methodologies followed by researchers in order to better assess the
generalization of their models.

##  3 Models

Figure 1: Shapes of the first convolution filters and input representations of
different models. (a) FCN, CRNN, self-attention, and short-chunk CNN with a
Mel spectrogram input (b) Musicnn with a Mel spectrogram input (c) sample-
level CNN with a raw audio input (d) Harmonic CNN with a stacked harmonic
tensor.

We here describe the architectures of all models considered in this study.
Essentia [18] and Librosa [19] libraries are used for loading audio files and
extracting Mel spectrograms, respectively. We re-sampled the audio to 16 kHz
sample rate. Table 1 shows different input lengths and the number of Mel
bands. Since some models have a smaller number of Mel bands in their original
implementations, we further experimented with a larger number of Mel bands for
the fair comparison (marked with *). All Mel-spectrogram-based approaches used
512-point FFT with a 50% overlap. For training the models, we used a unified
optimization method: a mixture of scheduled ADAM [20] and stochastic gradient
descent (SGD), introduced in [7]. The best model is selected based on the
validation loss. In the evaluation, chunk-level models average predictions
over 16 chunks to perform the final prediction. Note that FCN and CRNN are
song-level training models for MTAT and MSD but they perform as chunk-level
training models for the MTG-Jamendo dataset since the songs in the dataset are
longer than 29.1s which is the size of receptive fields of both models.

Model | input length | # of Mel | training
---|---|---|---
FCN | 29.1s | 96 | song-level
FCN (128)∗ | 29.1s | 128 | song-level
Musicnn | 3s | 96 | chunk-level
Musicnn (128)∗ | 3s | 128 | chunk-level
Sample-level CNN | 3.69s | - | chunk-level
CRNN | 29.1s | 96 | song-level
CRNN (128)∗ | 29.1s | 128 | song-level
Self-attention | 15s | 128 | chunk-level
Harmonic CNN | 5s | 128 | chunk-level
Short-chunk CNN | 3.69s | 128 | chunk-level
Table 1: Experimental setups of different models.

###  3.1 Fully Convolutional Network

A fully convolutional network (FCN) is a variant of CNN that consists of only
convolutional layers without any fully-connected layers. A FCN for music
tagging [1] uses Mel spectrogram inputs. In the preprocessing step, a 29.1s
audio segment is converted to a $96\times 1366$ Mel spectrogram. It is then
used as an input and is passed through 4 convolutional layers. Each
convolutional layer uses homogeneous $3\times 3$ 2D filters (Figure 1-(a))
followed by a max-pooling layer. Different sizes of strides are used for max-
pooling layers ((2, 4), (4, 5), (3, 8), (4, 8)) to increase the size of
receptive fields to cover the entire input Mel spectrogram ($96\times 1366$).
In the original paper, FCN was trained with a song-level training method since
the track durations in MTAT correspond to the size of the receptive field.
However, this is not the case for MTG-Jamendo containing longer tracks, where
chunk-level (29.1s) training is applied.

###  3.2 Musicnn

The Musicnn [2] model also uses Mel spectrograms as its inputs. Different from
previously proposed models, the architecture design choices in Musicnn rely on
some intuition from the music domain knowledge. The first convolutional layer
of Musicnn consists of vertical and horizontal filters. Vertical filters are
designed to capture pitch-invariant timbral features (bottom-left of Figure
1-(b)): e.g., $38\times 7$ filter captures sub-band information of short
period of time. To enforce the pitch-invariancy, the following max-pooling
layer pools the maximum value across the frequency axis. Horizontal filters,
on the other hand, capture temporal energy envelope of the audio. After the
mean-pooling across the frequency axis of input Mel spectrograms, horizontally
long filters (e.g., $1\times 165$) capture the temporal energy patterns (top-
right of Figure 1-(b)). The extracted timbral and temporal features are
concatenated in the channel, then the following 1D convolutional layers
summarize them to predict relevant tags. Different from FCN, the Musicnn only
uses short audio excerpts (3s) as its inputs during training, i.e., chunk-
level training.

###  3.3 Sample-level CNN

Sample-level CNN [3] tackles the automatic music tagging problem in an end-to-
end fashion. It takes raw audio waveforms as its inputs. Sample-level CNN is
simpler and deeper than Mel spectrogram-based approaches. It consists of ten
1D convolutional layers with $1\times 3$ filters and $1\times 3$ max-poolings
(Figure 1-(c)). Trained front-end filters perform similar to the process of
deriving Mel spectrograms and the back-end convolution layers summarize them.
We also considered a variation of sample-level CNN [4] with squeeze-and-
excitation (SE) [5] blocks. Sample-level CNN and its variant with SE blocks
also use short audio excerpts (3.69s) for the chunk-level training.

###  3.4 Convolutional Recurrent Neural Network

Convolutional recurrent neural network (CRNN) [6] uses Mel spectrogram inputs.
A CRNN can be described as a combination of CNNs and RNNs. The CNN front end
extracts local features and the RNN back end summarizes them temporally. Since
RNNs are more flexible than CNNs for summarizing sequential information, it
can be beneficial to use RNNs for predicting tags that may be affected by
global structures (e.g., moods/themes). Four convolutional layers with
$3\times 3$ 2D filters are used in the front end and two-layer RNNs with gated
recurrent units (GRU) are used in the back end. Long music excerpts (29.1s)
are used as inputs of CRNN.

###  3.5 Self-attention

The self-attention-based music tagging model [7] shares the same intuition as
CRNN to extract local features with CNNs and summarize them with sequence
models. The only difference is that the self-attention mechanism is used
instead of the RNNs for the temporal summarization back end. Motivated by its
huge success in natural language processing [8], the authors adapted the
Transformer encoder, which is a deep stack of self-attention layers, for
automatic music tagging. 15s-long audio excerpts are used for training the
self-attention model.

###  3.6 Harmonic CNN

Harmonic CNN [9] takes advantage of trainable band-pass filters and
harmonically stacked time-frequency representation inputs. Trainable filters
(mainly trainable bandwidths) bring more flexibility to the model. And
harmonically stacked representation preserves spectro-temporal locality while
keeping the harmonic structures through the channel of the input tensor in the
first convolution layer (Figure 1-(d)) as introduced in [21]. The number of
trainable frequency bands is set to 128 and the number of harmonics considered
for stacking is 6. Chunk-level training with 5s audio segments is performed.

###  3.7 Short-chunk CNN

According to the previous work [9], a simple 2D CNN with $3\times 3$ filters
can already claim exceptional results when it is trained with short chunks of
audio, i.e., chunk-level training. It is a very prevalent type of CNN
(sometimes referred as vgg-like) but, to the best of our knowledge, there are
no references for this architecture design in music tagging research. Hence,
we implemented a 7-layer CNN with a fully-connected layer and its extension
with residual connections [22]. Different from FCN, it uses a smaller size of
max-pooling ($2\times 2$) because the input segment is way shorter than the
song-level inputs (29.1s). We used 128 Mel bins so that 7 max-pooling layers
can summarize them into a single dimension ($2^{7}=128$). It uses 3.69s audio
excerpts, hence we call this model “short-chunk CNN” in this paper to
differentiate it from FCN.

##  4 Results

Methods | MTAT | MSD | MTG-Jamendo
---|---|---|---
ROC-AUC | PR-AUC | ROC-AUC | PR-AUC | ROC-AUC | PR-AUC
FCN [1] | 0.9005 | 0.4295 | 0.8744 | 0.2970 | 0.8255 | 0.2801
FCN (with 128 Mel bins) | 0.8994 | 0.4236 | 0.8742 | 0.2963 | 0.8245 | 0.2792
Musicnn [2] | 0.9106 | 0.4493 | 0.8803 | 0.2983 | 0.8226 | 0.2713
Musicnn (with 128 Mel bins) | 0.9092 | 0.4546 | 0.8788 | 3036 | 0.8275 | 0.2810
Sample-level [3] | 0.9058 | 0.4422 | 0.8789 | 0.2959 | 0.8208 | 0.2742
Sample-level + SE [4] | 0.9103 | 0.4520 | 0.8838 | 0.3109 | 0.8233 | 0.2784
CRNN [6] | 0.8722 | 0.3625 | 0.8499 | 0.2469 | 0.7978 | 0.2358
CRNN (with 128 Mel bins) | 0.8703 | 0.3601 | 0.8460 | 0.2330 | 0.7984 | 0.2378
Self-attention [7] | 0.9077 | 0.4445 | 0.8810 | 0.3103 | 0.8261 | 0.2883
Harmonic CNN [9] | 0.9127 | 0.4611 | 0.8898 | 0.3298 | 0.8322 | 0.2956
Short-chunk CNN | 0.9126 | 0.4590 | 0.8883 | 0.3251 | 0.8324 | 0.2976
Short-chunk CNN + Res | 0.9129 | 0.4614 | 0.8898 | 0.3280 | 0.8316 | 0.2951
Table 2: Performances of state-of-the-art models.

We report ROC-AUC and PR-AUC of all implemented models using three datasets in
Table 2. In general, models trained with short audio excerpts (Musicnn,
variants of sample-level CNN, Self-attention, Harmonic CNN, variants of short-
chunk CNN) outperform other models trained with relatively longer audio
segments (FCN, CRNN). Training with short chunks (instances) is noisier: e.g.,
an audio excerpt can have a tag guitar if a guitar appears in the song even
though the selected excerpt doesn’t include guitar sound in it. However, one
can expect a much larger number of examples during the training (e.g., 25,877
tracks $\times$ 16 chunks = 414,032 examples). We suspect this brings the
performance gain when the model is trained with short chunk-level examples.
Furthermore, most of the top 50 tags in the three datasets can be identified
only with a short audio excerpt (e.g., instruments, genres). Thus, the model
does not need a long sequence of audio to perform its binary classification
task. For the top 50 tags in each dataset we experimented with, it is more
beneficial to use chunk-level training with short audio excerpts than the
song-level training.

Short-chunk CNN, short-chunk CNN with residual connections, and Harmonic CNN
showed the best results for every dataset. These three models are trained on
short audio excerpts (3.69s or 5s) and they use $3\times 3$ convolutional
filters followed by $2\times 2$ max-poolings. FCN uses similar filters, but
larger max-poolings which increase its size of the receptive field to fit long
audio segments (29.1s). We conclude that smaller max-poolings with shorter
audio excerpts work better for CNNs with $3\times 3$ filters.

Musicnn shows competitive results in MTAT. However, other models (sample-level
+ SE, self-attention) outperform Musicnn on larger datasets (MSD and MTG-
Jamendo). This confirms an intuition that domain knowledge can be beneficial
for relatively small datasets, reported in [2]. However, the design choices
for parameters of Musicnn restricts the power of the model when it is trained
on larger datasets.

For the sequential models, self-attention outperforms the CRNN. Different from
self-attention mechanisms, RNNs with long sequence inputs suffer from
vanishing gradient problems. Self-attention mechanism alleviates the problems
by providing direct paths between all time steps. According to the reported
visualizations [7], self-attention performs well for pinpointing relevant
short-time acoustic features in the audio sequence, but it was difficult to
determine if the model learned long-time characteristics properly. To
determine such abilities, some tags related to a global structure have to be
cherry-picked and evaluated.

Since FCN, Musicnn, and CRNN use Mel spectrogram inputs with 96 Mel bands,
there can be relative disadvantages when they are compared with other models
using 128 Mel bands. For the fair comparison, we experimented with FCN,
Musicnn, and CRNN using 128 Mel bands. A larger number of Mel bands did not
show any significant impacts on the performances. Since each architecture
design was optimized for a smaller number of Mel bands, simply increasing the
size of input Mel bands cannot guarantee the optimized performance of the
models.

##  5 Robustness Studies

###  5.1 Input Deformations

To further investigate the performance of different state-of-the-art models,
we conducted robustness studies. If a pre-trained model has good
generalization abilities, the prediction of the model should not be sensitive
against small perturbations in the input audio. By applying four different
audio deformations to the test set (pitch shift, time stretch, dynamic range
compression, and addition of white noise), we intended to determine the
generalization abilities of the models. Note that we applied these four
deformations only to the test set, which means that the models have never been
exposed to the same deformations during training. All employed deformations
are based on an existing music data augmentation framework
(MUDA)888https://github.com/bmcfee/muda [10]:

Pitch shift by $n\in\\{-1,1\\}$ semitones.

Time stretch by $\gamma\in\\{2^{-1/2},2^{1/2}\\}$.

Dynamic range compression following speech and music (standard) settings of
Dolby E standards [23].

White noise addition $x_{mixed}=(1-\alpha)\cdot x+\alpha\cdot x_{noise}$ where
$\alpha\in\\{0.1,0.4\\}$.

###  5.2 Robustness results

Figure 2: Evaluations metrics with perturbed audio inputs. Dynamic range
compression is shortened as “drc” in the plot.

Figure 2 shows performances of each model under various input deformations.
Here we tested FCN, Musicnn, sample-level + SE, self-attention, Harmonic CNN,
and short-chunk CNN. We followed the original input preprocessing of each
model because a larger of Mel bands did not show significant effects in
Section 4. CRNN is not included due to its relatively low performance. Dynamic
range compression was the least influential and the white noise addition (0.4)
was the most critical among the four different perturbations considered.
Musicnn is robust against time stretching but it is relatively vulnerable
against pitch shift. We suspect the max pooling layer over frequency axis
hinders the Musicnn from learning generalized representations. Harmonic CNN
and short-chunk CNN were the two best models with original data. However,
Harmonic CNN showed better generalization abilities against input deformations
except for the white noise addition (0.4). Sample-level CNN with SE blocks
showed good performance with a small amount of noise (0.1), but it could not
generalize when this amount was increased (0.4).

##  6 Conclusion

In this paper, we revisit state-of-the-art automatic music tagging models and
report their performances with a consistent experimental setup. In general,
short-chunk-based approaches showed better results than models trained with
larger input segments (FCN, CRNN). The design choices followed by Musicnn
could show good performance on a small dataset, but it restricted the model
from learning more information on larger datasets. Sequential models (CRNN,
self-attention) showed competitive results but could not outperform other
models since most of tags in the datasets do not require long sequences for
their identification. Interestingly, the best performing model is a simple CNN
with $3\times 3$ filters trained on short audio excerpts (short-chunk CNN).
Although the original design choice of the CNN is from computer vision, it
outperformed other methods except for Harmonic CNN.

We further assessed generalization abilities of models by testing perturbed
inputs. We could observe a different ranking of the models in terms of their
performance on each deformation. This implies that the ROC-AUC and PR-AUC
scores are not enough to evaluate music tagging models. In our experiment,
Harmonic CNN and short-chunk CNN consistently report better scores than other
models. Specifically, Harmonic CNN showed the best generalization abilities
against every deformation types except for a heavy white noise addition.

We expect the reported results to be a useful reference for further research
in automatic music tagging. More detailed exploration of the deformations for
robustness tests should be done in the future work. Also, the efficacy of data
augmentation in music tagging has to be further investigated.

Acknowledgments
This work was funded by the predoctoral grant MDM-2015-0502-17-2 from the
Spanish Ministry of Economy and Competitiveness linked to the Maria de Maeztu
Units of Excellence Programme (MDM-2015-0502). Also we acknowledge this
research has been supported by Kakao Corp.

## References

  * [1] K. Choi, G. Fazekas, and M. Sandler, “Automatic tagging using deep convolutional neural networks,” _In Proc. of the 17th International Society for Music Information Retrieval Conference (ISMIR)_ , 2016.
  * [2] J. Pons, O. Nieto, M. Prockup, E. Schmidt, A. Ehmann, and X. Serra, “End-to-end learning for music audio tagging at scale,” _In Proc. of the 19th International Society for Music Information Retrieval Conference (ISMIR)_ , 2018.
  * [3] J. Lee, J. Park, K. L. Kim, and J. Nam, “Sample-level deep convolutional neural networks for music auto-tagging using raw waveforms,” _In Proc. of the 14th Sound and music computing (SMC)_ , 2017.
  * [4] T. Kim, J. Lee, and J. Nam, “Sample-level cnn architectures for music auto-tagging using raw waveforms,” in _Proc. of International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 366–370.
  * [5] J. Hu, L. Shen, and G. Sun, “Squeeze-and-excitation networks,” in _Proc. of the IEEE conference on computer vision and pattern recognition (CVPR)_ , 2018, pp. 7132–7141.
  * [6] K. Choi, G. Fazekas, M. Sandler, and K. Cho, “Convolutional recurrent neural networks for music classification,” in _Proc. of International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 2392–2396.
  * [7] M. Won, S. Chun, and X. Serra, “Toward interpretable music tagging with self-attention,” _arXiv preprint arXiv:1906.04972_ , 2019.
  * [8] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” _Proc. of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 2018.
  * [9] M. Won, S. Chun, , O. Nieto, and X. Serra, “Data-driven harmonic filters for audio representation learning,” _In Proc. of International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020.
  * [10] B. McFee, E. J. Humphrey, and J. P. Bello, “A software framework for musical data augmentation.” in _Proc. of the 16th International Society for Music Information Retrieval Conference (ISMIR)_ , vol. 2015, 2015, pp. 248–254.
  * [11] T. G. Dietterich, R. H. Lathrop, and T. Lozano-Pérez, “Solving the multiple instance problem with axis-parallel rectangles,” _Artificial intelligence_ , vol. 89, no. 1-2, pp. 31–71, 1997.
  * [12] B. McFee, J. Salamon, and J. P. Bello, “Adaptive pooling operators for weakly labeled sound event detection,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 11, pp. 2180–2193, 2018.
  * [13] J. Davis and M. Goadrich, “The relationship between precision-recall and roc curves,” in _Proc. of the 23rd international conference on Machine learning_ , 2006, pp. 233–240.
  * [14] E. Law, K. West, M. I. Mandel, M. Bay, and J. S. Downie, “Evaluation of algorithms using games: The case of music tagging.” in _In Proc. of International Society for Music Information Retrieval Conference (ISMIR)_ , 2009, pp. 387–392.
  * [15] T. Bertin-Mahieux, D. P. Ellis, B. Whitman, and P. Lamere, “The million song dataset,” _In Proc. of the 12th International Society for Music Information Retrieval Conference (ISMIR)_ , 2011.
  * [16] K. Choi, G. Fazekas, K. Cho, and M. Sandler, “The effects of noisy labels on deep convolutional neural networks for music classification,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , 2018.
  * [17] D. Bogdanov, M. Won, P. Tovstogan, A. Porter, and X. Serra, “The mtg-jamendo dataset for automatic music tagging,” _Machine Learning for Music Discovery Workshop, International Conference on Machine Learning (ICML)_ , 2019\.
  * [18] D. Bogdanov, N. Wack, E. Gómez Gutiérrez, S. Gulati, H. Boyer, O. Mayor, G. Roma Trepat, J. Salamon, J. R. Zapata González, X. Serra _et al._ , “Essentia: An audio analysis library for music information retrieval,” in _Proc. of 14th Conference of the International Society for Music Information Retrieval (ISMIR)_ , 2013.
  * [19] B. McFee, C. Raffel, D. Liang, D. P. Ellis, M. McVicar, E. Battenberg, and O. Nieto, “librosa: Audio and music signal analysis in python,” in _Proc. of the 14th python in science conference_ , vol. 8, 2015.
  * [20] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proc. of the International Conference on Learning Representations (ICLR)_ , 2015.
  * [21] R. M. Bittner, B. McFee, J. Salamon, P. Li, and J. P. Bello, “Deep salience representations for f0 estimation in polyphonic music.” in _ISMIR_ , 2017, pp. 63–70.
  * [22] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proc. of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [23] E. Dolby, “Standards and practices for authoring dolby digital and dolby e bitstreams,” _Dolby Labortories, Inc_ , 2002.