# Infrared full-Stokes polarimetry by parametric up-conversion

Zhanghang Zhu† The Key Laboratory of Weak-Light Nonlinear Photonics, Ministry
of Education, School of Physics and TEDA Applied Physics Institute, Nankai
University, Tianjin, 300071, P.R. China  Di Zhang† The Key Laboratory of Weak-
Light Nonlinear Photonics, Ministry of Education, School of Physics and TEDA
Applied Physics Institute, Nankai University, Tianjin, 300071, P.R. China  Fei
Xie  The Key Laboratory of Weak-Light Nonlinear Photonics, Ministry of
Education, School of Physics and TEDA Applied Physics Institute, Nankai
University, Tianjin, 300071, P.R. China  Jiaxin Chen  The Key Laboratory of
Weak-Light Nonlinear Photonics, Ministry of Education, School of Physics and
TEDA Applied Physics Institute, Nankai University, Tianjin, 300071, P.R. China
Shengchao Gong  The Key Laboratory of Weak-Light Nonlinear Photonics, Ministry
of Education, School of Physics and TEDA Applied Physics Institute, Nankai
University, Tianjin, 300071, P.R. China  Wei Wu  The Key Laboratory of Weak-
Light Nonlinear Photonics, Ministry of Education, School of Physics and TEDA
Applied Physics Institute, Nankai University, Tianjin, 300071, P.R. China  Wei
Cai  The Key Laboratory of Weak-Light Nonlinear Photonics, Ministry of
Education, School of Physics and TEDA Applied Physics Institute, Nankai
University, Tianjin, 300071, P.R. China  Xinzheng Zhang  The Key Laboratory of
Weak-Light Nonlinear Photonics, Ministry of Education, School of Physics and
TEDA Applied Physics Institute, Nankai University, Tianjin, 300071, P.R. China
Mengxin Ren<EMAIL_ADDRESS>The Key Laboratory of Weak-Light
Nonlinear Photonics, Ministry of Education, School of Physics and TEDA Applied
Physics Institute, Nankai University, Tianjin, 300071, P.R. China
Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan,
Shanxi 030006, P.R. China  Jingjun Xu<EMAIL_ADDRESS>The Key Laboratory
of Weak-Light Nonlinear Photonics, Ministry of Education, School of Physics
and TEDA Applied Physics Institute, Nankai University, Tianjin, 300071, P.R.
China

###### Abstract

Polarimetry aims to measure polarization information of an optical field,
providing a new freedom to enhance performance of optical metrology.
Furthermore, the polarimetry in infrared (IR) range holds a promise for a wide
range of academic and industrial applications because the IR light relates to
unique spectral signatures of various complex chemicals. However, a primary
challenge in IR applications remains lacking efficient detectors. Motivated by
such a constraint, we present in this paper a nonlinear up-conversion full-
Stokes IR Retrieval Polarimetry (IRP) that is able to decipher the IR
polarization using a well-commercial and high-performance visible light
detector. Assisted by a nonlinear sum-frequency generation (SFG) process in a
lithium niobate thin film, the polarization states of the IR light are encoded
into the visible SF wave. Based on a Stokes-Mueller formalism developed here,
the IR polarization is successfully retrieved from SF light with high
precision, and polarization imaging over the targets with either uniform or
non-uniform polarization distributions are demonstrated. Our results form a
fresh perspective for the design of novel advanced up-conversion polarimeter
for a wide range of IR metrological applications.

Infrared (IR) light is electromagnetic radiation with wavelength longer than
that of visible (VIS) light.[1] Although invisible to human naked eyes, the IR
light holds significant usefulness in various cutting-edge applications
including biological sensing,[2] astronomical observation,[3] and
environmental science,[4] owing to its capacity for unique chemical
specificity. Until now, the IR related applications are mainly challenged by
the lack of efficient detectors.[5] For example, becuase the IR sensitive
semiconductors have small band gap energy, their performance is easily
disturbed by ambient thermal fluctuation in terms of dark current and noise.
Furthermore, the available pixel number of the arrayed IR detectors is
limited. This renders the IR cameras a poor resolution compared with their
visible counterparts, such as complementary metal-oxide semiconductors (CMOS)
and charge-coupled devices (CCD).[6]

An appealing solution to above constraints is parametric frequency up-
conversion, such as sum-frequency generation (SFG).[7] In such architecture,
instead of directly detecting the IR fundamental frequency (FF) photons, the
IR signal (at frequency $\omega_{1}$) is firstly transferred to the visible
range with the help of another pump light ($\omega_{2}$) via a second order
nonlinear ($\chi^{(2)}$) material.[8] Then the generated SF photons
($\omega_{3}=\omega_{1}+\omega_{2}$, following the energy conservation law)
are detected by high-performance visible detectors. The first study of such
parametric up-conversion of IR light could be traced back to 1960s, when
Midwinter and Warner demonstrated a conversion of 1.7 $\mu$m radiation to 493
nm by a lithium niobate (LN) crystal.[9] One year later, such up-conversion
technique was further applied to realize IR up-conversion imaging.[10] Over
the past decades, the parametric up-conversion detection has received great
interests, and several technical improvements were made in various aspects,
such as finding new nonlinear materials for higher up-conversion
efficiency,[11] ameliorating sensitivity of the visible detectors,[12]
manufacturing high power IR sources at various wavelengths for target
illumination,[13, 14] optimizing system designs for better imaging field-of-
view and resolution.[15, 16]

Figure 1:  A schematic of IRP system. (a) Operation principle of the IRP
system. FF signal ($\omega_{1}$) and pump light ($\omega_{2}$) are up-
converted into SF wave ($\omega_{3}$) through a second order nonlinear
material. The nonlinear SFG process can be presented by a Mueller matrix
$\mathbf{M}\left(\chi^{(2)},\mathbf{S}^{\omega_{2}}\right)$, via which the
Stokes vectors of $\mathbf{S}^{\omega_{1}}$ can be calculated reversibly from
$\mathbf{S}^{\omega_{3}}$. (b) A sketch of the IRP setup. An optical pump-
probe system was built. The path difference between the pump and probe pulses
is compensated by a mechanical delay line. The pump polarization is vertically
polarized by a Glan-Taylor (GT) prism, and the polarization of the probe light
is prepared by a polarization state generator (PSG) consisting of a GT
polarizer, a half-wave plate (HWP) and a quarter-wave plate (QWP). The optical
beams are focused onto the LN film respectively afterwards by lens pairs with
focal lengths of 500 mm (L1, L2) and 100 mm (L3). The generated SF light is
collected by a lens (L4, 75 mm), followed by a polarization state analyzer
(PSA) comprising a rotating QWP, and a GT prism, and finally imaged onto a
visible CCD camera after a lens (L5, 200 mm). Physical orientation and the
shape of the polarization ellipse are defined by angles of azimuth ($\phi$)
and ellipticity ($\xi$), respectively. Positive values of $\phi$ and $\xi$
correspond to the clockwise rotation of the polarization azimuth and a right-
handed polarization ellipse, as observed against the propagation direction.
Intensity images by the signal light for ‘L’, ‘N’ and ‘house’ masks are given
in the top-right corner.

However, the current IR up-conversion techniques are still designed to acquire
the intensity and spectral information of the IR light, while the polarization
information is simply dropped. Complementary to the intensity and spectrum,
the polarization is another fundamental characteristic that describes how the
light’s electric field vector oscillates. While the intensity and spectrum
tell us about the geometries and material compositions, polarization
interrogation can bring numerous additional richness from the targets.
Specifically, by analyzing the light polarization changes by scattering or
transmitting through the objects, people can infer birefringence or stress
inside material volumes, and roughness or texture of a surface, moreover,
enhance contrast for objects that are difficult to distinguish otherwise.[17,
18] In this context, polarimetry, as a technique invented to measure the
polarization state of light, has become an emerging technique to enhance many
fields of optical metrology, optical sensing, and so on.[19]

Enlightened by the respective merits of the two techniques of the nonlinear IR
up-conversion and the polarimetry, it would be a good idea to further combine
them together. Such incorporation could inherit advantages of both techniques,
i.e. detecting IR polarization information while using the VIS detectors. This
would consequently benefit many fields of IR related applications. But to the
best of our knowledge, such combination still remains unexplored. In this
paper, we experimentally demonstrate an IR Retrieval Polarimetry (IRP) that
manages to decipher the polarization information of the IR signal from the
generated SF light. We develop a Stokes-Mueller formalism based on nonlinear
optics that correlates the input IR light polarization states with the SF
light. In such algorithmic framework, the incoming IR signals and outgoing SF
radiations are represented by 4$\times$1 Stokes vectors, and the nonlinear SFG
process is represented by a 4$\times$4 Mueller matrix. We further
experimentally achieve the IR polarization imaging of different targets with
uniform or non-uniform polarization distributions using our IRP system, which
reconstructs the polarization of the IR signal light with high precision. Such
IRP technique would pave the way to the design of novel IR polarimetry for
advanced remote sensing, astronomy and industrial inspection.

##  I Stokes-Mueller formalism for IRP

First, we develop a Stokes-Mueller algebra for retrieving the polarization of
the IR signal light based on that of the SF wave, which essentially forms the
mathematical basis for the operation of our IRP. We adopted here a four-
element Stokes vector $\mathbf{S}=[s_{0},s_{1},s_{2},s_{3}]^{T}$ (T denotes
matrix transpose) to describe the polarization state of light.[20] The Stokes
description of light has a chief advantage that it operates in the form of
intensity, which is real valued and is easily measurable using photo-detectors
in experiments. Furthermore, in contrast to other representations, such as the
Jones vector that is limited to the case of fully polarized light,[21] the
Stokes vector also spans the space of unpolarized and partially polarized
light. Each Stokes element is defined as a sum or difference between
intensities of different orthogonal polarization bases. Specifically, $s_{0}$
is the total intensity of the light, $s_{1}$ is the intensity difference
between linear polarization components along 0∘ and 90∘, $s_{2}$ is the
intensity difference between $+45^{\circ}$ and $-45^{\circ}$ linear
polarizations, and $s_{3}$ is the intensity difference between right- and
left-handed circular polarization components. In such representation, the
polarization dependent interaction between the light and matter is denoted by
a 4$\times$4 Mueller matrix ($\mathbf{M}$), and the incident polarization
state $\mathbf{S}$ relates to the output state $\mathbf{S^{\prime}}$ via
$\mathbf{S^{\prime}}=\mathbf{M}\cdot\mathbf{S}$. Such Stokes-Mueller formalism
acts as a basis for the traditional linear optical polarimetry to characterize
the signal light polarization.[21] Here we further develop a new nonlinear
Stokes-Mueller formalism for the IRP technique, which aims to reconstruct the
IR polarization from that of the SF light, rather than measuring that of the
IR light directly, as illustrated by Fig. 1(a).

Figure 2:  Calibration of IRP. (a) Poincaré sphere of FF signal light (at
$\omega_{1}$). 72 different signal polarization states (red dots) that
uniformly sample the sphere were chosen as input polarizations for calibrating
IRP. (b) Poincaré sphere of generated SF light. (c) The three-dimensional FF
Poincaré sphere is expanded into a planar map in a coordinate of
azimuth–ellipticity angles ($\phi-\xi$). The originally input FF polarization
states are presented by the dots in color, and the reconstructed results are
given by small gray dots. The well overlap between them demonstrates the high
precision of our IRP in reconstructing the FF polarizations. (d) The
$\phi-\xi$ map of the SF Poincaré sphere. The SF dots generated by the FF
light are shown using the same color as those in (c).

According to the nonlinear optics, the generated SF waves are closely related
to the second order nonlinear polarizations $\mathbf{P}^{(2)}$ excited inside
the materials. Specifically, the $\mathbf{P}^{(2)}$ for the SFG process
follows[22]

$\small
P^{(2)}_{i}(\omega_{3})=2\sum_{j,k}\epsilon_{0}\chi_{ijk}^{(2)}E_{j}^{\omega_{1}}E_{k}^{\omega_{2}},$
(1)

in which footnotes $i,j,k$ stand for the unit vector along the Cartesian $x$-,
$y$\- or $z$-axes, and $E_{j}^{\omega_{1}}$ and $E_{k}^{\omega_{2}}$ are the
electric fields of the signal and pump light, respectively. $\chi_{ijk}^{(2)}$
is the element of the nonlinear susceptibility tensor, which generates
$\mathbf{P}^{(2)}$ vector polarized along the $i$-direction under the combined
interaction of the fundamental electric fields along $j$\- and $k$-directions.
Such $P^{(2)}_{i}$ acts as a secondary source to radiate the nonlinear SF
waves $\mathbf{E}^{\omega_{3}}$ polarized along the $i$-direction, and
magnitude of the $\mathbf{E}^{\omega_{3}}$ is proportional to
$\mathbf{P}^{(2)}$. By further considering the mathematical relationship
between $\mathbf{S}^{\omega_{3}}$ and $\mathbf{E}^{\omega_{3}}$, we could
derive

$\mathbf{S}^{\omega_{3}}=\mathbf{M}\left(\chi^{(2)},\mathbf{S}^{\omega_{2}}\right)\cdot\mathbf{S}^{\omega_{1}},$
(2)

in which $\mathbf{M}\left(\chi^{(2)},\mathbf{S}^{\omega_{2}}\right)$ is a
4$\times$4 Mueller matrix describing the SFG process. Equation (2) has a
similar form as the linear Stokes-Mueller formalism. However, different from
the Mueller matrix in linear optics that only relates to medium’s properties
such as birefringence or optical chirality,[21] the SFG matrix $\mathbf{M}$
here is a function of not only properties of the material (i.e. nonlinear
$\chi^{(2)}$), but also the pump polarization $\mathbf{S}^{\omega_{2}}$. Once
the SFG matrix $\mathbf{M}$ and Stokes parameters of the SF light
($\mathbf{S}^{\omega_{3}}$) are determined, the polarization state of the
signal light could be calculated via

$\mathbf{S}^{\omega_{1}}=\mathbf{M}^{-1}\left(\chi^{(2)},\mathbf{S}^{\omega_{2}}\right)\cdot\mathbf{S}^{\omega_{3}}.$
(3)

##  II Setup of IRP

An $x$-cut LN thin film was used as the SFG material in our experiment. The
thickness of the LN film is about 200 nm, which resides on a 500 $\mu$m thick
silicon dioxide (SiO2) substrate (NANOLN, Jinan Jingzheng Electronics Co.,
Ltd). LN is a typical optical material showing ultra-broadband transparency
and strong second order nonlinearity covering the visible to mid-IR spectral
range, which renders it a feasible material platform for wide photonic and
optoelectronic applications. The LN crystal belongs to a point group of $3m$.
Resulting from the symmetry restrictions, $\chi^{(2)}$ tensor of the LN has 11
non-vanishing components, in which only 4 elements are independent:
$\chi^{(2)}_{eee}$, $\chi^{(2)}_{eoo}$, $\chi^{(2)}_{ooo}$,
$\chi^{(2)}_{ooe}$.[23] The experimental axes were defined to overlap with the
LN principal crystallographic coordinate. Optical axis of the LN is nominated
as $z(e)$-axis, and the light propagated along $-x$-direction, as shown in
Fig. 1(b).

Figure 3:  Polarization images of ‘L’ and ‘N’. (a) Images of Stokes parameters
$s_{1}$, $s_{2}$, and $s_{3}$. The first, second and third columns give the
experimentally measured results of SF waves, the retrieved FF signal, and the
directly measured FF results, respectively. Outlines of the letters are
sketched by dashed lines. Scale bar of 100 $\mu$m is given in bottom-left
conner of the first picture. (b) Images are given in forms of azimuth ($\phi$)
and ellipticity ($\xi$) of the polarization ellipse. It is clear that the
retrieved FF polarization parameters reproduce the directly measured FF
results very well, confirming the precision of our technique.

A schematic of the experimental setup is shown in Fig. 1(b). An optical pump-
probe system was built. The near-infrared (NIR) laser pulses at a wavelength
of 808 nm were generated by a Kerr lens mode-locked Ti:sapphire laser
oscillator (Spectra Physics, Maitai, pulse width of 230 fs, repetition rate of
80 MHz). The pulse train was divided into two paths using a beam splitter. A
delay line, consisting of a linear translation stage in the pump beam path,
was used to compensate the length difference between the pump and probe paths
and ensured the fine control of the pulse synchronization. The pulses were
adjusted to be well overlapped temporally to guarantee the strongest SFG. The
optical beams were focused onto the LN film afterwards by lens pairs with
focal lengths of 500 mm and 100 mm, forming focal points with the diameter of
135 $\mu$m. The intensity of the pump and probe beams were about 142 MW/cm2
and 130 MW/cm2 on the LN surface, respectively. The polarization of the pump
beam was fixed along the ordinary ($o$-) axis of the LN crystal. The generated
SF light at 404 nm was collected in the forward direction by a lens with the
focal length of 75 mm. Short-pass filters (BG40, not shown) were inserted in
the SF emission path to block the remaining 808 nm light. The SF light then
entered a home-built full-Stokes polarization state analyzer (PSA) comprising
a rotating quarter-wave plate (QWP), a Glan-Taylor (GT) calcite polarizer, and
afterwards was detected by a monochrome electron-multiplying charge-coupled
device (EMCCD, iXon Ultra 888). 91 frames of SF signal were recorded while the
QWP was rotated from 0 to 180∘ with a discrete step of 2∘. Through Fourier
transformation of the recorded SF intensity as a function of the QWP rotation
angle, the Stokes parameters of the SF wave could be obtained.[24]

##  III Calibration of IRP

In order to accurately reconstruct $\mathbf{S}^{\omega_{1}}$ of the incident
IR signal light via the measured SF Stokes parameters according to Eq. (3),
the first objective is to determine the Mueller matrix
$\mathbf{M}\left(\chi^{(2)},\mathbf{S}^{\omega_{2}}\right)$ of the IRP system.
Such process could be referred as calibration of the IRP, which could be
accomplished by resorting to multiple sets of known input and output
polarization pairs.

The polarization states of the signal beam incident on the LN were prepared by
a polarization state generator (PSG) consisting of a GT polarizer, a half-wave
plate (HWP) and a QWP. For the calibration end, 72 different signal
polarization states that uniformly sample the Poincaré sphere were chosen as
input signal polarizations, whose $\mathbf{S}^{\omega_{1}}$ parameters were
further experimentally characterized using the PSA [as illustrated by red dots
in Fig. 2(a)]. The corresponding generated $\mathbf{S}^{\omega_{3}}$ were
measured as well and labeled by blue dots on the SF Poincaré sphere in Fig.
2(b). The 16 elements of $\mathbf{M}$ linking $\mathbf{S}^{\omega_{1}}$ to
$\mathbf{S}^{\omega_{3}}$ can be numerically determined using a least squares
method. Afterwards, the FF polarization states can be retrieved through Eq.
(3). For the clearer comparison between the original and the retrieved IR
signal polarizations, the three-dimensional Poincaré spheres are expanded into
planar maps, as shown in Fig. 2(c,d), in forms of angles of azimuth
[$\phi=\frac{1}{2}\arctan(\frac{s_{2}}{s_{1}}$)] and ellipticity
[$\xi=\frac{1}{2}\arctan(\frac{s_{2}}{\sqrt{s_{1}^{2}+s_{2}^{2}}}$)], which
correspond to the physical orientation and the shape of the polarization
ellipse, respectively. The well overlap between the reconstructed polarization
states (small gray dots) and the directly measured ones (large color dots) in
Fig. 2(c) demonstrates the excellent precision of our IRP in reconstructing FF
signal polarization states from the SF polarizations.

Figure 4:  Polarization images of a ‘house’ mask with a non-uniform
birefringence distribution. The images are given in forms of Stokes parameters
($s$), azimuth ($\phi$) and ellipticity ($\xi$). The retrieved FF signal
polarization images (the second column) reproduce the measured FF values (the
third column) very well.

##  IV Stokes polarimetric imaging

We now apply our IRP in full-Stokes polarization imaging. Such polarization
imaging technique aims to map the polarization states across a scene of
interest. In the traditional IR imaginary polarimetry, the polarization
distribution over the cross-section of the signal beam is directly
characterized by the PSA equipped by a IR camera. However, here we would show
the retrieval of the spatial polarization information contained in the IR
signal beam via analyzing that of the SF wave. To demonstrate such
application, binary masks with transparent letters of ‘L’ and ‘N’ were used as
the targets, which were fabricated by focused ion beam milling through a 200
nm thick opaque gold film supported by a SiO2 substrate [the intensity images
of the masks are given in Fig. 1(b)]. The masks were inserted into the signal
beam, and were projected onto the LN film plane. The polarization states of
the signal light illuminating the mask were adjusted by the PSG in front. In
this configuration, the polarization states across the signal beam were
uniformly distributed. After interacting with the pump beam, both the
morphological and polarization information of the masks contained in the
signal beam were successfully transferred into the generated SF light and
recorded by the PSA. The SF Stokes parameters were analyzed at each pixel, and
the polarization images at the SF wavelength were constructed afterwards.
Figure 3 presents the results of the masks by our IRP system. In the
literatures of polarization imaging, different sets of polarization parameters
are commonly used. The first one is the Stokes parameters ($s$). In contrast
to the traditional intensity photography that gives the same images despite
the signal light polarization varies, the $s_{1}$, $s_{2}$ and $s_{3}$ present
a good capability in distinguishing different signal polarizations. The Stokes
parameters are uniformly distributed inside the letter area, while the
background is noisy because no reliable SF signal was detected resulting from
the opacity outside the letter areas. For each mask, the Stokes images at the
SF wavelength are given in the first column. The retrieved FF images following
Eq. (3) are given in the second column, and the directly measured FF results
are given in the third column. It is clear that the retrieved FF polarizations
reproduce the directly measured FF values very well, indicating the fidelity
and precision of our technique in reconstructing the FF polarization images.
Images of $\phi$ and $\xi$ are further plotted in Fig. 3(b) to intuitively
present the geometric features of the polarization ellipse, which also confirm
the good consistence between the retrieved and the directly measured FF signal
values.

We further implement our IRP technique to the target that contains an
inhomogeneous birefringent structure, which causes multiple polarization
components in the signal light. Such scenario holds a high application
significance in the fields of remote sensing for mineral exploration, as well
as microscopy for characterizing the anisotropic biological tissues. To
demonstrate such functionality, a commercial depolarizer was adopted as the
target. A ‘house’ shaped mask was adopted [as shown by the intensity image in
Fig. 1(b)] and put on the surface of the depolarizer. Such depolarizer is a
patterned microretarder array, which consists of a thin film of birefringent
liquid crystal polymer sandwiched between two glass plates. The fast-axis
orientation and phase retardation of the liquid crystal were designed to vary
with a periodical pattern inside the depolarizer. Thus, after transmitting
through such depolarizer, a spatially varied polarization distribution would
be produced over the cross section of the signal light. It is clear that the
birefringence feature of the depolarizer is not discernible by the intensity
image in Fig. 1(b), however, it is very well resolved in the polarization
images in Fig. 4. The pseudo-color plot essentially traces the contour of the
spatial birefringence variation of the liquid crystal, and each band of color
in the polarization images represents a region with the same polarization
state. Similar as Fig. 3, the retrieved FF polarization images here (the
second column of Fig. 4) reproduce the directly measured FF values (the third
column) very well, confirming the precision of our technique again.

##  V Conclusion

In conclusion, we build a nonlinear up-conversion full-Stokes IRP that manages
to retrieve the polarization information of the incident IR probe light from
that of the up-converted SF wave. The IR polarization images of targets with
either uniform or non-uniform polarization distributions are successfully
reconstructed using our IRP system with high precision. Despite we focus here
on the degenerate SFG process that both the pump and probe light have the same
wavelength at NIR, it is expected that our method also validates for the non-
degenerate SFG process where the signal photon can be in mid- or long-infrared
range and up-converted by visible pump light.[25, 26] Thus it is applicable to
vibrational spectroscopy for birefringent molecule identification, and thermal
imaging for night vision detection. Furthermore, utilizing the current state-
of-the-art scheme of nonlinear enhancement by nanostructures such as
metasurfaces[27, 28, 29, 30] or intracavity up-conversion configuration,[31,
32] the efficiency of the SFG process from the LN film could be furthermore
effectively improved. We could foresee a direct impact of our results on a
variety of IR polarimetric applications, such as optical crystallography,
chemical sensing, disease diagnosis, and explosive detection, etc.

###### Acknowledgements.

Z.Z. and D.Z. contribute equally to this work. This work was supported by
National Key R&D Program of China (2017YFA0305100, 2017YFA0303800,
2019YFA0705000); National Natural Science Foundation of China (92050114,
91750204, 61775106, 11904182, 12074200, 11774185); Guangdong Major Project of
Basic and Applied Basic Research (2020B0301030009); 111 Project (B07013);
PCSIRT (IRT0149); Open Research Program of Key Laboratory of 3D Micro/Nano
Fabrication and Characterization of Zhejiang Province; Fundamental Research
Funds for the Central Universities (010-63201003, 010-63201008, 010-63201009,
010-63211001); Tianjin Youth Talent Support Program. We thank Nanofabrication
Platform of Nankai University for fabricating samples.

## References

## References

  * [1] B. H. Stuart, _Infrared spectroscopy: fundamentals and applications_ (John Wiley & Sons, 2004).
  * [2] D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. G. De Abajo, V. Pruneri, and H. Altug, “Mid-infrared plasmonic biosensing with graphene,” Science 349, 165–168 (2015).
  * [3] B. Kloppenborg, R. Stencel, J. D. Monnier, G. Schaefer, M. Zhao, F. Baron, H. McAlister, T. Ten Brummelaar, X. Che, C. Farrington _et al._ , “Infrared images of the transiting disk in the $\varepsilon$ aurigae system,” Nature 464, 870–872 (2010).
  * [4] K. Pfeilsticker, A. Lotter, C. Peters, and H. Bösch, “Atmospheric detection of water dimers via near-infrared absorption,” Science 300, 2078–2080 (2003).
  * [5] R. J. Keyes, _Optical and infrared detectors_ , vol. 19 (Springer Science & Business Media, 2013).
  * [6] A. Rogalski, “Infrared detectors: an overview,” Infrared Phys. Technol. 43, 187–210 (2002).
  * [7] A. Barh, P. J. Rodrigo, L. Meng, C. Pedersen, and P. Tidemand-Lichtenberg, “Parametric upconversion imaging and its applications,” Adv. Opt. Photon. 11, 952–1019 (2019).
  * [8] Y.-R. Shen, _The principles of nonlinear optics_ (Wiley-Interscience, 1984).
  * [9] J. Midwinter and J. Warner, “Up-conversion of near infrared to visible radiation in lithium-meta-niobate,” J. Appl. Phys. 38, 519–523 (1967).
  * [10] J. E. Midwinter, “Image conversion from 1.6 $\mu$ to the visible in lithium niobate,” Appl. Phys. Lett. 12, 68–70 (1968).
  * [11] Z.-Y. Zhou, S.-L. Liu, Y. Li, D.-S. Ding, W. Zhang, S. Shi, M.-X. Dong, B.-S. Shi, and G.-C. Guo, “Orbital angular momentum-entanglement frequency transducer,” Phys. Rev. Lett. 117, 103601 (2016).
  * [12] W. S. Boyle and G. E. Smith, “Charge coupled semiconductor devices,” Bell Sys. Tech. J. 49, 587–593 (1970).
  * [13] A. Hugi, G. Villares, S. Blaser, H. Liu, and J. Faist, “Mid-infrared frequency comb based on a quantum cascade laser,” Nature 492, 229–233 (2012).
  * [14] A. Schliesser, N. Picqué, and T. W. Hänsch, “Mid-infrared frequency combs,” Nature Photon. 6, 440–449 (2012).
  * [15] A. Torregrosa, H. Maestre, and J. Capmany, “Intra-cavity upconversion to 631 nm of images illuminated by an eye-safe ase source at 1550 nm,” Opt. Lett. 40, 5315–5318 (2015).
  * [16] S. Junaid, S. C. Kumar, M. Mathez, M. Hermes, N. Stone, N. Shepherd, M. Ebrahim-Zadeh, P. Tidemand-Lichtenberg, and C. Pedersen, “Video-rate, mid-infrared hyperspectral upconversion imaging,” Optica 6, 702–708 (2019).
  * [17] J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45, 5453–5469 (2006).
  * [18] Y. P. Svirko and N. I. Zheludev, _Polarization of light in nonlinear optics_ (Wiley-VCH, 2000).
  * [19] N. A. Rubin, G. D’ Aversa, P. Chevalier, Z. Shi, W. T. Chen, and F. Capasso, “Matrix fourier optics enables a compact full-Stokes polarization camera,” Science 365, eaax1839 (2019).
  * [20] R. W. Schmieder, “Stokes-algebra formalism,” J. Opt. Soc. A 59, 297–302 (1969).
  * [21] M. Bass, C. DeCusatis, J. Enoch, V. Lakshminarayanan, G. Li, C. Macdonald, V. Mahajan, and E. Van Stryland, _Handbook of Optics, Third Edition Volume I: Geometrical and Physical Optics, Polarized Light, Components and Instruments(Set)_ (McGraw-Hill, Inc., New York, 2010), 3rd ed.
  * [22] R. W. Boyd, _Nonlinear Optics_ (Elsevier, 2003).
  * [23] J. Ma, J. Chen, M. Ren, W. Wu, W. Cai, and J. Xu, “Second-harmonic generation and its nonlinear depolarization from lithium niobate thin films,” Opt. Lett. 45, 145–148 (2020).
  * [24] M. Ren, M. Chen, W. Wu, L. Zhang, J. Liu, B. Pi, X. Zhang, Q. Li, S. Fan, and J. Xu, “Linearly polarized light emission from quantum dots with plasmonic nanoantenna arrays,” Nano Lett. 15, 2951–2957 (2015).
  * [25] J. Falk and J. Yarborough, “Detection of room-temperature blackbody radiation by parametric upconversion,” Appl. Phys. Lett. 19, 68–70 (1971).
  * [26] I. Kviatkovsky, H. M. Chrzanowski, E. G. Avery, H. Bartolomaeus, and S. Ramelow, “Microscopy with undetected photons in the mid-infrared,” Sci. Adv. 6, eabd0264 (2020).
  * [27] G. Li, S. Zhang, and T. Zentgraf, “Nonlinear photonic metasurfaces,” Nature Rev. Mater. 2, 17010 (2017).
  * [28] J. Ma, F. Xie, W. Chen, J. Chen, W. Wu, W. Liu, Y. Chen, W. Cai, M. Ren, and J. Xu, “Nonlinear lithium niobate metasurfaces for second harmonic generation,” Laser Photon. Rev. 15, 2000521 (2021).
  * [29] A. Fedotova, M. Younesi, J. Sautter, A. Vaskin, F. J. Löchner, M. Steinert, R. Geiss, T. Pertsch, I. Staude, and F. Setzpfandt, “Second-harmonic generation in resonant nonlinear metasurfaces based on lithium niobate,” Nano Lett. 20, 8608–8614 (2020).
  * [30] L. Carletti, A. Zilli, F. Moia, A. Toma, M. Finazzi, C. De Angelis, D. N. Neshev, and M. Celebrano, “Steering and encoding the polarization of the second harmonic in the visible with a monolithic LiNbO3 metasurface,” ACS Photon. 8, 731–737 (2021).
  * [31] J. S. Dam, P. Tidemand-Lichtenberg, and C. Pedersen, “Room-temperature mid-infrared single-photon spectral imaging,” Nature Photon. 6, 788–793 (2012).
  * [32] N. M. Israelsen, C. R. Petersen, A. Barh, D. Jain, M. Jensen, G. Hannesschläger, P. Tidemand-Lichtenberg, C. Pedersen, A. Podoleanu, and O. Bang, “Real-time high-resolution mid-infrared optical coherence tomography,” Light: Sci. Appl. 8, 11 (2019).


11institutetext: Harvard-MIT Health Sciences and Technology, Massachusetts
Institute of Technology, Cambridge, MA, USA
22institutetext: Computer Science and Artificial Intelligence Laboratory,
Massachusetts Institute of Technology, Cambridge, MA, USA
22email<EMAIL_ADDRESS>
# Fast Auto-Differentiable Digitally Reconstructed Radiographs for Solving
Inverse Problems in Intraoperative Imaging

Vivek Gopalakrishnan 1122 Polina Golland 1122

###### Abstract

The use of digitally reconstructed radiographs (DRRs) to solve inverse
problems such as slice-to-volume registration and 3D reconstruction is well-
studied in preoperative settings. In intraoperative imaging, the utility of
DRRs is limited by the challenges in generating them in real-time and
supporting optimization procedures that rely on repeated DRR synthesis. While
immense progress has been made in accelerating the generation of DRRs through
algorithmic refinements and GPU implementations, DRR-based optimization
remains slow because most DRR generators do not offer a straightforward way to
obtain gradients with respect to the imaging parameters. To make DRRs
interoperable with gradient-based optimization and deep learning frameworks,
we have reformulated Siddon’s method, the most popular ray-tracing algorithm
used in DRR generation, as a series of vectorized tensor operations. We
implemented this vectorized version of Siddon’s method in PyTorch, taking
advantage of the library’s strong automatic differentiation engine to make
this DRR generator fully differentiable with respect to its parameters.
Additionally, using GPU-accelerated tensor computation enables our vectorized
implementation to achieve rendering speeds equivalent to state-of-the-art DRR
generators implemented in CUDA and C`++`. We illustrate the resulting method
in the context of slice-to-volume registration. Moreover, our simulations
suggest that the loss landscapes for the slice-to-volume registration problem
are convex in the neighborhood of the optimal solution, and gradient-based
registration promises a much faster solution than prevailing gradient-free
optimization strategies. The proposed DRR generator enables fast computer
vision algorithms to support image guidance in minimally invasive procedures.
Our implementation is publically available at https://github.com/v715/DiffDRR.

###### Keywords:

DRRs Differentiable programming Inverse problems.

##  1 Introduction

Digitally reconstructed radiographs (DRRs) are simulated 2D X-ray images
generated from 3D computational tomography (CT) volumes using a variety of
ray-tracing techniques. While DRRs are widely used in preoperative settings
(e.g., optimizing dose delivery in radiation oncology), many potentially
valuable intraoperative use cases (e.g., real-time multimodal registration for
image-guided procedures) are infeasible due to computational bottlenecks in
generating DRRs and using them in slice-to-volume registration. Most open-
source CPU-based implementations take about $1\text{\,}\sec$ to generate a
single DRR [9], which is not fast enough for intraoperative imaging systems
with sampling rates of about 7.5 frames per second [8]. Numerous GPU-
accelerated DRR generators have been proposed with run times on the order of
$100\text{\,}\mathrm{m}\sec$ [1, 5, 7], but to the best of our knowledge, no
publically available implementation exists.

The second limitation of currently available DRR generators is that it is
challenging to efficiently compute derivatives using these renderers because
they are implemented in low-level languages such as CUDA and C`++`. DRR
generators are often used in conjunction with numerical optimization schemes
to solve fundamental medical imaging problems (e.g., slice-to-volume
registration), and the difficulty in computing derivatives means that
gradient-based optimization techniques are often infeasible [13]. While many
end-to-end deep learning approaches can solve X-ray to CT registration
problems with high accuracy [2, 3], these methods often require large amounts
of training data, which can make them impractical for specialized
interventional problems. Instead, many applications use iterative gradient-
free methods, such as the Nelder-Mead method [11], to optimize an image
similarity metric with respect to the parameters of the DRR generator [3, 13].
While these methods are effective for optimizing highly nonlinear loss
functions, we show that the loss landscapes for slice-to-volume registration
in particular is convex in a large region around the optimum, making this
problem better suited for gradient-based optimization methods.

We present a fast vectorized renderer that generates DRRs and their
derivatives with respect to image geometry parameters automatically. We
utilize PyTorch as a GPU-accelerated tensor algebra library with robust
source-to-source automatic differentiation to implement differentiable DRRs.
That is, using our implementation, DRR generation can be used as a
differential operator to train deep learning algorithms for fast
reconstruction and registration algorithms. We analyze the performance of our
implementation and the correctness of the automatically obtained derivatives,
and demonstrate an experiment where our differentiable DRR generator solves a
slice-to-volume registration problem. Our hope is that this open-source
package will be useful for translating computer vision algorithms to real-time
implementations for interventional applications.

##  2 Methods

We start by summarizing Siddon’s method [10], commonly used for ray-tracing in
DRR synthesis, and its extensions that accelerate rendering speed. We then
describe our vectorized implementation of Siddon’s method, which achieves
rendering speeds equivalent to those of existing GPU-accelerated methods while
also being fully differentiable.

Figure 1: DRR synthesis. (a) We assume an idealized model of a projectional
radiography imaging system: X-ray beams are emitted with a fixed initial
energy from a point source $\mathbf{s}\in\mathbb{R}^{3}$, beam energy
diminishes as the X-rays travel through the CT volume $\mathbf{V}$, and energy
in the attenuated beams is measured when the X-rays hit a point on the
detector $\mathbf{p}\in\mathbb{R}^{3}$, producing the DRR. (b) In Siddon’s
method, the image location value at $\mathbf{p}$ is a weighted average of the
intensities of the voxels through which the ray passes, where the weight is
the length of the ray’s intersection with the voxel. The values $\alpha_{m}$
and $\alpha_{m+1}$ parameterize the intersection of the ray with two adjacent
planes, and the midpoint $\frac{\alpha_{m+1}+\alpha_{m}}{2}$ identifies the
current voxel through which the ray is passing. (c) Our vectorized Siddon’s
method generates a $200\times 200$ DRR in $72.7\text{\,}\mathrm{m}\sec$ $\pm$
$10\text{\,}\mathrm{\SIUnitSymbolMicro}\sec$ on an NVIDIA GeForce RTX 2080 Ti.

###  2.1 DRR Generation

The process of generating a DRR models the geometry of an idealized
projectional radiography system (Fig. 1a). Let $\mathbf{s}\in\mathbb{R}^{3}$
be the X-ray source, and $\mathbf{p}\in\mathbb{R}^{3}$ be a target pixel on
the detector plane. Then $R(\alpha)=\mathbf{s}+\alpha(\mathbf{p}-\mathbf{s})$
is a ray that originates from $\mathbf{s}$ ($\alpha=0$), passes through the
imaged volume, and hits the detector plane at $\mathbf{p}$ ($\alpha=1$). The
total energy attenuation experienced by the X-ray by the time it reaches pixel
$\mathbf{p}$ is given by the following line integral:

$E(R)=\|\mathbf{p}-\mathbf{s}\|_{2}\int_{0}^{1}\mathbf{V}\left(\mathbf{s}+\alpha(\mathbf{p}-\mathbf{s})\right)\mathrm{d}\alpha,$
(1)

where $\mathbf{V}:\mathbb{R}^{3}\mapsto\mathbb{R}$ is the imaged volume. The
term $\|\mathbf{p}-\mathbf{s}\|_{2}$ endows the unit-free $\mathrm{d}\alpha$
with the physical unit of length. For DRR synthesis, $\mathbf{V}$ is
approximated by a discrete 3D CT volume, and Eq. (1) becomes

$E(R)=\|\mathbf{p}-\mathbf{s}\|_{2}\sum_{m=1}^{M-1}(\alpha_{m+1}-\alpha_{m})\mathbf{V}\left[\mathbf{s}+\frac{\alpha_{m+1}+\alpha_{m}}{2}(\mathbf{p}-\mathbf{s})\right],$
(2)

where $\alpha_{m}$ parameterizes the locations where ray $R$ intersects one of
the orthogonal planes comprising the CT volume, and $M$ is the number of such
intersections (Fig. 1b). Note that this model does not account for patterns of
reflection and scattering that are present in real X-ray systems. While these
simplifications preclude synthesis of realistic X-rays, the model in Eq. (2)
has been widely and successfully used in slice-to-volume registration [13].
Additionally, our approach of vectorizing DRR generation might also be
interoperable with more sophisticated image synthesis models, an extension we
examine further in the Discussion.

###  2.2 Siddon’s Method and Its GPU Extensions

Siddon’s method [10] provides a parametric method to identify the plane
intersections $\\{\alpha_{m}\\}_{m=1}^{M}$. Let $\Delta X$ be the CT voxel
size in the $x$-direction and $b_{x}$ be the location of the $0$-th plane in
this direction. Then the intersection of ray $R$ with the $i$-th plane in the
$x$-direction is given by

$\alpha_{x}(i)=\frac{b_{x}+i\Delta
X-\mathbf{s}_{x}}{\mathbf{p}_{x}-\mathbf{s}_{x}},$ (3)

with analogous expressions for $\alpha_{y}(\cdot)$ and $\alpha_{z}(\cdot)$. We
can use Eq. (3) to compute the values $\bm{\alpha}_{x}$ for all the
intersections between $R$ and the planes in the $x$-direction:

$\bm{\alpha}_{x}=\\{\alpha_{x}(i_{\min}),\dots,\alpha_{x}(i_{\max})\\},$

where $i_{\min}$ and $i_{\max}$ denote the first and last intersections of $R$
with the $x$-direction planes. Defining $\bm{\alpha}_{y}$ and
$\bm{\alpha}_{z}$ analogously, we construct the array

$\bm{\alpha}=\mathrm{sort}(\bm{\alpha}_{x},\bm{\alpha}_{y},\bm{\alpha}_{z}),$
(4)

which contains $M$ values of $\alpha$ parameterizing the intersections between
$R$ and the orthogonal $x$-, $y$-, and $z$-directional planes. We substitute
values in the sorted set $\bm{\alpha}$ into Eq. (2) to evaluate $E(R)$, which
corresponds to the intensity of pixel $\mathbf{p}$ in the synthesized DRR.

A faster variant determines consecutive intersecting planes iteratively [4].
For example, the value of $\alpha$ at the second plane intersected by $R$ is
given by
$\alpha_{2}=\min\\{\alpha_{x}(i_{\min}+1),\alpha_{y}(j_{\min}+1),\alpha_{z}(k_{\min}+1)\\}$.
The algorithm iteratively finds the next value of $\alpha$ until we reach the
edge of the CT volume, making this approach more memory efficient by requiring
fewer intermediate values to be stored. This modified algorithm, known as
Siddon-Jacobs’ method, is commonly implemented in CUDA and C`++` to create
multi-threaded GPU-accelerated DRR generators that exploit data parallelism by
assigning each thread to trace an independent ray intersecting the detector
plane [1, 5, 7].

###  2.3 Vectorizing Siddon’s Method

While Siddon-Jacobs’ method is more memory efficient, the iterative loop it
relies on is not amenable to implementations in vectorized tensor algebra
libraries. Thus we vectorize the original Siddon’s method as follows. Let
$\mathbf{P}\in\mathbb{R}^{H\times W\times 3}$ contain the 3D pixel coordinates
of a DRR with dimension $H\times W$. We compute the $\alpha$ values for
intersections with all of the $x$-, $y$-, and $z$-planes for all
$\mathbf{p}\in\mathbf{P}$ in parallel:

$\mathbf{A}=\left(\begin{pmatrix}b_{x}\\\ b_{y}\\\
b_{z}\end{pmatrix}+\begin{pmatrix}i\\\ j\\\
k\end{pmatrix}\otimes\begin{pmatrix}\Delta X\\\ \Delta Y\\\ \Delta
Z\end{pmatrix}-\mathbf{s}\right)\oslash(\mathbf{P}-\mathbf{s})\in\mathbb{R}^{H\times
W\times(n_{x}+n_{y}+n_{z})},$ (5)

where $(n_{x},n_{y},n_{z})$ are the dimensions of the CT volume $\mathbf{V}$,
$(i,j,k)$ are the CT voxel indices, $(\Delta X,\Delta Y,\Delta Z)$ are the CT
voxel sizes, and $\otimes$ and $\oslash$ are the Hadamard product and division
operators, respectively. Rather than explicitly compute the indices
$(i_{\min},i_{\max})$, $(j_{\min},j_{\max})$, and $(k_{\min},k_{\max})$ for
each ray, as is done in Siddon’s original method, we instead compute the
minimum and maximum values of $\bm{\alpha}$, corresponding to when each ray
enters and exits the volume:

$\displaystyle\bm{\alpha}_{\min}$
$\displaystyle=\max\big{\\{}\min\\{\bm{\alpha}_{x}(0),\bm{\alpha}_{x}(n_{x})\\},\min\\{\bm{\alpha}_{y}(0),\bm{\alpha}_{y}(n_{y})\\},\min\\{\bm{\alpha}_{z}(0),\bm{\alpha}_{z}(n_{z})\\}\big{\\}}$
$\displaystyle\bm{\alpha}_{\max}$
$\displaystyle=\min\big{\\{}\max\\{\bm{\alpha}_{x}(0),\bm{\alpha}_{x}(n_{x})\\},\max\\{\bm{\alpha}_{y}(0),\bm{\alpha}_{y}(n_{y})\\},\max\\{\bm{\alpha}_{z}(0),\bm{\alpha}_{z}(n_{z})\\}\big{\\}},$

where $\bm{\alpha}_{\min},\bm{\alpha}_{\max}\in\mathbb{R}^{H\times W}$. We
filter $\mathbf{A}$ to include only values in the range
$[\bm{\alpha}_{\min},\bm{\alpha}_{\max}]$ and sort each row
$\mathbf{A}(h,w,\cdot)$ for $h\in\\{1,\dots,H\\},w\in\\{1,\dots,W\\}$.
Finally, we evaluate Eq. (2) with this sorted tensor to compute the intensity
for each pixel in the DRR, completing a chain of vectorized tensor operations.

Because we reformulated the original Siddon’s method as a series of tensor
operations, our vectorized version benefits from the mature GPU compilers and
memory allocators developed for optimizing large-scale deep learning models.
For empirical evaluation of our method, we also implemented a partially-
vectorized version of Siddon-Jacobs’ method in which the $\alpha$ updates are
still computed iteratively (i.e., with a loop), but the updates are applied in
a vectorized form to every target pixel in the detector plane.

###  2.4 Differentiating DRRs with Respect to Imaging Parameters

We specify the 3D position of the X-ray source and detector plane relative to
the CT volume with the following seven geometric parameters: radius $\rho$
that acts as a scaling factor; three rotational degrees of freedom (DoF)
$(\theta,\varphi,\gamma)$; and three translational DoF $(b_{x},b_{y},b_{z})$.
Using spherical coordinates, we express the position of the X-ray source as
$\mathbf{s}=(\rho,\theta,\varphi)$, where $\rho$ is half of the source-to-
detector distance, and $\theta$ and $\varphi$ are the azimuthal and polar
angles, respectively. We assume the detector plane is tangent to this implied
sphere at the point opposite $\mathbf{s}$. The orientation of this plane is
determined by a rotation about the $x$-axis by the angle $\gamma$. We add the
translation $(b_{x},b_{y},b_{z})$ to the coordinates of the X-ray source and
detector plane to create a reference frame wherein the patient is not
perfectly centered relative to the X-ray scanner.

Since every step of our pipeline, from the generation of the pixels on the
detector plane to the computation of the pixel intensities, is performed in
PyTorch’s tensor framework, the resulting DRRs are differentiable with respect
to the parameters described above. That is, the gradient
$\nabla_{\bm{\eta}}\mathcal{L}(I(\bm{\eta}))$ of the loss function
$\mathcal{L}(\cdot)$ evaluated for the DRR $I(\bm{\eta})$ with respect to the
parameters $\bm{\eta}=(\rho,\theta,\varphi,\gamma,b_{x},b_{y},b_{z})$ is
obtained automatically for any differentiable $\mathcal{L}(\cdot)$.

##  3 Experiments

We evaluate the proposed efficient implementation of the algorithm on a
reference chest CT scan from Slicer3D [6]. We compare the performance of the
method to two baseline approaches and illustrate its application for gradient-
based slice-to-volume registration.

###  3.1 Performance Analysis

We compare our vectorized GPU version of Siddon’s method (VGS) to two baseline
approaches: a widely used CPU implementation in the Plastimatch package (CP)
[9], and our vectorized GPU implementation of Siddon-Jacobs’ method (VGSJ)
which we described in Section 2.3. Note that the CUDA-accelerated DRR
generator in Plastimatch is not working at the time of publication. GPU
benchmarks are run on an NVIDIA GeForce RTX 2080 Ti, and CPU benchmarks are
run on an 18-core Linux computer with Intel(R) Xeon(R) CPU E5-2697 v4 @
2.30GHz processors.

####  3.1.1 Results.

Table 1 summarizes statistics for the run time and accuracy of our method and
the two baseline approaches, as well as intensity differences between our
method and Plastimatch, and gradient differences between our method and FFD.
Our implementation is much faster than Plastimatch, which is to be expected as
Plastimatch is executed on the CPU. Numerically, the two implementations are
very similar with an average root-mean-square error (RMSE) of (8.3$\pm$1.9)
$\times 10^{-4}$ where the images are normalized to the range of $[0,1]$. For
DRRs smaller than $H=500$, our method is faster than the vectorized version of
Siddon-Jacobs’ despite our method’s high memory requirements. However, at
$H=500$, they have roughly equivalent run times. For smaller image sizes, our
DRR generator achieves equivalent run times to previously reported GPU-
accelerated implementations [1, 5, 7]. The gradients obtained via PyTorch
auto-differentiation for our method are within 0.05$\pm$0.01 of those computed
via forward finite differences with a step size of $10^{-6}$ and are an order
of magnitude more efficient to generate (35.1 ms ± 73.3 µs vs 400.5 ms ± 821.4
µs).

Table 1: Benchmark results. The dimension of the DRRs is $H\times W$. Each metric is averaged over $20$ runs. (VGS = Vectorized GPU Siddon’s method, VGSJ = Vectorized GPU Siddon-Jacobs’ method, and CP = CPU Plastimatch.) | Timing (ms) | RMSE | Autograd vs FFD
---|---|---|---
$H=W$ | VGS | VGSJ | CP | VGS vs CP | VGS
100 | 17.6 $\pm$ 0.05 | 380 $\pm$ 20.9 | 1028 $\pm$ 132 | (6.9$\pm$2.2)$\times 10^{-4}$ | 0.03 ± 0.02
200 | 72.7 $\pm$ 0.01 | 424 $\pm$ 4.2 | 1784 $\pm$ 488 | (8.7$\pm$2.7)$\times 10^{-4}$ | 0.06 ± 0.01
300 | 165 $\pm$ 0.13 | 432 $\pm$ 19.2 | 2941 $\pm$ 821 | (6.4$\pm$1.2)$\times 10^{-4}$ | 0.08 ± 0.02
400 | 296 $\pm$ 0.06 | 425 $\pm$ 2.8 | 6472 $\pm$ 643 | (9.0$\pm$1.9)$\times 10^{-4}$ | 0.03 ± 0.005
500 | 453 $\pm$ 41.2 | 425 $\pm$ 4.6 | 8472 $\pm$ 478 | (11.7$\pm$4.3)$\times 10^{-4}$ | 0.07 ± 0.006

###  3.2 DRR-based Gradient Descent for Slice-to-Volume Registration

We use our auto-differentiable DRR generator to implement slice-to-volume
registration with synthetic DRRs. Specifically, we generate a fixed DRR from a
set of ground truth parameters
$\bm{\eta}^{*}=(\theta,\varphi,\gamma,b_{x},b_{y},b_{z})$, and generate a
second moving DRR from a set of random initial parameters $\bm{\eta}_{0}$. We
use basic gradient descent to minimize the negative Zero-Normalized Cross-
Correlation (ZNCC) between the fixed DRR and the moving DRR.

####  3.2.1 Image Similarity Metrics are Locally Convex.

First, we conduct a simulation study to show that the loss landscape generated
by negative ZNCC is convex in a neighborhood around $\bm{\eta}^{*}$. We
generate moving DRRs by sampling rotational and positional displacements. We
sample all parameters uniformly from ranges of $90^{\circ}$ for $\theta$ and
$\varphi$, $45^{\circ}$ for $\gamma$, and $30\text{\,}\mathrm{mm}$ for
$b_{x}$, $b_{y}$, and $b_{z}$ around the ground truth parameters
$\bm{\eta}^{*}$. Although our generator provides gradients with respect to
other model parameters like the source-to-detector distance ($2\rho$) and the
DRR’s dimensions and pixel spacing ($H,W,\Delta x,\Delta y$), we assume that
those parameters are fixed (e.g., provided in the DICOM header), and instead
focus our analysis on this 6 DoF registration problem. We observe that
negative ZNCC is locally convex (Fig. 2), suggesting that it would be an apt
loss function to optimize with gradient descent. We observed similar loss
landscapes for the $L^{2}$ norm.

Figure 2: Negative ZNCC is convex around the optimal DRR parameters.
Rotational and positional displacements were sampled uniformly from ranges of
$90^{\circ}$ for $\theta$ and $\varphi$, $45^{\circ}$ for $\gamma$, and
$30\text{\,}\mathrm{mm}$ for $b_{x}$, $b_{y}$, and $b_{z}$. Figure 3:
Differentiable DRRs can be used to perform slice-to-volume registration. We
generated a moving DRR from randomly initialized parameters and used gradient
descent to maximize similarity with a fixed DRR. Convergence was achieved for
745/1000 simulated DRRs in an average 66 iterations ($1.92\text{\,}\sec$).
Examples of the optimization process are visualized at initial, intermediate
(the 20th and 40th iterations), and final steps. DRRs for which convergence
fails to occur get stuck in a local minimum with low negative ZNCC.

####  3.2.2 Differentiable DRR Registration Converges Quickly.

Given a fixed DRR and a moving DRR, we optimized the parameters of the moving
DRR with a basic implementation of gradient descent. We used different update
rates for the rotational and translational parameters because they have
different units ($\beta_{\theta\varphi\gamma}=5.3\times 10^{-2}$ and
$\beta_{xyz}=7.5\times 10^{1}$), and momentum $\lambda=0.9$. Additionally, to
investigate the local minima in which our gradient descent algorithm could get
stuck, we expanded the space of possible initializations beyond a convex
neighborhood to $120^{\circ}$ for $\theta$, $\varphi$, and $\gamma$, and
$60\text{\,}\mathrm{mm}$ for $b_{x}$, $b_{y}$, and $b_{z}$. The size of the CT
volume is $360\text{\,}\mathrm{mm}$ $\times$ $360\text{\,}\mathrm{mm}$
$\times$ $332.5\text{\,}\mathrm{mm}$.

For each randomly initialized DRR, we ran 250 iterations of gradient descent
and determined that the moving DRR had converged to the fixed DRR if the
negative ZNCC between the two images was less than $-0.999$. If this did not
occur within 250 iterations, we treated the run as having failed to converge.
Of the 1,000 randomly initialized DRRs we generated, 745 converged and 255
failed to converge (Fig. 3a). The initializations that converged solved the 6
DoF slice-to-volume registration problem $65.48\pm 14.27$ iterations
($1.92\text{\,}$ $\pm$ $0.43\text{\,}\sec$).

We visualize multiple optimization steps for our gradient descent registration
algorithm (Fig. 3c). From the examples that converged, our model successfully
recovers the true pose parameters from challenging initializations, reaching a
more reasonable estimate by the intermediate 20th iteration (Fig. 3b). In one
example of an initialization that failed to converge, we see that the model
gets stuck in a local minimum (Fig. 3c). In the final iteration in this
example, the estimated DRR is orthogonal to the sagittal plane. The geometry
of this scene looks similar to the coronal DRR (ground truth), giving the
illusion of two lungs. We emphasize that our goal is not to propose a novel
registration algorithm, but to provide an efficient DRR synthesis procedure
that can support numerous downstream optimization applications, including
registration.

##  4 Discussion

We present a fast auto-differentiable DRR generator that can be used to solve
inverse problems in intraoperative imaging. By reformulating Siddon’s method
for ray-tracing through a CT volume as series of vectorized tensor operations,
we obtain gradients of image loss functions with respect to DRR generating
parameters while achieving rendering speeds equivalent to multi-threaded GPU-
accelerated generators written in low-level languages such as CUDA and C`++`.
This approach promises to enable fast solutions to DRR-based optimization
problems with gradient methods, a strategy which was previously infeasible due
to the inefficiency of generating gradients with finite differences. We
demonstrate the effectiveness of this approach by solving a 6 DoF slice-to-
volume registration problem using a locally convex image loss function.

Our future research on auto-differentiable DRRs will investigate their use in
a deep learning framework to pre-train for specific intraoperative imaging
tasks (e.g., spatiotemporal registration in interventional cardiology). Such
pre-training could yield even faster solutions to inverse problems by finding
better initializations or powering end-to-end models. One issue that might
limit the effectiveness of such pre-training is that DRRs generated using
Siddon’s method are not realistic because they do not model any form of
scattering. We will investigate fusing our method with DeepDRR [12], a deep
learning framework that estimates scattering effects using Monte Carlo
simulations, to produce DRRs that are simultaneously realistic, fast, and
differentiable.

Acknowledgements. The authors thank Ruby Liu and Allie Forman for helpful
feedback and support. This work was supported, in part, by NIH NIBIB
5T32EB1680, NIH NIBIB NAC P41EB015902, and NIH NINDS U19NS115388.

## References

  * [1] De Greef, M., Crezee, J., Van Eijk, J., Pool, R., Bel, A.: Accelerated ray tracing for radiotherapy dose calculations on a gpu. Medical physics 36(9Part1), 4095–4102 (2009)
  * [2] Esteban, J., Grimm, M., Unberath, M., Zahnd, G., Navab, N.: Towards fully automatic x-ray to ct registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 631–639. Springer (2019)
  * [3] Hou, B., Alansary, A., McDonagh, S., Davidson, A., Rutherford, M., Hajnal, J.V., Rueckert, D., Glocker, B., Kainz, B.: Predicting slice-to-volume transformation in presence of arbitrary subject motion. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 296–304. Springer (2017)
  * [4] Jacobs, F., Sundermann, E., De Sutter, B., Christiaens, M., Lemahieu, I.: A fast algorithm to calculate the exact radiological path through a pixel or voxel space. Journal of computing and information technology 6(1), 89–94 (1998)
  * [5] Mori, S., Kobayashi, M., Kumagai, M., Minohara, S.: Development of a gpu-based multithreaded software application to calculate digitally reconstructed radiographs for radiotherapy. Radiological physics and technology 2(1), 40–45 (2009)
  * [6] Pieper, S., Halle, M., Kikinis, R.: 3d slicer. In: 2004 2nd IEEE international symposium on biomedical imaging: nano to macro. pp. 632–635. IEEE (2004)
  * [7] Ruijters, D., ter Haar Romeny, B.M., Suetens, P.: Gpu-accelerated digitally reconstructed radiographs. BioMED 8, 431–435 (2008)
  * [8] Sadamatsu, K., Nakano, Y.: The effect of low frame rate fluoroscopy on the x-ray dose during coronary intervention. Internal medicine 55 15, 1943–6 (2016)
  * [9] Sharp, G.C., Li, R., Wolfgang, J., Chen, G., Peroni, M., Spadea, M.F., Mori, S., Zhang, J., Shackleford, J., Kandasamy, N.: Plastimatch: an open source software suite for radiotherapy image processing. In: Proceedings of the XVI’th International Conference on the use of Computers in Radiotherapy (ICCR) (2010)
  * [10] Siddon, R.L.: Fast calculation of the exact radiological path for a three-dimensional ct array. Medical physics 12(2), 252–255 (1985)
  * [11] Singer, S., Nelder, J.: Nelder-mead algorithm. Scholarpedia 4(7), 2928 (2009)
  * [12] Unberath, M., Zaech, J.N., Lee, S.C., Bier, B., Fotouhi, J., Armand, M., Navab, N.: Deepdrr–a catalyst for machine learning in fluoroscopy-guided procedures. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 98–106. Springer (2018)
  * [13] Van Der Bom, I., Klein, S., Staring, M., Homan, R., Bartels, L.W., Pluim, J.P.: Evaluation of optimization methods for intensity-based 2d-3d registration in x-ray guided interventions. In: Medical Imaging 2011: Image Processing. vol. 7962, pp. 657–671. SPIE (2011)