# Efficient Modeling of Depolarizing Mueller BRDFs

Quinn Jarecki and Meredith Kupinski

###### Abstract

Light-matter interactions within indoor environments are significantly
depolarizing. Nonetheless, the relatively small polarization attributes are
informative. To make use of this information, polarized-BRDF (pBRDF) models
for common indoor materials are sought. Fresnel reflection and diffuse partial
polarization are popular terms in pBRDF models, but the relative contribution
of each is highly material-dependent and changes based on scattering geometry
and albedo. An efficient pBRDF would describe these dependencies with as few
parameters as possible while retaining physical significance and task-relevant
information. This work compares a triply-degenerate (TD)-Mueller matrix (MM)
model to measurements of 3D printed objects. In this TD-MM model, the
radiometric, polarimetric, and depolarization attributes are decoupled to
reduce the number of parameters. The depolarization is quantified by a single
geometry-dependent parameter, four geometry-independent material constants
describe the polarization properties, and our TD-MM model is normalized to
unit radiance so that the BRDF is decoupled. To test an application of the TD-
MM model the material constants are assumed and the geometry-dependent
depolarization parameter for a red 3D printed sphere is estimated from linear
Stokes images. The geometry-averaged error of the depolarization parameter is
4.2% at 662 nm (high albedo) and 11.7% at 451 nm (low albedo). Since the error
is inversely proportional to albedo and depolarization, the TD-MM model is
referred to as appropriate for depolarization-dominant materials. The
robustness of the TD-MM model is also tested by comparing ground-truth Mueller
images to extrapolations of a red 3D printed Stanford bunny under arbitrary
polarized illumination.

###### Index Terms:

Polarization, partial polarimetry, Mueller matrix imaging, depolarization,
pBRDF

##  I Introduction

As light propagates through a given path connecting the source and observer,
each light-matter interaction can change the orientation of electric field
oscillations. The polarization changes between an incident and exitant ray
direction after interaction are quantified by a polarized bidirectional
reflectance distribution function (pBRDF) [1, 2]. The unpolarized, or
radiometric, scattering behavior of a material is characterized by the scalar-
valued BRDF. A BRDF is the ratio of the differential reflected radiance to the
differential input irradiance and, in general, is a function of the ray
directions [3]. The pBRDF is a Mueller matrix (MM)-valued function of the ray
directions. A MM is a 4$\times$4 real-valued matrix that describes the linear
transformation of incident to exitant Stokes parameters [4, 5, 6].

In high-precision metrology applications, such as defect detection for silicon
wafer manufacturing or ellipsometry, the pBRDF is modeled by Fresnel
interactions at multiple boundaries, e.g. thin film stacks [7]. Some work has
even been done in efficiently parameterizing polarized scattering from rough
surfaces using Fresnel reflections, though depolarization is not accounted for
ignored[8]. However, in the regime of everyday materials found within common
indoor environments, the polarization attributes are dominated by
depolarization. In these scenes, the polarization attributes are relatively
small, but nonetheless informative and have demonstrated utility.

Interest in polarization for computer graphics and vision algorithms has grown
since commercial off-the-shelf (COTS) polarization cameras became available in
2018. The ability to track polarization in physics-based rendering (PBR)
engines increased in relevance as polarimetric measurements became more
accessible [9, 10, 11, 12, 13]. Very recently, data-driven models successfully
used COTS polarimeters for shape from polarization (SfP) and material texture
discrimination [14, 15, 16, 17, 18]. These studies confirm that polarimetry
offers unique information about a scene; however, analytic pBRDF models for
common indoor materials remain a barrier to exploring new applications using
physics-based approaches.

Both BRDF and pBRDF models typically contain a specular component that
describes light scattered from the surface of a material and a diffuse
component attributed to light scattered from within the material [19]. The
polarization of the specular component, more accurately called a first-surface
reflection component, is commonly modeled as Fresnel reflection from
hypothetical sub-resolution microfacets[20, 21]. Early work on polarized light
scattering assumed that the diffuse component becomes completely depolarized
such that diffuse contributions to reflection can be eliminated by using a
polarizer [22]. Although first-surface reflection tends to produce greater
polarimetric modulation, it is generally incorrect to assume that diffuse
reflection is completely unpolarized [23]. Many SfP approaches consider both
first-surface and diffuse polarized light scattering but use either a purely
specular or purely diffuse pBRDF for a given material[23, 24, 25]. Data driven
models have had success with a decomposition into diffuse and specular terms
as input into neural networks [13, 15].

More rigorous pBRDF models characterize polarimetric light scattering in terms
of a combination of first-surface and diffuse polarized components. The pBRDF
model introduced by Baek et al. consists of a first-surface microfacet term
and a polarized diffuse term[26]. The diffuse term describes Fresnel
transmission into, depolarizing scattering within, and Fresnel transmission
out of the material. In this model, depolarization arises from the diffuse
reflection term and from the opposing polarizance orientations between first-
surface microfacet and diffuse terms. Kondo et al. extended this model by
including an ideal depolarizer as an independent third term [27].

An empirical BRDF is generated by measuring the reflected radiance
distribution over the hemisphere of output directions for every point in the
hemisphere of illumination directions. A popular way to expedite this process
is to image a homogeneous spherical object with multiple illumination
geometries[28, 29]. Imaging a sphere captures many scattering geometries
simultaneously, and intermediate scattering geometries can be interpolated. A
pBRDF can be obtained in a similar way by performing MM imaging of homogeneous
spheres as in work by Baek et al.[30, 26].

Complete MM imaging requires at least sixteen polarimetric linearly-
independent measurements to constrain the sixteen degrees of freedom of a
MM[31, 6]. For measuring general MMs, it is an established practice to use
more than sixteen measurements because the pseudoinverse solution becomes more
robust to noise as additional linearly-dependent measurements are
performed[32, 33, 34, 35]. Dual rotating retarder (DRR) polarimeters are
commonly used to perform MM imaging and a DRR polarimeter is used in this work
as a baseline for polarimetric accuracy[36, 37, 38]. The design of MM imaging
measurement protocols are well-studied [32, 33, 34, 35]. DRR polarimeters (and
similar polarimeters which modulate the polarized signal with respect to time)
can be considered ill-suited for measuring objects with polarization that
changes in time and applications where a large number of objects or geometries
must be measured. In these cases, partial polarimetry can be employed with the
goal of maximizing polarimetric information in a reduced period of time[39].

Partial MM polarimeters perform fewer than sixteen linearly independent
measurements and therefore cannot reconstruct the full MM. Such polarimeters
and techniques to analyze their data have been developed for use in fields
such as biomedicine and computer vision. Qi et al. demonstrated the use of a
linear partial MM polarimeter for endoscopy[40]. For polarimetric imaging of
biological tissue near normal incidence, Novikova et al. showed the analysis
of partial MM data can yield comparable results to complete MM data[41].
Gonzalez et al. introduced a decomposition technique for MMs lacking the
fourth row[42]. Kupinski et al. optimized the polarization state generator
(PSG) and polarization state analyzer (PSA) states to perform a binary
classification task[43].

This work makes use of the triple-degenerate (TD) MM model introduced by Li
and Kupinski[44]. The TD-MM model decouples depolarization from the largest
contributing MJM, given the symbol $\widehat{\mathbf{m}}_{0}$, in a Cloude
decomposition[45, 46]. Brosseau showed a similar decomposition of a MM into a
MJM and an ideal depolarizer which performs an equivalent transformation of an
input Stokes vector to an output Stokes vector as a depolarizing MM, where the
MJM component was dependent on the input Stokes vector[47]. In a TD-MM model,
$\widehat{\mathbf{m}}_{0}$ does not depend on the incident Stokes vector.
Ossikovski and Arteaga defined the integral decomposition of a MM, where a
depolarizing MM is written as the sum of the largest contributing MJM and a
residual depolarizing MM[48]. The residual term is the sum of the three
remaining MJMs, using their corresponding coherency eigenvalues as weights.
The integral decomposition is similar to the TD model, but the TD model is the
special case where the three smaller eigenvalues are equal. Since these three
eigenvalues are equal, a fraction of the contribution of
$\widehat{\mathbf{m}}_{0}$ is combined with the remaining three MJMs to yield
$\mathbf{m}_{ID}$ as shown in Fig. 1, simplifying the description of
depolarization.

Prior work has demonstrated the broad applicability of this TD-MM model and
the use of the largest eigenvalue, $\xi_{0}$, as an important attribute that
depends on a material’s texture and albedo, and is a function of scattering
geometry. Omer and Kupinski demonstrated the efficiency of a TD representation
by truncating measured MM data to compress tabulated pBRDFs by 50%, with the
additional benefit of enforcing physicality when interpolating to unmeasured
geometries[49]. Jarecki and Kupinski demonstrated using an assumption of a TD
coherency eigenspectrum and modeling $\widehat{\mathbf{m}}_{0}$ as Fresnel
reflection to extrapolate MM images from linear Stokes measurements[50, 51].
They performed MM extrapolations on a set of two colors of plastic LEGO bricks
treated to have varying surface roughness. This dataset represented an
ensemble of many albedo and texture cases for a single material, and the
observed results obeyed Umov’s effect (_i.e._ that depolarization magnitude
should trend positively with albedo) and that depolarization should also
depend on texture and geometry[52, 53]. They used a microfacet-inspired
Fresnel reflection model for $\widehat{\mathbf{m}}_{0}$ and compared DRR
polarimeter measurements to extrapolations from a COTS linear Stokes camera.
However, the range of measured scattering geometries was small due to imaging
(macroscopically) flat samples, so purely first-surface reflection model based
on Fresnel reflection was sufficient to describe polarization scattering
behavior for that ensemble of materials. The material studied in this work is
better described using a model which combines Fresnel reflection and diffuse
polarization. The objective of this work is to introduce a TD-MM model that
efficiently describes measurements of spherical objects in a DRR Mueller
polarimeter. The accuracy of this model is assessed as performance in the task
of estimating $\xi_{0}$ from an assumed form of $\widehat{\mathbf{m}}_{0}$ and
partial polarimetric measurements.

The details of the model are described in Sect. III, the method for estimating
$\xi_{0}$ is described in Sect. IV, and quantitative results for a 3D printed
red sphere and qualitative results for a 3D printed Stanford bunny of the same
material are described in Sect. V

##  II A Triple Degenerate pBRDF for Depolarization-Dominant Materials

(a) Cloude decomposition

(b) Equal Weights

(c) TD Model

Figure 1: Visualization of Cloude decomposition where
$\xi_{1}=\xi_{2}=\xi_{3}=\overline{\xi}$ is denoted a triple degenerate (TD)
model. As shown in Eq. 1, the relative weight between
$\widehat{\mathbf{m}}_{0}$ and $\mathbf{m}_{ID}$ is controlled by the largest
coherency eigenvalue $\xi_{0}$, making it a single-parameter depolarization
model. In (a), the Cloude decomposition of a depolarizing MM into a convex sum
of four MJMs where the fractional size of each wedge corresponds to the
coherency eigenvalues and the weight of each MJM. In (b), the TD weights allow
the four MJMs to be separated in equal contribution. In (c), the equal
contributions combine to form an ideal depolarizer $\mathbf{m}_{ID}$.

A MM maps an input polarization state to an output polarization state[6, 31].
A MM has sixteen degrees of freedom: one for average throughput, three for
diattenuation, three for retardance, and nine for depolarization[6].
Diattenuation is the polarization-dependent reflectance which is encoded in
the top row of the MM and denoted by the vector$\mathbf{D}$. Retardance is the
polarization-dependent phase. Depolarization is the randomization of
polarization with respect to position, angle, time, and/or wavelength faster
than the particular detector in use can resolve[54]. Depolarization is of
particular interest in this work because it is the dominant polarization
property in most every-day environments. Several different metrics exist for
depolarization such as the depolarization index and the coherency
eigenvalues[55, 56].

The eigenvectors of a MM do not strictly correspond to physical Stokes
parameters. Instead, eigenanalysis is performed on the linearly related
coherency matrix[45, 46]. The eigenvectors of the coherency matrix can be used
to decompose a depolarizing MM into a convex sum of up to four non-
depolarizing MMs $\mathbf{m}=\sum_{n=0}^{3}\xi_{n}\widehat{\mathbf{m}}_{n}$.
Each of these non-depolarizing MM has an equivalent Jones matrix
representation $\mathbf{J}$ and therefore is also called a Mueller-Jones
matrix (MJM), and is notated with a hat $\widehat{\cdot}$. The $\xi_{n}$ in
the convex sum are called the coherency matrix eigenvalues, and are normalized
such that $\sum_{n=0}^{3}\xi_{n}=1$. With this normalization, each eigenvalue
$\xi_{n}$ represents the fraction of light which is transformed by the MJM
$\widehat{\mathbf{m}}_{n}$ as shown in Fig. 1 (a). The set of four MJM in this
decomposition are orthonormal in the sense that
$\frac{1}{2}\mathrm{tr}(\mathbf{J}_{n}^{\dagger}\mathbf{J}_{m})=\delta_{nm}$
and $\frac{1}{4}\sum_{n=0}^{3}\widehat{\mathbf{m}}_{n}=\mathbf{m}_{ID}$.

A useful special case is a triply-degenerate (TD) eigenspectrum
($\xi_{1}=\xi_{2}=\xi_{3}=\overline{\xi}$), where
$\overline{\xi}=(1-\xi_{0})/3$ as in Fig. 1 (b). In the TD case, the
depolarizing MM can be written as the convex sum of a MJM and an ideal
depolarizer

$\mathbf{M}=\frac{4M_{00}}{3}\left[\left(\xi_{0}-\frac{1}{4}\right)\widehat{\mathbf{m}}_{0}+\left(1-\xi_{0}\right)\mathbf{m}_{ID}\right],$
(1)

where $\xi_{0}$ is the largest eigenvalue, $\widehat{\mathbf{m}}_{0}$ is the
dominant MJM, $\mathbf{m}_{ID}$ which is zero everywhere except the top-left
element equals one, and $M_{00}$ is the throughput for unpolarized light. The
model in Eq. 1 is represented by the pie chart in Fig. 1 (c).

In a TD-MM model, the degrees of freedom are reduced from 16 to eight: one for
throughput, one for depolarization $\xi_{0}$, and six for the dominant MJM
$\widehat{\mathbf{m}}_{0}$. The single degree of freedom for depolarization is
$\xi_{0}$ which controls the relative weight between an ideal depolarizer and
a MJM. The largest eigenvalue is bounded in the range $0.25\leq\xi_{0}\leq
1.0$, where $\xi_{0}=0.25$ means the MM is the ideal depolarizer and
$\xi_{0}=1.0$ means the MM is a MJM. In a TD-MM model, the depolarization
index $DI$ is monotonically related to $\xi_{0}$ by the equation
$DI=\frac{1}{3}\sqrt{16\xi_{0}^{2}-8\xi_{0}+1}$. Diattenuation and retardance
orientations match those of $\widehat{\mathbf{m}}_{0}$ and are invariant to
$\xi_{0}$. The maximum diattenuation and retardance magnitudes in a TD-MM
match those of $\widehat{\mathbf{m}}_{0}$ when $\xi_{0}=1$, but are reduced as
$\xi_{0}$ approaches $0.25$ where $\mathbf{m}_{ID}$ dominates the sum.

(a) $\xi_{1}-\overline{\xi}$ 451 nm

(b) $\xi_{1}-\overline{\xi}$ 662 nm

(c) $\xi_{2}-\overline{\xi}$ 451 nm

(d) $\xi_{2}-\overline{\xi}$ 662 nm

(e)  $\xi_{3}-\overline{\xi}$ 451 nm

(f)  $\xi_{3}-\overline{\xi}$ 662 nm

Figure 2: Deviation of each of the smaller three coherency eigenvalues
$\xi_{1}$, $\xi_{2}$, and $\xi_{3}$ from $\overline{\xi}=(1-\xi_{0})/3$ for a
red 3D printed bunny. Values shown in (a), (c), (e) are at 451 nm and values
shown in (b), (d), (f) are at 662 nm. Each column shares a colorbar. When
$\xi_{1}$, $\xi_{2}$, and $\xi_{3}$ are all equal to $\overline{\xi}$, then
the coherency eigenspectrum is TD and the depolarizing MM can be exactly
described as a sum of a dominant MJM and an ideal depolarizer (see Fig. 1).
The eigenspectrum for the red 3D printed bunny under 451 nm illumination (low
albedo, low depolarization per Umov’s effect) has larger deviation from TD
than 662 nm (high albedo, high depolarization). This trend with wavelength
suggests that even if the dominant MJM model $\widehat{\mathbf{m}}_{0}$ is
accurate, the error in assuming the MM fits Eq. 1 would be greater at 451 nm
than 662 nm. In particular, regions of high curvature deviate more from TD as
shown on the left side of the bunny’s hind leg.

To demonstrate the appropriateness of a TD assumption for the red 3D printing
material used in this work, Fig. 2 shows the deviation of the coherency
eigenspectrum from TD of a printed Stanford bunny at two wavelengths. Since
the eigenvalues are normalized to sum to unity, the difference of a given
eigenvalue from $\overline{\xi}$ can be thought of as the fraction of light
that is not described by a TD model. The low albedo case of 451 nm
illumination deviates from TD for as much as 9% of the light while the high
albedo case of 662 nm illumination only deviates as much as 3%. For both
wavelengths, the deviation from TD is larger in regions of high curvature.
This could be due to increased likelihood of multiple ray bounces resulting in
more complex polarized light scattering behavior. In regions with high
deviation from a TD eigenspectrum, tasks such as estimation of $\xi_{0}$ will
be less accurate even when $\widehat{\mathbf{m}}_{0}$ is perfectly known.

##  III Mueller-Jones Model of First-Surface and Diffuse Polarization
Attributes

###  III-A Rusinkiewicz Coordinate System

The measurement configuration for a spherical object and the parameters used
to define the scattering geometry are shown in Figure 3 (a). The green light
bulb indicates the source which is connected to a point on the object by the
vector $\widehat{\boldsymbol{\omega}}_{i}$. The red vector from this object
point to a given camera element is $\widehat{\boldsymbol{\omega}}_{o}$. The
angle between the source and the camera is denoted $\Omega$. Figure 3 (b)
shows a zoomed in view where these two vectors and the object surface normal
are parameterized in a Rusinkiewicz coordinate system.

(a) Measurement configuration

(b) Rusinkiewicz coordinate system

Figure 3: In (a), a spherical object by a camera and light source separated by
an angle $\Omega$. The point on the sphere imaged by a given pixel is observed
along the red vector $\widehat{\boldsymbol{\omega}}_{o}$ and is illuminated
along the green vector $\widehat{\boldsymbol{\omega}}_{i}$. This point on the
sphere has a surface normal along the black vector $\widehat{\mathbf{n}}$. The
blue halfway vector $\widehat{\mathbf{h}}$ bisects
$\widehat{\boldsymbol{\omega}}_{o}$ and $\widehat{\boldsymbol{\omega}}_{i}$
and is used in (b) to parameterize the scattering geometry in terms of
Rusinkiewicz coordinates. This work considers isotropic surfaces where the
pBRDF is invariant to $\phi_{h}$ since $\widehat{\mathbf{t}}$ is arbitrary.

Rusinkiewicz coordinates are an established parameterization that reduces
dimensionality for isotropic surfaces and recasts the scattering geometry in
terms of variables with useful physical interpretations[57]. The
parameterization defines a halfway vector $\widehat{\mathbf{h}}$ (so named for
bisecting $\widehat{\boldsymbol{\omega}}_{i}$ and
$\widehat{\boldsymbol{\omega}}_{o}$) in terms of spherical coordinates
$\theta_{h}$ and $\phi_{h}$ about the surface normal $\widehat{\mathbf{n}}$.
The incoming and outgoing ray directions, $\widehat{\boldsymbol{\omega}}_{i}$
and $\widehat{\boldsymbol{\omega}}_{o}$, are then defined with spherical
coordinates $\theta_{d}$ and $\phi_{d}$ about $\widehat{\mathbf{h}}$. The
angle between the surface normal and the halfway vector, $\theta_{h}$,
describes how a particular scattering geometry deviates from specular
reflection. $\phi_{h}$ orients the scattering geometry relative to an axis of
anisotropy $\widehat{\mathbf{t}}$. For isotropic surfaces the BRDF/pBRDF is
invariant to $\phi_{h}$, reducing the functional dependence on four variables
to three. $\theta_{d}$ can be thought of as the angle of incidence onto a
hypothetical microfacet as further described below. Last, $\phi_{d}$ describes
the deviation of a scattering geometry from being “in-plane.” “Out-of-plane”
scattering produces changes in polarization due to geometric effects rather
than polarization-dependent reflectance or optical path length.

(a) Sphere $\theta_{h}$

(b) Bunny $\theta_{h}$

(c) Sphere $\theta_{d}$

(d) Bunny $\theta_{d}$

(e) Sphere $\phi_{d}$

(f) Bunny $\phi_{d}$

Figure 4: Rusinkiewicz coordinates $\theta_{h}$, $\theta_{d}$, and $\phi_{d}$
are shown in (a), (c), and (e), respectively, for a spherical object where
$\Omega=35^{\circ}$ and in (b), (d), and (f), respectively, a Stanford bunny.
The Rusinkiewicz definitions are shown in Fig. 3. The regions in (a) and (b)
where $\theta_{h}$ is small correspond to near-specular geometries where
Fresnel reflection dominates.

###  III-B First-Surface Reflection Component

The first-surface reflection, frequently referred to as specular reflection,
is commonly modeled as Fresnel reflection from a hypothetical sub-resolution
feature called a microfacet (the term first-surface reflection will be used
going forward so that specular may be reserved to describe scattering
configurations where $\theta_{h}=0^{\circ}$)[21]. The microfacet is oriented
such that it satisfies the Law of Reflection for a given pair of input and
output ray directions.

(a) First surface reflection MJM

(b) Diffuse polarization MJM

Figure 5: MM images of (a) the first-surface reflection model in Eq. 2 and (b)
the diffuse model in Eq. 3 evaluated over a sphere in the same scattering
configuration as in Fig. 4. The first-surface reflection component is modeled
as Fresnel reflection from microfacets, so the polarized scattering behavior
does not depend directly on the surface normal. For this reason, the
diattenuation and polarizance orientations do not follow the surface of the
sphere. The non-zero values present in the $m_{1,2}$ and $m_{2,1}$ elements of
(a) and (b) represent changes in the polarization state due to geometric
effects rather than retardance.

This model is conceptually simple, however, implementation is complicated by
the need to consider rotations in the plane transverse to the direction of
propagation. The orientations of $s$\- and $p$-polarizations in space vary
based on the scattering geometry, and the orientations of horizontal and
vertical for both the input and output coordinate systems are implicit
parameters of the Jones matrix (JM) or MM. The JM for Fresnel reflection from
a microfacet $\mathbf{F}$ is

$\\!\\!\\!\\!\\!\\!\mathbf{F}_{n_{\lambda}}\\!(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o})\\!=\\!\\!\begin{bmatrix}\widehat{\mathbf{x}}_{o}\\\
\widehat{\mathbf{y}}_{o}\end{bmatrix}\\!\\!\\!\begin{bmatrix}\widehat{\mathbf{s}}_{o}\\\
\widehat{\mathbf{p}}_{o}\end{bmatrix}^{\mathrm{T}}\\!\\!\\!\begin{bmatrix}r_{s}(n_{\lambda},\theta_{d})&0\\\
0&r_{p}(n_{\lambda},\theta_{d})\end{bmatrix}\\!\\!\\!\begin{bmatrix}\widehat{\mathbf{s}}_{i}\\\
\widehat{\mathbf{p}}_{i}\end{bmatrix}\\!\\!\\!\begin{bmatrix}\widehat{\mathbf{x}}_{i}\\\
\widehat{\mathbf{y}}_{i}\end{bmatrix}^{\mathrm{T}}\\!\\!\\!\\!\\!.$ (2)

Here $\widehat{\boldsymbol{\omega}}_{i}$ and
$\widehat{\boldsymbol{\omega}}_{o}$ are the input (i) and output (o) ray
directions, respectively; $\widehat{\mathbf{x}}_{i/o}$,
$\widehat{\mathbf{y}}_{i/o}$ are the horizontal and vertical directions for
the input and output coordinate system, respectively;
$\widehat{\mathbf{s}}_{i/o}$ and $\widehat{\mathbf{p}}_{i/o}$ are the
directions of $s$\- and $p$-polarization for the incoming and outgoing rays
with respect to the microfacet, respectively. The directions of $x$\- and
$y$-polarization are found
$\widehat{\mathbf{x}}_{i/o}=\widehat{\mathbf{n}}\times\widehat{\boldsymbol{\omega}}_{i/o}$
and
$\widehat{\mathbf{y}}_{i/o}=\widehat{\boldsymbol{\omega}}_{i/o}\times\widehat{\mathbf{x}}_{i/o}$.
The directions of $s$\- and $p$-polarization are found
$\widehat{\mathbf{s}}_{i}=\widehat{\mathbf{s}}_{o}=\widehat{\mathbf{h}}\times\widehat{\boldsymbol{\omega}}_{i}$
and
$\widehat{\mathbf{p}}_{i/o}=\widehat{\boldsymbol{\omega}}_{i/o}\times\widehat{\mathbf{s}}_{i/o}$.
The Fresnel coefficients $r_{s}$ and $r_{p}$ depend on the refractive index
ratio $n_{\lambda}$ and the angle of incidence onto the microfacet
$\theta_{d}$[21]. The values used for $n_{\lambda}$ are presented in Tab. I.

Eq. 2 is invariant to the surface normal since Fresnel reflection is only
defined for a normal vector that is halfway between the illumination and view
vectors (i.e. a microfacet). Instead, the Fresnel term only depends on the
incident and scattered ray directions $\widehat{\boldsymbol{\omega}}_{i}$ and
$\widehat{\boldsymbol{\omega}}_{o}$. The diattenuation magnitude depends only
on $\theta_{d}$, which is half of the angle between
$\widehat{\boldsymbol{\omega}}_{i}$ and $\widehat{\boldsymbol{\omega}}_{o}$.
$\theta_{d}$ varies slowly over the sphere as shown in the example geometry in
Fig. 4, so the diattenuation magnitude also varies slowly. Variation in the
diattenuation orientation is caused by variation in the orientation of $s$\-
and $p$-polarizations relative to the input and output coordinate systems.
Both the $s$\- and $p$-directions and the input and output coordinate systems
vary over the field of view. These coordinate transforms lead to coupling
between linear polarization states which can be observed faintly in the center
$2\times 2$ elements of the MM image in Fig. 6 (a). This should not be
confused with true circular retardance, as it is the result of geometric
effects rather than a relative phase difference between polarization states.
When the microfacet surface normal $\widehat{\mathbf{h}}$ and the macroscopic
surface normal $\widehat{\mathbf{n}}$ are coplanar (and therefore also
coplanar with $\widehat{\boldsymbol{\omega}}_{i}$ and
$\widehat{\boldsymbol{\omega}}_{o}$), there is no coupling between linear
states caused by geometric effects. This can be noted in the MM image of a
purely Fresnel reflection $\widehat{\mathbf{m}}_{0}$ is shown in Fig. 5 (a),
where there is a dark band through the center of the $m_{1,2}$ and $m_{2,1}$
elements. The geometries in this region are referred to as “in-plane”
scattering.

It is more correct to think of
$\mathbf{J}_{s}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o})$
as a first-surface reflection model rather than a microfacet model. This is
because the distribution of microfacet normals and adjacency effects are not
explicitly included. Conventional microfacet models use distributions on the
orientations of the microfacets to model surface roughness [58, 59, 60, 61].
Adjacency effects, specifically shadowing and masking, and the microfacet
distribution can be interpreted as dictating the likelihood of an incident ray
undergoing first-surface reflection toward the detector[62]. This likelihood
corresponds to the relative contribution of first-surface reflection to other
terms in the pBRDF model. In this work, the adjacency effects and microfacet
distribution are absorbed into the dominant eigenvalue
$\xi_{0}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o})$
from Eq. 1 which controls the weight between and the first-surface to diffuse
weighting function $f_{\lambda}(\theta_{h})$.

###  III-C Diffuse Reflection Component

The diffuse reflection component is modeled as

$\displaystyle\mathbf{S}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o},\widehat{\mathbf{n}})$
$\displaystyle=\begin{bmatrix}1&0\\\
0&-1\end{bmatrix}\mathbf{R}(\phi_{d})\begin{bmatrix}1&0\\\
0&0\end{bmatrix}\mathbf{R}(-\phi_{d}),$ (3)

where the rightmost three JMs represent a polarizer with unit diattenuation
and a transmission axis oriented at $\phi_{d}$. This matrix has a
diattenuation orientation that is in the plane spanned by the surface normal
$\widehat{\mathbf{n}}$ and $\widehat{\boldsymbol{\omega}}_{o}$, a pattern
which is centered about the point on the sphere where $\theta_{h}=0^{\circ}$
as shown in Fig. 4 (a). When left multiplied by the reflection matrix, the
$45^{\circ}$ and $135^{\circ}$ polarizations are reversed. This pattern being
centered about $\theta_{h}=0^{\circ}$ is consistent with the
$\widehat{\mathbf{m}}_{0}$ images observed from DRR measurements of spheres
(see Fig. 6 (c)), however, it does mark a departure from many other diffuse
polarization models. In models such as the one proposed by Atkinson and
Hancock, the orientation of the diffuse term is centered about the central
camera axis[23].

###  III-D Mixed Dominant Mueller-Jones Matrix Model

Previous work demonstrated extrapolating MM images using a purely Fresnel
reflection model for the dominant MJM, but this work was performed on flat
objects which have a smaller range of scattering geometries[50, 51]. The
objects and scattering geometries described in this work require models for
$\widehat{\mathbf{m}}_{0}$ that introduce a polarized diffuse reflection term.
This term is characterized by a diattenuation orientation that depends on the
surface normal. However, the dominant MJM necessarily must be non-
depolarizing, and simply summing different MMs, in general, introduces
depolarization.

To create a mixed polarization model without introducing depolarization, the
individual first-surface and diffuse components are combined as Jones matrices
as in

$\mathbf{J}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o},\widehat{\mathbf{n}}|n_{\lambda},a_{\lambda},b_{\lambda})=\mathbf{F}_{n_{\lambda}}+a_{\lambda}\mathrm{sin}^{b_{\lambda}}(\theta_{h})\mathbf{S}(\widehat{\mathbf{n}}).$
(4)

Here, the functional dependencies of $\mathbf{J}$ are separated into the
scattering geometry and material parameters, with the separation notated by
the vertical bar $|$. There are four material-dependent parameters which are
constant for all scattering geometries: the real and imaginary components of
the refractive index $n_{\lambda}$, and two parameters $a_{\lambda}$ and
$b_{\lambda}$ which control how the weight of the diffuse polarization term
depends on scattering geometry. The four material parameters efficiently
reduce the original six degrees of freedom for the diattenuation and
retardance of the MJM. Both matrices on the right-hand-side depend on
$\widehat{\boldsymbol{\omega}}_{i}$ and $\widehat{\boldsymbol{\omega}}_{o}$,
but this has been dropped for brevity.

The surface normal of the material $\widehat{\mathbf{n}}$ (see Fig. 3) affects
the diffuse but not Fresnel term.
$\mathbf{F}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o},n_{\lambda})$
is the first-surface reflection component modeled as Fresnel reflection based
on refractive index ratio $n_{\lambda}$ at wavelength $\lambda$,
$\mathbf{S}(\widehat{\boldsymbol{\omega}}_{i},\widehat{\boldsymbol{\omega}}_{o},\widehat{\mathbf{n}})$
is the diffuse polarization term and is defined with unit throughput, so its
relative weight as a function of scattering geometry is given by the sine
function. $\mathbf{J}$ is converted to a MJM, normalized, and then used as
$\widehat{\mathbf{m}}_{0}$ for a TD-MM model as in Eq. 1.

TABLE I: Material constants determined _ad hoc_ to match the models to the observed MM measurements. These parameters are constant with respect to scattering geometry. Material constants | 451 nm | 662 nm
---|---|---
$n_{\lambda}$ | $1.20+0.25i$ | $1.30+0.08i$
$a_{\lambda}$ | 0.03 | 0.17
$b_{\lambda}$ | 2.5 | 2.0

The values for the material properties are found in Tab. I. These values, as
well as the choice to use a sine function for the geometric dependence of the
weight, were determined _ad hoc_ to match the measurements. The weight
function is modeled as only depending on $\theta_{h}$ because $\theta_{h}$ is
a measure of deviation from a specular scattering configuration. The sine
function was chosen for the weight functions based on the observation that
Fresnel reflection-like behavior dominates for near-specular scattering
geometries (where $\theta_{h}$ is small) and diffuse polarization behavior
tends to become dominant away from specular geometries. The coefficients and
powers in these equations were empirically chosen to match the diattenuation
trends observed in measurements.

(a) Measured $\widehat{\mathbf{m}}_{0}$ 451 nm

(b) Modeled $\widehat{\mathbf{m}}_{0}$ 451 nm

(c) Measured $\widehat{\mathbf{m}}_{0}$ 662 nm

(d) Modeled $\widehat{\mathbf{m}}_{0}$ 662 nm

(e) Measured $\widehat{\mathbf{m}}_{0}$ 451 nm

(f) Modeled $\widehat{\mathbf{m}}_{0}$ 451 nm

(g) Measured $\widehat{\mathbf{m}}_{0}$ 662 nm

(h) Modeled $\widehat{\mathbf{m}}_{0}$ 662 nm

Figure 6: Comparisons of dominant MJM from DRR measurements at
$\Omega=35^{\circ}$ and from the models (see Eq. 4) for a red sphere (top row)
and a red bunny (bottom row) of the same material. At 451 nm (left side) the
low albedo creates lower diffuse scattering compared to 662 nm (right side).
Properties of the capture system are not included in the MJM images from the
models, specifically the finite spatial and polarimetric resolutions.
Additionally, the model bunny MJM images are based on the geometry captured
for the MM extrapolation from Stokes data experiment detailed in Sect. V C.
There is therefore some inherent visual disagreement due to the geometry.

##  IV Eigenvalue Estimation Method

A collection of noise-free polarimetric measurements $\mathbf{P}$ of a general
MM is modeled as $\mathbf{P}=\mathbf{W}_{N}\mathbf{M}$, where $\mathbf{W}_{N}$
is the polarimetric measurement matrix with $N$ number of measurements and the
MM is expressed as a 16 $\times$ 1 vector $\mathbf{M}$. If $\mathbf{M}$ is
accurately described by a TD-MM model as in Eq. 1, then the noise-free
polarimetric measurements can be modeled as
$\mathbf{P}=\bm{\Phi}\boldsymbol{\alpha}.$ Here $\bm{\Phi}$ is an $N\times 2$
matrix

$\bm{\Phi}=\begin{bmatrix}\mathbf{p}_{0}&\mathbf{p}_{ID}\end{bmatrix}=\mathbf{W}_{N}\begin{bmatrix}\widehat{\mathbf{m}}_{0}&\mathbf{m}_{ID}\end{bmatrix},$
(5)

where $\boldsymbol{\alpha}$ is a vector whose elements are the coefficients on
each MM in the TD model

$\boldsymbol{\alpha}=\begin{bmatrix}\alpha_{0}\\\
\alpha_{ID}\end{bmatrix}=\frac{4M_{00}}{3}\begin{bmatrix}\xi_{0}-\frac{1}{4}\\\
1-\xi_{0}\end{bmatrix}.$ (6)

The estimate for $\widetilde{\boldsymbol{\alpha}}$ is given by applying the
pseudoinverse of $\bm{\Phi}$ to noisy measurements $\widetilde{\mathbf{P}}$ as
in $\widetilde{\boldsymbol{\alpha}}=\bm{\Phi}^{+}\widetilde{\mathbf{P}}$. The
estimates for the dominant eigenvalue $\widetilde{\xi}_{0}$ and for the
average reflectance $\widetilde{M}_{00}$ are calculated from the elements of
$\widetilde{\boldsymbol{\alpha}}$

$\widetilde{\xi}_{0}=\frac{\frac{1}{4}+\widetilde{\alpha}_{0}/\widetilde{\alpha}_{ID}}{1+\widetilde{\alpha}_{0}/\widetilde{\alpha}_{ID}},$
(7)

and

$\widetilde{M}_{00}=\frac{3\widetilde{\alpha}_{0}}{4\widetilde{\xi}_{0}-1}.$
(8)

In this work, the estimated average reflectance $\widetilde{M}_{00}$ is not
considered because the polarimeters used are not calibrated to produce results
that are scaled as fractions of the illumination radiance. By normalizing the
results, or setting $M_{00}=1$, the polarimetric properties can be directly
compared.

For a general MM, $\mathbf{W}$ must be rank sixteen to perform a complete
reconstruction. In the case of a TD-MM when $\widehat{\mathbf{m}}_{0}$ is
known _a priori_ or provided by a model, the two unknown degrees of freedom
can be estimated with a rank two measurement matrix. The complete DRR Mueller
imaging polarimeter used in this work performs 40 polarimetric measurements,
$\mathbf{W}_{40}$, making it an overdetermined system [38]. The partial
polarimeter used in this work is a Sony Triton 5.0MP Polarization Camera. This
is a division of the focal plane linear Stokes camera. The light source from
the DRR polarimeter is used to illuminate the object with horizontally
polarized light. The measurement matrix for the linear Stokes camera,
$\mathbf{W}_{4}$, has four rows but is rank three. This is underdetermined for
reconstructing a general MM, but is overdetermined for estimating the two
remaining degrees of freedom in Eq. 6.

##  V Results

###  V-A Dominant MJM Model Comparison

The measured and modeled $\widehat{\mathbf{m}}_{0}$ are compared by their
diattenuation magnitudes $|\mathbf{D}|$ and orientations $\psi$. Diattenuation
is represented in the top row of a MM. In Fig. 7 $\psi$ and in Fig. 8
$|\mathbf{D}|$ are compared at $\Omega=35^{\circ}$ based on the Rusinkiewicz
coordinates from the DRR measurements rather than with the LUT method used in
Sect. V-B because the models can be evaluated at arbitrary scattering
geometries. In Fig. 7 (a) and (c), $\psi$ is primarily vertical due to the low
albedo case of 451 nm illumination being dominated by Frensel reflection. In
Fig. 7 (b) and (d), $\psi$ is vertical in the center where Fresnel reflection
dominates the near-specular geometries but the diffuse term dominates away
from specular, which is consistent with high albedo behavior.

(a) Measured $\psi$ 451 nm

(b) Measured $\psi$ 662 nm

(c) Modeled $\widetilde{\psi}$ 451 nm

(d) Modeled $\widetilde{\psi}$ 662 nm

(e) $\psi$ deviation 451 nm

(f) $\psi$ deviation 662 nm

Figure 7: Diattenuation orientation $\psi$ images of the
$\widehat{\mathbf{m}}_{0}$ models at $\Omega=35^{\circ}$ used in the
eigenvalue estimation at (a) 451 nm and at (b) 662 nm, and of the
$\widehat{\mathbf{m}}_{0}$ calculated from DRR MM polarimeter data at (c) 451
nm and at (d) 662 nm. The angle between the measured and modeled diattenuation
orientations is shown in (e) and (f) for 451 nm and 662 nm respectively. The
diattenuation orientation in (a) is primarily vertically oriented which is
consistent with Fresnel reflection-dominated scattering as expected from a low
albedo case. The orientation in (b) processes smoothly around the edge of the
sphere which is consistent with diffuse polarization. This procession is
radially oriented but with the 45∘/135∘ components flipped due to the
measurement being a reflection configuration. The center of the pattern shows
vertically oriented diattenuation which suggests that Fresnel reflection
dominates in that region. The $\widetilde{\psi}$ values in the extrapolations
(c) and (d) are invariant to the estimation of $\widetilde{\xi}_{0}$ and only
depend on the model $\widehat{\mathbf{m}}_{0}$. Agreement between the
measurements and extrapolations is summarized according to Eq. 9 in Fig. 9
(a).

(a) Model $|\mathbf{D}|$ 451 nm

(b) Model $|\widetilde{\mathbf{D}}|$ 662 nm

(c) True $|\mathbf{D}|$ 451 nm

(d) True $|\widetilde{\mathbf{D}}|$ 662 nm

(e) $|\mathbf{D}|$ deviation 451 nm

(f) $|\widetilde{\mathbf{D}}|$ deviation 662 nm

Figure 8: Diattenuation magnitude $|\mathbf{D}|$ images for the red sphere at
$\Omega=35^{\circ}$. The ground truth from the DRR polarimeter measurements in
(a) and the extrapolation from linear Stokes data in (b) at 451 nm are the low
albedo case. The ground truth from the DRR measurements in (c) and the
extrapolation from linear Stokes data in (d) at 662 nm are the high albedo
case. The agreement is summarized with respect to $\Omega$ in Fig. 9 (b). The
two minima in (c) and (d), referred to as “neutral points,” indicate the
presence of at least two polarimetric processes which have comparable
diattenuation magnitude and opposing orientations.

The agreement in $\psi$ and $|\mathbf{D}|$ are summarized numerically in Fig.
9 according to their root-mean-squared deviation (RMSD) values over the sphere
at each acquisition geometry $\Omega$ as in

$\begin{split}\epsilon(\widetilde{\psi},\psi)&=\sqrt{\frac{1}{K}\sum_{n=1}^{K}||(\widetilde{\psi}_{k}-\psi_{k})||^{2}}\\\
&=\sqrt{\frac{1}{K}\sum^{K}_{n=1}\left[\frac{1}{2}\text{acos}\left(\widetilde{\mathbf{d}}_{k}\cdot\mathbf{d}_{k}\right)\right]^{2}},\end{split}$
(9)

and

$\epsilon(|\widetilde{\mathbf{D}}|,|\mathbf{D}|)=\sqrt{\frac{1}{K}\sum_{n=1}^{K}\left|\left||\widetilde{\mathbf{D}}|_{k}-|\mathbf{D}|_{k}\right|\right|^{2}}$
(10)

where $K$ is the total number pixels, and $k$ is the pixel index. The dot
product of the normalized linear diattenuation parameters,
$\mathbf{d}=[D_{1},D_{2}]/|\mathbf{D}|$, is used to avoid angular phase
wrapping effects caused by the fact that $\psi=0^{\circ}$ and
$\psi=180^{\circ}$ describe the same orientation.

In Fig. 9 (a), the error in diattenuation orientation trends downward with
$\Omega$. The MJM term is dominated by first-surface reflection which has a
simpler, easier to model orientation (primarily vertical). This is especially
true for 451 nm, where first-surface reflection is overall more dominant due
to being low albedo. The error for 451 nm hovers around $2^{\circ}$ for large
$\Omega$. The more complicated orientation pattern for diffuse reflection is
matched less well as shown in the generally larger error for 662 nm. In Fig. 9
(b), the error in diattenuation magnitude is strictly lower for 451 nm than
for 662 nm for all $\Omega$. This is likely due to the simpler first-surface
model dominating at 451 nm whereas the more complex diffuse model has a
stronger contribution at 662 nm.

(a) Diattenuation orientation error

(b) Diattenuation magnitude error

Figure 9: Error (a) in the diattenuation orientation $\widetilde{\psi}$
according to Eq. 9 and (b) in the diattenuation magnitude
$|\widetilde{\mathbf{D}}|$ according to Eq. 10 of the
$\widehat{\mathbf{m}}_{0}$ models at 451 nm (blue squares) and 662 nm (red
triangles). RMSD for fifteen acquisition geometries $\Omega$. The orientation
error for 451 nm hovers around $2^{\circ}$ for large $\Omega$. The more
complicated orientation pattern for diffuse reflection, which is more
prominent at 662 nm, is matched less well as shown in the generally larger
error for this high albedo case. The error in diattenuation magnitude is
strictly lower for 451 nm than for 662 nm for all $\Omega$. Again, this is
likely due to the simpler first-surface model dominating at 451 nm whereas the
more complex diffuse model has a stronger contribution at 662 nm.

###  V-B $\widetilde{\xi}_{0}$ Estimation from Stokes for a Sphere

To measure polarized behavior at many scattering geometries in fewer
acquisitions, the measured object was a red 3D printed sphere of 1-inch
diameter. The sphere was positioned at the center of rotation of a goniometric
arm which was used to perform polarimetric measurements at fifteen angles
between the camera and the light source denoted with $\Omega$. Polarimetric
measurements were performed at 662 nm and 451 nm. These two wavelengths
represent high and low albedo cases, respectively, for the red sphere. Per
Umov’s effect, these correspond to an expectation of high depolarization (low
$\xi_{0}$) and low depolarization (high $\xi_{0}$), respectively.

The DRR MM polarimeter and linear Stokes camera were both used with 8 mm focal
length lenses but have different detector sizes and resolutions. Additionally,
the reliably positioning the two cameras to have the same view proved
challenging. A pixel-to-pixel comparison between DRR and Stokes camera images
would therefore mean comparing polarimetric behavior at different scattering
geometries. Instead, the estimated $\widetilde{\xi}_{0}$ from the linear
Stokes camera is compared to the values from the DRR polarimeter by tabulating
a lookup table (LUT) of $\xi_{0}$ values with respect to scattering geometry
for each polarimeter. These LUTs are binned according to three angles in the
Rusinkiewicz coordinate system, $\theta_{h}$, $\theta_{d}$, and $\phi_{d}$
with 91, 91, and 361 bins assigned respectively. The Rusinkiewicz coordinate
system is explained further in Sect. III-A. $\theta_{d}$ and $\theta_{h}$ are
binned in steps of 1∘ from $0^{\circ}$ to $90^{\circ}$, and $\phi_{d}$ is
binned in steps of 1∘ from $-180^{\circ}$ to $180^{\circ}$. The LUT is then
evaluated at the scattering geometries of a virtual sphere matching the linear
Stokes camera measurement to compare the true and estimated $\xi_{0}$ at
identical points.

(a) Estimated $\widetilde{\xi}_{0}$ 451 nm

(b) Estimated $\widetilde{\xi}_{0}$ 662 nm

(c) True $\xi_{0}$ 451 nm

(d) True $\xi_{0}$ 662 nm

(e)  $\xi_{0}$ deviation 451 nm

(f)  $\xi_{0}$ deviation 662 nm

Figure 10: Results for the estimation of $\widetilde{\xi}_{0}$ for angle
between the camera and light source $\Omega=35^{\circ}$. The estimated values
from the linear Stokes camera in (a) and the ground truth values from the DRR
polarimeter in (b) at 451 nm are compared in (e) by taking their difference.
The estimated values from the linear Stokes camera in (c) and the ground truth
values from the DRR polarimeter in (d) at 662 nm are compared in (f) by taking
their difference. 662 nm is high albedo for the red sphere, so depolarization
is high and $\xi_{0}$ is low. 451 nm is low albedo, so $\xi_{0}$ is larger. In
(e) and (f), positive values indicate an overestimation of the contribution of
$\widehat{\mathbf{m}}_{0}$ while negative values indicate an overestimation of
the contribution of $\mathbf{m}_{ID}$. Because the eigenvalues are normalized
to sum to unity, the errors shown in (e) and (f) can be thought of as related
to the fraction of the light that is attributed to the wrong MM component.

The estimations of $\widetilde{\xi}_{0}$ from the linear Stokes camera over
the sphere at $\Omega=35^{\circ}$ are shown in Fig. 10 (a) and (b), the ground
truth values measured in the DRR polarimeter are shown in Fig. 10 (c) and (d),
and the difference between estimate and ground truth are shown in Fig. 10 (e)
and (f). The white region in the center consists of specific Rusinkiewicz
coordinates that the instruments do not have in common. 662 nm is high albedo
for the red sphere, so the expectation according to Umov’s effect is that
depolarization is high and $\xi_{0}$ is low. Conversely, 451 nm is low albedo
so $\xi_{0}$ is expected to be larger. Both of these trends are observed in
the ground truth measurements as well as in the estimations. In (e) and (f),
positive values indicate an overestimation of the contribution of
$\widehat{\mathbf{m}}_{0}$ while negative values indicate an overestimation of
the contribution of $\mathbf{m}_{ID}$. Because the eigenvalues are normalized
to sum to unity, the errors shown in (e) and (f) can be thought of as related
to the fraction of the light that is attributed to the wrong MM component.

Figure 11: RMSD in $\widetilde{\xi_{0}}$ estimate from Stokes imaging as
compared to ground truth from DRR MM imaging according to Eq. 11 at 451 nm
(blue squares) and 662 nm (red triangles) at the fifteen acquisition
geometries $\Omega$. RMSD values at 451 nm, the low albedo case, are strictly
greater than those for 662 nm. This is potentially explained by the material
being less depolarizing according to Umov’s effect than at 662 nm and
therefore more sensitive to disagreement between the ground truth and modeled
$\widehat{\mathbf{m}}_{0}$. Additionally, the eigenspectrum at 451 nm is
further from TD than for 662 nm, so a TD model may be less appropriate. The
lower throughput due to the low albedo could exacerbate the effects of noise
and degrade the estimation due to the smaller dynamic range of the linear
Stokes camera.

The error averaged over the image of the sphere is summarized quantitatively
by the RMSD

$\epsilon(\widetilde{\xi}_{0},\xi_{0})=\sqrt{\frac{1}{K}\sum_{k=1}^{K}||\widetilde{\xi}_{0,k}-\xi_{0,k}||^{2}},$
(11)

where $K$ is the total number pixels, and $k$ is the pixel index. Figure 11
shows the trend of $\epsilon(\widetilde{\xi}_{0},\xi_{0})$ as a function of
the angle between the camera and light source $\Omega$. The error at 451 nm is
larger than at 662 nm for every $\Omega$. According to Umov’s effect, the
depolarization at 451 nm for a red object should be lower, suggesting that the
estimation is more sensitive to disagreement between the ground truth and
model $\widehat{\mathbf{m}}_{0}$. Additionally, the relatively lower overall
reflectance at 451 nm also makes the linear Stokes camera measurements more
susceptible to noise.

###  V-C Mueller Extrapolation from Stokes for a Stanford Bunny

To extend the demonstration of the efficient TD-MM model to non-spherical
objects, the MM image of a Stanford bunny was extrapolated from linear Stokes
camera measurements. This Stanford bunny was printed from the same material as
the sphere in Sect. V-B, so the same pBRDF models can be used. These models
are detailed in Sect. III.

Linear Stokes measurements are performed using the COTS camera and one
polarized illumination state from the DRR light source as in Sect. V-B. The
$\widehat{\mathbf{m}}_{0}$ models are evaluated at the Rusinkiewicz angles
calculated from the .stl file used to 3D print the Stanford bunny. The
$\widehat{\mathbf{m}}_{0}$ model and Stokes data are used to estimate
$\widetilde{\xi}_{0}$ following Sect. IV. MM extrapolation is performed by
plugging the estimated $\widetilde{\xi}_{0}$ back into Eq. 1.

Figure 12 shows a visual comparison of the Stanford bunny between various
polarizers using the DRR measurements (left) and extrapolated MM image from
linear Stokes data (right). Each pair of images is on the same color scale.
Crossed linear, aligned linear, and crossed circular measurements are
simulated. The MM images are normalized, so variation in the throughput of the
second polarizer is due to the diattenuation, depolarization, and geometric
transformation of the incident polarizer. The overall appearance between the
results from DRR measurement and Stokes extrapolation are in agreement,
particularly with respect to wavelength. At 451 nm, shown in Fig. 12 (a-c),
the smaller amount of depolarization expected for a low albedo material
appears as greater polarimetric modulation as compared to 662 nm, shown in
Fig. 12 (d-f). The overall trend with geometry tends to match as well: regions
of the bunny which appear brighter or darker are in the results from DRR
measurement correspond to brighter or darker regions in the results
extrapolation from Stokes data. Notably, the region with the most significant
difference is the left side of the bunny’s leg. This is the region where the
TD assumption is weakest due to a higher angle of scattering, as shown in Fig.
2. Additionally, this region of high curvature is where misregistration
between the modeled geometry shown in Fig. 4 (b,d,f) and the geometry measured
with the Stokes camera would cause large deviation in the modeled and measured
$\widehat{\mathbf{m}}_{0}$, reducing the quality of the extrapolation. There
is unavoidable disagreement between the DRR simulations and the linear Stokes
images due to even slight differences in viewing geometry, resolution, and
camera properties of the two polarimeters. An example of such disagreement due
to different camera systems can clearly be seen on the left ear which has a
significantly different shadow between DRR and Stokes results.

(a) Crossed linear polarizers 451 nm

(b) Aligned linear polarizers 451 nm

(c) Crossed circular polarizers 451 nm

(d) Crossed linear polarizers 662 nm

(e) Aligned linear polarizers 662 nm

(f) Crossed circular polarizers 662 nm

Figure 12: Each image pair in (a-f) is arranged so that the simulated
measurements are on the left the normalized MM image of the Stanford bunny
from the DRR polarimeter and (right) the normalized extrapolated MM image from
linear Stokes data at 451 nm. In (a) vertically polarized illumination and a
horizontal analyzer, in (b) vertically polarized illumination and a vertical
polarizer, and (c) between left circular polarization illumination and a right
circular polarizer denoted between each image pair. In (d-f) the same
experiments is repeated at 662 nm. The MM images were normalized, so variation
in the intensity transmitted through the second polarizer is due to
diattenuation, depolarization, and geometric transformation of the
illumination polarization. If the estimate of $\xi_{0}$ were too close to 0.25
(_i.e._ extrapolating to the ideal depolarizer), then there would be no
variation and all three polarizer conditions would look flat. Overall, there
is good visual agreement, with regions of high and low intensity in the DRR
results tending to correspond to high and low intensity in the Stokes camera
extrapolated results. It is interesting to note that, despite the MJM models
having only residual circular polarization properties from the complex
refractive indices, the crossed circular polarizer images show some agreement.
This is because of the variation in the depolarization due to the $\xi_{0}$
estimate. There is inherent disagreement between the measurement simulated
from DRR data and from the linear Stokes data due to the slightly different
viewing geometry, resolution, and camera properties between the two
polarimeters. Qualitatively, the agreement between DRR and Stokes results is
lowest on the side of the leg. This is the region where the TD assumption is
weakest, as shown in Fig. 2. Additionally, this region of high curvature is
where misregistration between the modeled geometry shown in Fig. 4 (b,d,f) and
the geometry measured with the Stokes camera would cause large deviation in
the modeled and measured $\widehat{\mathbf{m}}_{0}$, reducing the quality of
the extrapolation.

##  VI Conclusions and Discussion

In contrast with precision optical components, depolarization is the dominant
polarimetric property of everyday materials. The TD-MM pBRDF model lends
itself to efficient descriptions for such depolarization-dominant materials.
This work presents a closed-form Mueller-matrix valued pBRDF that is
efficiently represented by a reduced parameter set: four wavelength-dependent
material constants that define the polarization properties and a parameter for
depolarization that depends on both geometry and wavelength. The material
constants determine how the six degrees of freedom for the diattenuation and
retardance vary with geometry. This reduced parameterization assumes a TD
coherency eigenspectrum and specifies two distinct non-depolarizing MJMs for
first-surface and diffuse polarized light scattering. A microfacet Fresnel MJM
is a popular polarized-light scattering model that is simple enough to
describe first-surface scattering in many cases. However, polarization from
diffuse scattering becomes non-negligible when a wider range of scattering
geometries are considered.

In this work, the TD-MM pBRDF model is demonstrated over a large range of
scattering geometries by combining polarized first-surface and diffuse
contributions into a single MJM model. The first-surface MJM is Fresnel
reflection based on microfacet theory, but purely radiometric effects such as
shadowing and masking are ignored. The polarimetric properties of microfacet
Fresnel reflections are invariant to the surface normal. However, the diffuse
MJM is a diattenuator with orientation that strongly depends on the surface
normal in a pattern centered where $\theta_{h}=0^{\circ}$ (_i.e._ specular
reflection). This orientation is parallel to the plane of incidence with the
caveat that 45∘ and 135∘ are flipped due to reflection. The magnitude of the
diffuse contribution is determined by material constants.

The four material constants depend chiefly on the spectrally-dependent albedo
of the material. The depolarization parameter depends on both the albedo and
scattering geometry. The error in the TD assumption is inversely proportional
to depolarization. For example, the red object used in this work is higher
albedo at 662 nm compared to 451nm and the TD assumption is more appropriate
for the high albedo case, as shown in Fig. 2.

Multi-angle MM measurements at two illumination wavelengths (for both high and
low albedo conditions) of a red 3D printed material are compared to our pBRDF
model. The RMSD in diattenuation orientation for the modeled versus measured
dominant MJM (see Fig. 9 (a)) averaged over acquisition geometry was 7.49∘ and
14.29∘ at 451 nm (low albedo) and 662 nm (high albedo), respectively. The RMSD
in diattenuation magnitude (see Fig. 9 (b)) averaged over acquisition geometry
was 4.96% at 451 nm and 11.73% at 662 nm. Next, as an example application, the
MJM term of our pBRDF model is evaluated at assumed values of the material
constants, see Table I, and the geometry-dependent depolarization parameter
and the average reflectance are estimated from only linear Stokes images.
Despite the MJM model diattenuation properties matching the measurement better
at 451 nm than at 662 nm, using the MJM models to estimate $\xi_{0}$ resulted
in RMSD values (see Fig. 11) averaged over acquisition geometry of 11.11% at
451 nm and 4.24% at 662 nm. This performance difference can be understood by
noting that the depolarization magnitude at 662 nm is higher and therefore the
estimation is less sensitive to MJM model accuracy.

To demonstrate the generalization of the TD-MM model to more complicated
geometry, MM images of a 3D printed Stanford bunny were extrapolated from
linear Stokes images. To compare extrapolated MM images to ground truth from
DRR polarimetry, Fig. 12 shows simulated measurements of the bunny between
various polarizers. Overall agreement is good, with regions of high and low
transmitted intensity in the extrapolated results corresponding to those of
the DRR results. The largest deviations occur on the left side of the hind leg
where Fig. 2 shows the TD eigenspectrum assumption is weakest as well as being
the region where registration between the geometry modeled (shown in Fig. 4
(b,d,f)) and the actual measured geometry is most sensitive.

Quantitative results are presented for the sphere and visual agreement is
assessed for the Stanford bunny, but the success of the model is application-
dependent.

The quantitative results for the sphere and the visual assessment performed on
the Stanford bunny highlight an important point: the accuracy of a
polarimetric model does not need to exceed the accuracy of the polarimeter in
use. For example, the polarimetric accuracy of a model required to approach
photorealism in a polarization-aware physics-based rendering is lower than the
accuracy required in silicon wafer metrology. In the former application, a
COTS Stokes polarimeter might be used to characterize a material based on an
efficient, reduced parameterization with sufficient accuracy. Facilitating the
characterization of many materials according to a simple model in a reduced
length of time may present a significant benefit over the cost of more
rigorous characterization.

The TD-MM pBRDF model is intended to generalize to materials other than the
red 3D printing material used for demonstration purposes in this work. By
combining first-surface and diffuse polarized reflection in amounts that are
material-dependent the model is designed to capture the most relevant features
of polarized light scattering. To apply the model in Eq. 4 to a given
material, the albedo-dependent material constants $n_{\lambda}$,
$a_{\lambda}$, and $b_{\lambda}$ must be determined. One method to estimate
these values is MM imaging at a range of scattering geometries and least-
squares fitting to the dominant process. This has the advantage of also
providing a ground truth of the eigenspectrum to gauge the appropriateness of
the TD assumption. Another approach, based on the _a priori_ TD assumption, is
to maximize the linear correlation coefficient between polarized measurements
and measurements simulated from the MJM model defined by the material
constants. Modulation in the polarimetric measurements of a TD MM, see Eq. 5,
is determined only by the $\widehat{\mathbf{m}}_{0}$ term, since the ideal
depolarizer $\mathbf{m}_{ID}$ is an offset. Therefore, polarimetric
measurements of the full depolarizing MM at a given scattering geometry are
linearly related to polarimetric measurements of the dominant MJM. The slope
depends on $\xi_{0}$ and $M_{00}$, which in turn each depend on scattering
geometry, but the material constants $n_{\lambda}$, $a_{\lambda}$, and
$b_{\lambda}$ are geometry-independent. Therefore, $n_{\lambda}$,
$a_{\lambda}$, and $b_{\lambda}$ could be globally tuned to optimize the
linear fit across the pixel-dependent scattering geometry.

## References

  * [1] W. Boerner, “Use of polarization in electromagnetic inverse scattering,” _Radio science_ , vol. 16, no. 6, p. 1037–1045, 1981, articleID:1S0304.
  * [2] K. I. Hopcraft, K. I. Smith, and P. R. Smith, _An introduction to electromagnetic inverse scattering_ , 1st ed. Dordrecht: Springer Netherlands, 1992.
  * [3] F. O. Bartell, E. L. Dereniak, and W. L. Wolfe, “The theory and measurement of bidirectional reflectance distribution function (BRDF) and bidirectional transmittance distribution function (BTDF),” in _Radiation Scattering in Optical Systems_ , G. H. Hunt, Ed., vol. 0257, International Society for Optics and Photonics. SPIE, 1981, pp. 154 – 160. [Online]. Available: https://doi.org/10.1117/12.959611
  * [4] H. Mueller, “Proceedings of the winter meeting of the optical society of america,” _J. Opt. Soc. Am._ , vol. 38, no. 7, p. 661, Jul 1948. [Online]. Available: https://opg.optica.org/abstract.cfm?URI=josa-38-7-657
  * [5] N. G. Parke, “Optical algebra,” _Journal of Mathematics and Physics_ , vol. 28, no. 1-4, p. 131–139, 1949.
  * [6] R. Chipman, W. Lam, and G. Young, _Polarized Light and Optical Systems_ , ser. Optical Sciences and Applications of Light. CRC Press, 2018. [Online]. Available: https://books.google.com/books?id=KwFlDwAAQBAJ
  * [7] T. A. Germer, “pySCATMECH: a Python interface to the SCATMECH library of scattering codes,” in _Reflection, Scattering, and Diffraction from Surfaces VII_ , L. M. Hanssen, Ed., vol. 11485, International Society for Optics and Photonics. SPIE, 2020, p. 114850J. [Online]. Available: https://doi.org/10.1117/12.2568578
  * [8] I. G. E. Renhorn, T. Hallberg, and G. D. Boreman, “Efficient polarimetric brdf model,” _Opt. Express_ , vol. 23, no. 24, pp. 31 253–31 273, Nov 2015\. [Online]. Available: https://opg.optica.org/oe/abstract.cfm?URI=oe-23-24-31253
  * [9] A. Wilkie and A. Weidlich, “How to write a polarisation ray tracer,” in _SIGGRAPH Asia 2011 Courses_ , ser. SA ’11. New York, NY, USA: Association for Computing Machinery, 2011. [Online]. Available: https://doi.org/10.1145/2077434.2077442
  * [10] S.-H. Baek, R. Raskar, J. Ye, A. Dave, A. Kadambi, and H. Chen, “Polarization-based visual computing,” in _SIGGRAPH 2023 Courses_. New York, NY, USA: Association for Computing Machinery, 2023.
  * [11] A. Jarabo and D. Gutierrez, “Bidirectional rendering of polarized light transport,” in _Proceedings of the XXVI Spanish Computer Graphics Conference_ , ser. CEIG ’16. Goslar, DEU: Eurographics Association, 2016, p. 15–22.
  * [12] K. Berger and M. Blanchon, “Introduction to polarization for rendering and vision,” _SIGGRAPH Asia 2020 Courses_ , 2019.
  * [13] X. Wu, P. Li, X. Zhang, J. Chen, and F. Huang, “Three dimensional shape reconstruction via polarization imaging and deep learning,” _Sensors_ , vol. 23, no. 10, 2023. [Online]. Available: https://www.mdpi.com/1424-8220/23/10/4592
  * [14] S.-H. Baek, D. S. Jeon, X. Tong, and M. H. Kim, “Simultaneous acquisition of polarimetric SVBRDF and normals,” _ACM Trans. Graph._ , vol. 37, no. 6, dec 2018. [Online]. Available: https://doi.org/10.1145/3272127.3275018
  * [15] Y. Ba, A. Gilbert, F. Wang, J. Yang, R. Chen, Y. Wang, L. Yan, B. Shi, and A. Kadambi, “Deep shape from polarization,” in _Computer Vision – ECCV 2020_ , A. Vedaldi, H. Bischof, T. Brox, and J.-M. Frahm, Eds. Cham: Springer International Publishing, 2020, pp. 554–571.
  * [16] A. Dave, Y. Zhao, and A. Veeraraghavan, “PANDORA: Polarization-Aided Neural Decomposition Of Radiance,” _arXiv:2203.13458 [cs]_ , Mar. 2022, arXiv: 2203.13458. [Online]. Available: http://arxiv.org/abs/2203.13458
  * [17] C. Lei, C. Qi, J. Xie, N. Fan, V. Koltun, and Q. Chen, “Shape from polarization for complex scenes in the wild,” _arXiv: 2112.11377_ , 2021\.
  * [18] I. Hwang, D. S. Jeon, A. Muñoz, D. Gutierrez, X. Tong, and M. H. Kim, “Sparse ellipsometry: Portable acquisition of polarimetric SVBRDF and shape with unstructured flash photography,” _ACM Trans. on Graphics (Proc. SIGGRAPH 2022)_ , vol. 41, no. 4, 2022.
  * [19] L. B. Wolff, “Diffuse-reflectance model for smooth dielectric surfaces,” _J. Opt. Soc. Am. A_ , vol. 11, no. 11, pp. 2956–2968, Nov 1994. [Online]. Available: https://opg.optica.org/josaa/abstract.cfm?URI=josaa-11-11-2956
  * [20] K. E. Torrance and E. M. Sparrow, “Theory for off-specular reflection from roughened surfaces$\ast$,” _J. Opt. Soc. Am._ , vol. 57, no. 9, pp. 1105–1114, Sep 1967. [Online]. Available: http://www.osapublishing.org/abstract.cfm?URI=josa-57-9-1105
  * [21] R. Priest and T. Germer, “Polarimetric BRDF in the microfacet model: Theory and measurements.” Proceedings of the Military Sensing Symposia (MSS) Specialty Group Meeting on Passive Sensors, Undefined, 2000-03-01 00:03:00 2000.
  * [22] L. Wolff, “Using polarization to separate reflection components,” in _1989 IEEE Computer Society Conference on Computer Vision and Pattern Recognition_. Los Alamitos, CA, USA: IEEE Computer Society, jun 1989, pp. 363,364,365,366,367,368,369. [Online]. Available: https://doi.ieeecomputersociety.org/10.1109/CVPR.1989.37873
  * [23] G. Atkinson and E. Hancock, “Recovery of surface orientation from diffuse polarization,” _IEEE transactions on image processing_ , vol. 15, no. 6, p. 1653–1664, 2006.
  * [24] Z. Cui, J. Gu, B. Shi, P. Tan, and J. Kautz, “Polarimetric multi-view stereo,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 369–378.
  * [25] A. Kadambi, V. Taamazyan, B. Shi, and R. Raskar, “Polarized 3d: High-quality depth sensing with polarization cues,” in _2015 IEEE International Conference on Computer Vision (ICCV)_ , 2015, pp. 3370–3378.
  * [26] S.-H. Baek, D. S. Jeon, X. Tong, and M. H. Kim, “Simultaneous acquisition of polarimetric SVBRDF and normals,” _ACM Trans. Graph._ , vol. 37, no. 6, dec 2018. [Online]. Available: https://doi-org.ezproxy3.library.arizona.edu/10.1145/3272127.3275018
  * [27] Y. Kondo, T. Ono, L. Sun, Y. Hirasawa, and J. Murayama, “Accurate polarimetric BRDF for real polarization scene rendering,” in _European Conference on Computer Vision_. Springer, 2020, pp. 220–236.
  * [28] S. R. Marschner, S. H. Westin, E. P. Lafortune, and K. E. Torrance, “Image-based bidirectional reflectance distribution function measurement,” _Applied Optics_ , vol. 39, no. 16, p. 2592–2600, 2000.
  * [29] W. Matusik, H. Pfister, M. Brand, and L. McMillan, “A data-driven reflectance model,” _ACM Transactions on Graphics_ , vol. 22, no. 3, p. 759–769, 2003.
  * [30] S.-H. Baek, T. Zeltner, H. Ku, I. Hwang, X. Tong, W. Jakob, and M. Kim, “Image-based acquisition and modeling of polarimetric reflectance,” _ACM transactions on graphics_ , vol. 39, no. 4, p. 139:1–139:14, 2020.
  * [31] r. Ossikovski and J. Gil, _Polarized Light and the Mueller Matrix Approach_. CRC Press, 2016.
  * [32] N. Hagen, K. Oka, and E. L. Dereniak, “Snapshot Mueller matrix spectropolarimeter,” _Optics Letters_ , vol. 32, no. 15, p. 2100, Aug 2007.
  * [33] M. Dubreuil, S. Rivet, B. L. Jeune, and J. Cariou, “Snapshot Mueller matrix polarimeter by wavelength polarization coding,” _Opt. Express_ , vol. 15, no. 21, pp. 13 660–13 668, Oct 2007. [Online]. Available: http://opg.optica.org/oe/abstract.cfm?URI=oe-15-21-13660
  * [34] P. Lemaillet, S. Rivet, and B. L. Jeune, “Optimization of a snapshot Mueller matrix polarimeter,” _Opt. Lett._ , vol. 33, no. 2, pp. 144–146, Jan 2008\. [Online]. Available: http://opg.optica.org/ol/abstract.cfm?URI=ol-33-2-144
  * [35] A. S. Alenin and J. S. Tyo, “Task-specific snapshot Mueller matrix channeled spectropolarimeter optimization,” in _Polarization: Measurement, Analysis, and Remote Sensing X_ , D. B. Chenault and D. H. Goldstein, Eds., vol. 8364, International Society for Optics and Photonics. SPIE, 2012, p. 836402. [Online]. Available: https://doi.org/10.1117/12.921911
  * [36] R. M. A. Azzam, “Photopolarimetric measurement of the Mueller matrix by fourier analysis of a single detected signal,” _Optics Letters_ , vol. 2, no. 6, p. 148, Jun 1978\.
  * [37] D. H. Goldstein, “Mueller matrix dual-rotating retarder polarimeter,” _Applied Optics_ , vol. 31, no. 31, p. 6676–6683, Nov 1992.
  * [38] J. M. López-Téllez, R. A. Chipman, L. W. Li, S. C. McEldowney, and M. H. Smith, “Broadband extended source imaging Mueller-matrix polarimeter,” _Opt. Lett._ , vol. 44, no. 7, pp. 1544–1547, Apr 2019. [Online]. Available: http://ol.osa.org/abstract.cfm?URI=ol-44-7-1544
  * [39] M. Kupinski, M. Boffety, F. Goudail, R. Ossikovski, A. Pierangelo, J. Rehbinder, J. Vizet, and T. Novikova, “Polarimetric measurement utility for pre-cancer detection from uterine cervix specimens,” _Biomed. Opt. Express_ , vol. 9, no. 11, pp. 5691–5702, Nov 2018. [Online]. Available: http://opg.optica.org/boe/abstract.cfm?URI=boe-9-11-5691
  * [40] J. Qi, M. Ye, M. Singh, N. T. Clancy, and D. S. Elson, “Narrow band 3 × 3 mueller polarimetric endoscopy,” _Biomed. Opt. Express_ , vol. 4, no. 11, pp. 2433–2449, Nov 2013. [Online]. Available: https://opg.optica.org/boe/abstract.cfm?URI=boe-4-11-2433
  * [41] T. Novikova and J. C. Ramella-Roman, “Is a complete mueller matrix necessary in biomedical imaging?” _Opt. Lett._ , vol. 47, no. 21, pp. 5549–5552, Nov 2022. [Online]. Available: https://opg.optica.org/ol/abstract.cfm?URI=ol-47-21-5549
  * [42] M. Gonzalez, R. Ossikovski, T. Novikova, and J. C. Ramella-Roman, “Introduction of a 3 × 4 Mueller matrix decomposition method,” _Journal of Physics D: Applied Physics_ , vol. 54, no. 42, p. 424005, aug 2021\. [Online]. Available: https://doi.org/10.1088/1361-6463/ac1622
  * [43] M. K. Kupinski, J. Bankhead, A. Stohn, and R. Chipman, “Binary classification of mueller matrix images from an optimization of poincaré coordinates,” _J. Opt. Soc. Am. A_ , vol. 34, no. 6, pp. 983–990, Jun 2017. [Online]. Available: https://opg.optica.org/josaa/abstract.cfm?URI=josaa-34-6-983
  * [44] L. Li and M. Kupinski, “Merit functions and measurement schemes for single parameter depolarization models,” _Opt. Express_ , vol. 29, no. 12, pp. 18 382–18 407, Jun 2021. [Online]. Available: http://www.opticsexpress.org/abstract.cfm?URI=oe-29-12-18382
  * [45] S. R. Cloude, “Group theory and polarisation algebra,” _Optik_ , vol. 75, pp. 26–36, 1986.
  * [46] S. R. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering,” _Optical Engineering_ , vol. 34, no. 6, pp. 1599 – 1610, 1995\. [Online]. Available: https://doi.org/10.1117/12.202062
  * [47] C. Brosseau, _Fundamentals of polarized light: a statistical optics approach_. New York: Wiley, 1998, includes bibliographical references and indexes.
  * [48] R. Ossikovski and O. Arteaga, “Integral decomposition and polarization properties of depolarizing mueller matrices,” _Opt. Lett._ , vol. 40, no. 6, pp. 954–957, Mar 2015. [Online]. Available: https://opg.optica.org/ol/abstract.cfm?URI=ol-40-6-954
  * [49] K. Omer and M. Kupinski, “Compression, interpolation, and importance sampling for polarized brdf models,” _Opt. Express_ , vol. 30, no. 14, pp. 25 734–25 752, Jul 2022. [Online]. Available: http://opg.optica.org/oe/abstract.cfm?URI=oe-30-14-25734
  * [50] Q. Jarecki and M. Kupinski, “Extrapolating Mueller matrices from linear Stokes images,” in _Polarization: Measurement, Analysis, and Remote Sensing XV_ , D. B. Chenault and M. K. Kupinski, Eds., vol. 12112, International Society for Optics and Photonics. SPIE, 2022, p. 121120D. [Online]. Available: https://doi.org/10.1117/12.2619055
  * [51] Q. Jarecki and M. K. Kupinski, “Underdetermined polarimetric measurements for Mueller extrapolations,” _Optical Engineering_ , vol. 61, no. 12, p. 123104, 2022. [Online]. Available: https://doi.org/10.1117/1.OE.61.12.123104
  * [52] v. N. Umow, “Chromatische depolarisation durch lichtzerstreuung,” _Physikalische Zeitschrift_ , vol. 6, no. 20, pp. 674–676, 1905.
  * [53] M. K. Kupinski, C. L. Bradley, D. J. Diner, F. Xu, and R. A. Chipman, “Angle of linear polarization images of outdoor scenes,” _Optical Engineering_ , vol. 58, no. 8, pp. 1 – 12, 2019. [Online]. Available: https://doi.org/10.1117/1.OE.58.8.082419
  * [54] B. H. Billings, “A monochromatic depolarizer$\ast$†,” _J. Opt. Soc. Am._ , vol. 41, no. 12, pp. 966–975, Dec 1951. [Online]. Available: http://www.osapublishing.org/abstract.cfm?URI=josa-41-12-966
  * [55] J. J. Gil and E. Bernabeu, “Depolarization and polarization indices of an optical system,” _Optica Acta: International Journal of Optics_ , vol. 33, no. 2, pp. 185–189, 1986. [Online]. Available: https://doi.org/10.1080/713821924
  * [56] R. Ossikovski and J. Vizet, “Eigenvalue-based depolarization metric spaces for Mueller matrices,” _J. Opt. Soc. Am. A_ , vol. 36, no. 7, pp. 1173–1186, Jul 2019. [Online]. Available: http://josaa.osa.org/abstract.cfm?URI=josaa-36-7-1173
  * [57] S. M. Rusinkiewicz, _A New Change of Variables for Efficient BRDF Representation_. Vienna: Springer Vienna, 2011, p. 11–22.
  * [58] B. Walter, S. R. Marschner, H. Li, and K. E. Torrance, “Microfacet models for refraction through rough surfaces.” _Rendering techniques_ , vol. 2007, p. 18th, 2007.
  * [59] F. Bréon, “An analytical model for the cloud-free atmosphere/ocean system reflectance,” _Remote sensing of environment_ , vol. 43, no. 2, pp. 179–192, 1993.
  * [60] M. Ashikmin, S. Premože, and P. Shirley, “A microfacet-based brdf generator,” in _Proceedings of the 27th annual conference on Computer graphics and interactive techniques_ , 2000, pp. 65–74.
  * [61] F.-M. Breon and F. Maignan, “A brdf–bpdf database for the analysis of earth target reflectances,” _Earth System Science Data_ , vol. 9, no. 1, pp. 31–45, 2017.
  * [62] E. Heitz, “Understanding the masking-shadowing function in microfacet-based BRDFs,” _Journal of Computer Graphics Techniques Vol_ , vol. 3, no. 2, 2014\.

##  VII Biography Section

|  Meredith Kupinski received the M.S. and Ph.D. degrees in optical science
from The University of Arizona, Tucson, AZ, USA, in 2003 and 2008,
respectively. She is currently an Assistant Professor with the Wyant College
of Optical Sciences, The University of Arizona. Her research interests include
task-relevant metrics for polarimeter design, stochastic systems analysis and
information theory.
---|---
|  Quinn Jarecki received the B.S. and B.M. degrees in optical science and
music performance, respectively, from The University of Arizona, Tucson, AZ,
USA, in 2015. He is currently a Ph.D candidate at the Wyant College of Optical
Sciences, The University of Arizona. His research interests include full and
partial Mueller polarimetric imaging and polarized light scattering.
---|---