April 15, 2012 / Vol. 37, No. 8 / OPTICS LETTERS

1367

Snapshot imaging Mueller matrix polarimeterusing polarization gratings
Michael W. Kudenov,1,* Michael J. Escuti,2 Nathan Hagen,3 Eustace L. Dereniak,1 and Kazuhiko Oka4
1

College of Optical Science, The University of Arizona, 1630 E. University Boulevard, Tucson, Arizona 85721, USA
2
Department of Electrical & Computer Engineering, North Carolina State University,
2410 Campus Shore Drive, Raleigh, North Carolina 27606, USA
43
Rebellion Photonics, 7457 South Freeway, Houston, Texas 77021, USA
Division of Applied Physics, Hokkaido University, N-13, W-8, Sapporo 060-8628, Japan*Corresponding author: mkudenov@optics.arizona.edu

Received December 15, 2011; revised February 24, 2012; accepted February 28, 2012;
posted February 29, 2012 (Doc. ID 160054); published April 12, 2012
A snapshot imaging Mueller matrix polarimeter (SIMMP) is theoretically described and empirically demonstrated
through simulation. Spatial polarization fringes are localized onto a sample by incorporating polarization gratings(PGs) into a polarization generator module. These fringes modulate the Mueller matrix (MM) components of thesample, which are subsequently isolated with PGs in an analyzer module. The MM components are amplitude
modulated onto spatial carrier frequencies which, due to the PGs, maintain high visibility in spectrally broadband
illumination. An interference model of the SIMMP is provided, followed by methods of reconstruction andcalibration. Lastly, a numerical simulation is used to demonstrate the system’s performance in the presence ofnoise. © 2012 Optical Society of America
OCIS codes: 110.5405, 120.5410, 260.3160.

SIMMP Analyzer

P 1 L1 L2

L3 L4

L6

L5

L7

L8 P 2

FPA

SourceSIMMP Generator

SampleM(x, y)
x f1
y z

t

t
y

f2

f3

x z

2t

2t
QWP2

0146-9592/12/081367-03$15.00/0

T represents the transpose operation and S 0;in is the total
incident power. Transmission of this light through PGs L1
through L4 enables the second re-imaging lens (f 2 ) to localize polarization interference fringes [10,11] onto thesample. These fringes modulate the Mueller matrixMx; y of the sample by two spatial carrier frequencies,
the field of which is then collimated by lens f 3 into theanalyzing optics. It should be mentioned that Mx; y isassumed to have no wavelength dependence over the
spectral range of the measurement. The 4 beams fromthe generator are then sheared into 16 beams by the analyzing PGs L5 through L8 . Lastly, the linear polarizer (P 2 )enables intensity fringes to be localized onto the focal
plane array (FPA) by the re-imaging lens f 4 .
The functional form of the intensity distribution on
the FPA was calculated using a combination of Mueller
calculus and scalar diffraction theory [8,12]. The irradiance at the FPA, assuming that f 1  f 2  f 3  f 4 ,becomes

QWP1

Mueller matrix (MM) imaging polarimetry is capableof characterizing the diattenuation, retardation, and depolarization properties of a surface, optical element, or
specimen [1]. As such, it has applications in ground-truth
for remote sensing [2], quality control [3], and biomedical
imaging [4]. A sample’s MM is often characterized bytime-sequentially changing the polarization state of both
the generated and analyzed light using rotating polarization elements. By measuring the intensity under at least
16 unique combinations of generated and analyzed
polarization states, the 16 Mueller matrix components
can be calculated. However, more measurements are often taken to increase the measured signal-to-noise ratio(SNR) and to optimize the condition number of the system’s measurement matrix [4,5].
Since a conventional MM imaging polarimeter uses
temporal scanning, intensity variations caused by sample
motion must be minimized and often require an image
registration scheme [6]. Misregistration is also causedby systematic errors, such as beam-wander from the rotating optical elements. This further increases the complexity in applications where MM imaging polarimetry
may be beneficial, especially for in vivo biomedical imaging of scenes with low inherent contrast [6].
To eliminate the need for image registration andcapture all information in a single measurement (i.e., asnapshot), we propose to combine a technique forchanneled Mueller matrix spectropolarimetry [7] with
an imaging channeled polarimeter based on polarization
gratings (PGs) [8,9]. The layout of the snapshot imaging
Mueller matrix polarimeter (SIMMP) is depicted inFig. 1. It contains both a generator and analyzer module.
In the generator module, the source is either spatiallyincoherent quasi-monochromatic light or spectrally
broadband illumination that is collimated by lens f 1 .
Transmission of this light through the linear polarizer(P 1 ) makes the incident normalized Stokes vector
Sin  S 0;in x; y∕2 1 0 1 0 T , where the superscript

y

f4

x z

Fig. 1. (Color online) SIMMP optical configuration. PGs L1 , L2 ,L5 and L6 shear the beam along x while L3 , L4 , L7 and L8 shear
along y. P 1 and P 2 are linear polarizers at 45º while two quarterwave-plates, QWP1 and QWP2 , have fast axes oriented at 45º
and 0º, respectively. All PGs have identical grating periods Λand the generator’s and analyzer’s PGs are separated by a distance t and 2t, respectively.
© 2012 Optical Society of America

1368

OPTICS LETTERS / Vol. 37, No. 8 / April 15, 2012

In Table 1, the MM elements are represented as mij ,where the subscripts i and j represent the row and column of the 4 × 4 element Mueller matrix M. All MM elements are implicitly dependent on x and y for clarity.
Fourier transformation of Ix; y produces 33 channels
in I 0 ξ; η, the Fourier domain representation of the measured irradiance, where each channel is proportional toone of the coefficients in Table 1. A map of the channels
in the Fourier domain is depicted in Fig. 2. Extractionof these channels is identical to [8, 12–14], in which a2-dimensional filter is placed at each channel to isolate it.
Once a channel is isolated, its content is inverse Fourier transformed to produce k filtered channels, C k , thatare proportional to Ak modulated by a complex exponential phase. To demodulate the phase from Ak , each channel’s phase is quantified using a previously recorded
reference taken of one or more known Mueller matrices
Mr . From Mr , the reference light’s Ak coefficients (Ark )are calculated. The kth channel is demodulated and calibrated byIx; y  jA1 j cosκαx  3y  Aa1  jA2 j cosκα3x − 3y  Aa2  jA3 j cosκαx  y  Aa3  jA4 j cosκα3x − y  Aa4  jA5 j cosκαx − y  Aa5  jA6 j cosκα3x  y  Aa6  jA7 j cosκαx − 3y  Aa7  jA8 j cosκα3x  3y  Aa8  jA9 j cosκαx  2y  Aa9 
 jA10 j cosκα3x − 2y  Aa10  jA11 j cosκαx − 2y  Aa11  jA12 j cosκα3x  2y  Aa12  jA13 j cosκα2x − 2y  Aa13  jA14 j cosκα2x  2y  Aa14 

C k;cal  Ark C sk ∕C rk C r17 ∕C s17 ; jA15 j cosκαx  Aa15  jA16 j cosκα2y  Aa16   A17 ;(1)

where κ  2π∕λf 4 is the phase constant, λ is the freespace wavelength, and the generator’s and analyzer’s
shear are α  2tλ∕Λ and β  4tλ∕Λ, respectively, whereΛ is the PG’s period. Note that κα and κβ are constantwith wavelength meaning that the carrier frequencies
are independent of the illumination’s temporal coherence
length. This is directly analogous to [8] in that the shear
depends linearly on the wavelength. Meanwhile, thecomplex coefficients Ak x; y represent summations ofMueller matrix elements while Aak  argAk x; y. The
Ak coefficients are shown in Table 1 and implicitly
depend on x and y for clarity.
Table 1. Coefficient Definitions for the Intensity
Pattern
Ak 

Coefficient × S 0;in x; y∕16

A1 
A2 
A3 

−m23 − m32 − im22  im33  2m42  2im43−m23  m32  im22  im33
−m23  m32  im22  im33 − 2m42
− 2im43 − 4m13  4im12−m23 − m32 − im22  im33
−m23 − m32 − im22  im33 − 2m42
 2im43 − 4m13  4im12−m23  m32  im22  im33
−m23  m32  im22  im33  2m42  2im43−m23 − m32 − im22  im33
−2m34 − 2im24 − 4m44
−2m34  2im24−2m34 − 2im24  4m44
2m34  2im24
4m21 − 4im31
4m21 − 4im31
8im14
8im41
8m11

A4 =
A5 
A6 
A7 
A8 
A9 
A10 
A11 
A12 
A13 
A14 
A15 
A16 
A17 

(2)

where C sk , C rk are the kth sample and reference channels, respectively, and C r17 , C s17 are included to remove
any illumination dependencies in S 0;in between the reference and sample. Once demodulated, the individually calibrated channels must be algebraically combined toisolate the MM components per Table 2.
The intensity pattern of Eq. (1) was simulated assuming a spatially uniform input where all mij  0.25. A Poisson-noise-limited 1024 × 1024 pixel array was simulated
for use as the imaging sensor. The noise was modeledusing laboratory measurements taken from an 8-bit
Imaging Source 41BU02 machine vision camera at room
temperature (gain  456, exposure  1∕83 s). The system’s simulated parameters are t  5.1 mm, f 1  f 2 
f 3  f 4  50 mm, λ  0.45 − 0.65 μm, and Λ  8 μm.
These parameters enable a cutoff frequency, for eachchannel, of 10.3 cycles∕mm. Sample measurements involved no time-averaging while simulated calibration
images were averaged over 20 frames. To ensure a
non-zero output within each channel, calibration data
were taken assuming two reference targets: (1) a linear
polarizer at 0 degrees followed by a quarter-wave retarder at 22.5 degrees (R1) and (2) a quarter-wave retarderFig. 2. (Color online) Fourier domain of a channeledimage obtained from the SIMMP. Channel numbers correspondto the k subscripts of the Ak coefficients per Table 1. Only the
non-conjugated channels are numbered for clarity.

April 15, 2012 / Vol. 37, No. 8 / OPTICS LETTERS
Table 2. Mueller Matrix Solutions from theFourier domain
mij

Channel Combinations × jC r17 j∕jC s17 j

m11
m12

jC s17 j∕jC r17 j−2 ImAr2 C s2 ∕C r2 − Ar3 C s3 ∕C r3 − Ar5 C s5 ∕C r5
Ar8 C s8 ∕C r8 
−2 ReAr1 C s1 ∕C r1  Ar3 C s3 ∕C r3  Ar5 C s5 ∕C r5
Ar7 C s7 ∕C r7  − m23 ∕2−2 ImAr15 C s15 ∕C r15 
4 ReAr13 C s13 ∕C r13  Ar14 C s14 ∕C r14 
4 ImAr2 C s2 ∕C r2 − Ar4 C s4 ∕C r4  Ar6 C s6 ∕C r6−Ar8 C s8 ∕C r8 
−4 ReAr2 C s2 ∕C r2  Ar4 C s4 ∕C r4  Ar6 C s6 ∕C r6
Ar8 C s8 ∕C r8 
−4 ImAr9 C s9 ∕C r9  Ar11 C s11 ∕C r11 −2 ImAr13 C s13 ∕C r13  Ar14 C s14 ∕C r14 −4 ReAr1 C s1 ∕C r1 − Ar3 C s3 ∕C r3  Ar5 C s5 ∕C r5−Ar7 C s7 ∕C r7 
4 ImAr2 C s2 ∕C r2  Ar4 C s4 ∕C r4  Ar6 C s6 ∕C r6
Ar8 C s8 ∕C r8 
4 ReAr10 C s10 ∕C r10  Ar12 C s12 ∕C r12 
2 ImAr16 C s16 ∕C r16 
4 ReAr1 C s1 ∕C r1 − Ar4 C s4 ∕C r4 − Ar6 C s6 ∕C r6
Ar7 C s7 ∕C r7 
−4 ImAr1 C s1 ∕C r1 − Ar4 C s4 ∕C r4 − Ar6 C s6 ∕C r6
Ar7 C s7 ∕C r7 
−2 ReAr9 C s9 ∕C r9 − Ar11 C s11 ∕C r11 

m13
m14
m21
m22
m23
m24
m31
m32
m33
m34
m41
m42
m43
m44

at 22.5 degrees followed by a linear polarizer at 0 degrees(R2). Here, reference R1 is used to calibrate channels
1–15 and 17 while R2 is used to calibrate channel 16.
First, the MTF was quantified by calculating thecontrast of spatial frequencies spanning 0.1 to10.3 cycles∕mm. The MTF is provided in Fig. 3(a) andis primarily influenced by the shape of a band-limiting
window that was used to prevent cross talk between adjacent channels. Such band-limiting can be physicallyrealized by defocusing the sample while maintaining focus on the fringe localization field from the generator [8].
However, for the purposes of this simulation we used a
Hamming window with a full-width at half-maximum frequency of 5.4 cycles∕mm.
Each reconstructed MM element’s SNR was also calculated and is depicted in Fig. 3(b). Generally, MM components that occupy their own channel (e.g., m14 and m41 )have larger SNRs than MM components in highly multiplexed channels (e.g., m22 , m23 , m33 , m42 and m43 ).
The spatial output of a quarter-wave vortex retarder,
on a 4.8 × 4.8 mm square substrate, was also simulated[15]. A depiction of the ideal Mueller matrix and the
80
70

1

SNR

MTF

0.8
0.6

60
50

0.4

40
30

0.2

20
10

0
2

3

4

5

6

7

8

9

Spatial Frequency (cycles/mm)

(a)

10

m44
m43
m42
m41
m34
m33
m32
m31
m24
m23
m22
m21
m14
m13
m12

0
1

(b)

Fig. 3. (Color online) (a) Average MTF response of all 17 channels where error bars represent one standard deviation, (b) the
SNR of each Mueller matrix element.

Simulated Input

1369

Simulated Measurement

m11

m12

m13

m14

m11

m12

m13

m14

m21

m22

m23

m24

m21

m22

m23

m24

m31

m32

m33

m34

m31

m32

m33

m34

m41

m42

m43

m44

m41

m42

m43

m44

-1

-0.8 -0.6 -0.4 -0.2

0

0.2

0.4

0.6

0.8

1

Fig. 4. (Color online) Simulated input (left) and measured(right) Mueller matrix of the quarter-wave vortex retarder.

simulated output are provided in Figs. 4(a) and 4(b), respectively. From the simulated measured data, it is clearthat spatial frequencies near the center of the vortex are
not discernable, while lower spatial frequency information is preserved. This relates to the spatial band limiting
discussed previously.
The SIMMP approach takes advantage of the rapid advance in detector array technology, allowing users to uselarge pixel count arrays in order to achieve fast readout
rates. This allows a 16-fold improvement in measurement
speed at the cost of a 7-fold reduction in image resolution(in each x and y direction). Such advantages may be beneficial in biomedical imaging applications for real-time
diagnostics or for imaging a large number of samples
in rapid succession for quality control applications. Use
of PGs offers a lower cost, a more compact size, andlarger apertures over birefringent crystal optics withthe capability of using spectrally broadband illumination.
References
1. D. H. Goldstein, Polarized Light (CRC, 2011).
2. B. D. Bartlett, M. G. Gartley, D. W. Messinger, C. Salvaggio,
and J. R. Schott, J. Appl. Remote Sens. 4, 043552 (2010).
3. F. Delplancke, Appl. Opt. 36, 5388 (1997).
4. K. M. Twietmeyer, R. A. Chipman, A. E. Elsner, Y. Zhao, andD. VanNesdale, Opt. Express 16, 21339 (2008).
5. S. Murali, “Analysis of an interferometric Stokes imaging
polarimeter,” Ph. D. dissertation (University of Arizona,2010).
6. S. Guyot, M. Anastasiadou, E. Delechelle, and A. D. Martino,
Opt. Express 15 7393 (2007).
7. N. Hagen, K. Oka, and E. L. Dereniak, Opt. Lett. 32, 2100
(2007).
8. M. W. Kudenov, M. J. Escuti, E. L. Dereniak, and K. Oka,
Appl. Opt. 50, 2283 (2011).
9. C. Oh and M. J. Escuti, Opt. Lett. 33, 2287 (2008).
10. S. Pancharatnam, Proc. Indian Acad. Sci. A 44, 247 (1956).
11. M. Mujat, A. Dogariu, and E. Wolf, J. Opt. Soc. Am. A 21,
2414 (2004).
12. K. Oka and N. Saito, Proc. SPIE 6295, 629508 (2006).
13. K. Oka and T. Kaneko, Opt. Express 11, 1510 (2003).
14. M. W. Kudenov, M. E. L. Jungwirth, E. L. Dereniak, and G. R.
Gerhart, Opt. Express 17, 22520 (2009).
15. S. C. McEldowney, D. M. Shemo, R. A. Chipman, and P. K.
Smith, Opt. Lett. 33, 134 (2008).