# Diabatic Hamiltonian matrix elements made simple

K. Hagino  Department of Physics, Kyoto University, Kyoto 606-8502, Japan
G.F. Bertsch  Department of Physics and Institute of Nuclear Theory, Box
351560, University of Washington, Seattle, Washington 98915, USA

###### Abstract

With a view to applying the Generator Coordinate Method to large configuration
spaces, we propose a simple approximate formula to compute diabatic many-body
matrix elements without having to evaluate two-body interaction matrix
elements. The method is illustrated with two analytically solvable
Hamiltonians based on the harmonic oscillator.

##  I Introduction

In nuclear physics we use transparent models not only to understand the
observed characteristics of nuclear spectra and dynamics but also to provide
checks and guidance on the complex computer codes that produce the
quantitative theory of their properties. The Lipkin model is a typical example
RS .

In this spirit we derive and apply a simple formula, Eq. (10) below, for
calculating matrix elements between configurations in the generator coordinate
method (GCM). In GCM, many Slater determinants are linearly superposed to
describe low-energy dynamics of many-body systems RS ; RG87 ; RRR18 . Since
those Slater determinants are not orthogonal, one needs to calculate their
overlaps and Hamiltonian matrix elements. To this end, the Hamiltonian matrix
elements are computed using the generalized Wick’s theorem (GWT)BB  as was
first done in Ref. bo90 .

However, the GWT requires a full knowledge of the Hamiltonian in a Fock-space
representation. In contrast, Eq. (10) below only needs the overlap matrix
elements of the orbitals and the orbital energies calculated from the single-
particle Hamiltonian. It is an approximate formula that can be used with any
of the theories in which an effective single-particle Hamiltonian can be
defined. This includes theories based on energy functionals such as Skyrme,
Gogny and the BCPMba13 , as well as those making use of the Strutinsky
treatment of the total energy. In this paper, we shall apply the formula to
two simple models, namely, the translational motion in a harmonic oscillator
potential and a simplified Nilsson model.

##  II Orbital energy formula

We first summarize the usual computational procedure. The starting point is a
configuration space of many-particle wave functions that have the form of
Slater determinants of orbitals in some common basis.111The extension to
include pairing by the Hartree-Fock-Bogoliubov equation is beyond the scope of
this paper. The Hartree-Fock method reduces a Hamiltonian with two-particle
interactions to a one-particle Hamiltonian with the orbitals as its
eigenstates.222 The method can easily be extended to Hamiltonians with three-
particle and higher terms by the same reduction of the higher-order density
matricesro12  to what has been done in the Hartree-Fock approximation. The
orbitals are obtained by minimizing the single-particle Hamiltonian or energy
functional. The space of configurations in the GCM is expanded by introducing
a collective coordinate. This is often implemented by adding a term in the
Hamiltonian to constrain the expectation value of some single-particle
operator, but the coordinate can also be defined in other ways. In the
examples below, it is convenient to introduce the coordinate directly into the
orbital wave functions.

Assuming that the constraint is described by a single parameter $q$, the task
is to calculate the Hamiltonian matrix element between a configuration
$|\Psi_{q}\rangle$ and some other configuration $|\Psi_{q^{\prime}}\rangle$.
The GWT method can be derived from Thouless’s theorem thouless1960  which is
used to express the second wave function the orbital basis of the first one.
That is,

$|\Psi_{q^{\prime}}\rangle={\cal
N}\prod_{h}\left(1+\sum_{p}c_{ph}a^{\dagger}_{p}a_{h}\right)|\Psi_{q}\rangle.$
(1)

Here $h$ and $p$ are occupied and unoccupied orbitals in the basis of
$|\Psi_{q}\rangle=\prod_{h}a^{\dagger}_{h}|\rangle$. The Fock-space operators
are expressed in the $q$ orbital basis and $\cal N$ is a normalization
constant given by

${\cal N}=(\det M)^{-1/2}$ (2)

where $M_{h,h^{\prime}}=\delta_{h,h^{\prime}}+\sum_{p}c_{ph}c_{ph^{\prime}}$.
With all of the operators expressed in the same basis it is easy to apply the
ordinary Wick’s theorem to calculate arbitrary many-body matrix elements.

To derive the orbital energy formula we place the reference configuration
midway between the two states $|\Psi_{q}\rangle$ and
$|\Psi_{q^{\prime}}\rangle$. Changing the notation somewhat, the many-particle
configurations can be written

$|\Psi_{\pm q}\rangle={\cal
N^{\pm}}\prod_{h}\left(1+\sum_{p}c_{ph}^{\pm}a^{\dagger}_{p}a_{h}\right)|\Psi_{0}\rangle,$
(3)

where $|\Psi_{0}\rangle$ is the reference configuration and $|\Psi_{\pm
q}\rangle$ are the end configurations. The overlap between the two end
configurations is given by

$\langle\Psi_{-q}|\Psi_{+q}\rangle={\cal N}^{-}{\cal
N}^{+}\det\left|\delta_{h,h^{\prime}}+\sum_{p}c^{+}_{ph}c^{-}_{ph^{\prime}}\right|.$
(4)

The object is to calculate the matrix element of the Hamiltonian between them.

We write the effective Hamiltonian in Fock space as

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle
E_{0}+\sum_{p}\varepsilon_{p}a^{\dagger}_{p}a_{p}-\sum_{h}\varepsilon_{h}a_{h}a^{\dagger}_{h}+$
(5) $\displaystyle+\sum_{i\neq
j}H_{ij}a^{\dagger}_{i}a_{j}+\sum_{ij,k\ell}v_{ij,k\ell}a^{\dagger}_{i}a^{\dagger}_{j}a_{\ell}a_{k}.$

Here $E_{0}=\langle\Psi_{0}|\hat{H}|\Psi_{0}\rangle$ is the expectation value
of the Hamiltonian in the reference configuration. The effective single-
particle Hamilton $H_{ij}$ has been split into diagonal contributions
$H_{ii}=\varepsilon_{i}$ and off-diagonal ones appearing on the second line.
The last term is the effective two-particle interaction. It contains all the
interaction matrix elements not present in the other terms. We will argue that
its contribution is small when the matrix elements are calculated between
states related by a diabatic transformation.

In a diabatic evolution, each orbit is transformed independently according to
a local single-particle operator such as the quadrupole operator
$\hat{Q}=z^{2}-(x^{2}+y^{2})/2$. Under those conditions the transition density
$\langle\Psi_{q}|\hat{\rho}|\Psi_{q^{\prime}}\rangle=\delta\rho$ is collective
with a strong local component proportional to $d\rho(\mbox{\boldmath$r$})/dq$
in the coordinate space representation. The mean field $V_{0}$ behaves
similarly. The interaction can then be parameterized in the separable form

$v_{ij,k\ell}=\frac{\kappa}{2}\sum_{{\lambda}\mu}O^{{\lambda}\mu*}_{ik}O^{{\lambda}\mu}_{j\ell}$
(6)

with single-particle matrix elements

$O^{{\lambda}\mu}_{ij}=\int
d^{3}r\,\phi^{*}_{i}(\mbox{\boldmath$r$})\phi_{j}(\mbox{\boldmath$r$})\frac{dV_{0}}{dr}Y_{{\lambda}\mu}(\hat{\mbox{\boldmath$r$}}).$
(7)

The strength $\kappa$ is determined by the self-consistency condition bohr-
mottelson ; bertsch1983

$\kappa^{-1}=\int_{0}^{\infty}r^{2}dr\,\frac{d\rho_{0}(r)}{dr}\frac{dV_{0}(r)}{dr}.$
(8)

In calculating the interaction between many-particle configurations, one
includes only the direct term in the interaction; exchange effects are
implicit in the single-particle potential $V_{0}$. Except for $E_{0}$, the
resulting Hamiltonian can be generated entirely in terms of the effective
single-particle mean field and its variation under the diabatic
transformation. Interactions of this form have been applied not only to
collective excitations su81 ; tagami2013  but even to treat nuclear spectra as
a whole alhassid .

We now estimate many-body Hamiltonian matrix elements using Eq. (5-8). Several
assumptions in the derivation are justified by considering the limit of a
large number of particles $N_{h}$ together with small displacements in the
collective coordinate. This leads to a simplification of the Löwdin formula
lowdin1955 ; brink1965  for matrix elements requiring the minors of the
normalization determinant $\det(M)$. The off-diagonal minors are much smaller
than the diagonal ones when the $c_{ph}$ amplitudes are small. Also, the
diagonal ones are approximately equal to $\det(M)$. The Hamiltonian matrix
element can then be evaluated as

$\displaystyle\frac{\langle\Psi_{-q}|\hat{H}|\Psi_{+q}\rangle}{\langle\Psi_{-q}|\Psi_{+q}\rangle}$
$\displaystyle\approx$ $\displaystyle
E_{0}+\sum_{ph}c^{+}_{ph}c^{-}_{ph}(\varepsilon_{p}-\varepsilon_{h})$
$\displaystyle+\sum_{ph}H_{ph}(c^{+}_{ph}+c^{-}_{ph})$
$\displaystyle+\frac{\kappa}{2}\sum_{ph,p^{\prime}h^{\prime}}\left(c^{+}_{ph}c^{-}_{p^{\prime}h^{\prime}}+c^{-}_{ph}c^{+}_{p^{\prime}h^{\prime}}+c^{+}_{ph}c^{+}_{p^{\prime}h^{\prime}}+c^{-}_{ph}c^{-}_{ph}\right)\sum_{{\lambda}\mu}O^{{\lambda}\mu*}_{ph}O^{{\lambda}\mu}_{p^{\prime}h^{\prime}}.$

Eq. (II) can be further simplified in several ways. First, we limit the scope
to axially symmetric diabatic transformations. Then the last sum can be
restricted to $\mu=0$ operators, and the amplitudes can be treated as real
quantities. Next, we can assume that $c^{-}=-c^{+}$ if $q$ is not too large.
Then the last term vanishes. Note that the cancellation requires that the
separable terms in the interaction are even under the time-reversal
transformation. Finally, if the state $|0\rangle$ is a self-consistent
solution of the mean-field equation, the off-diagonal elements $H_{ph}$ vanish
as well.

Changing the notation back to the original, the resulting formula reads

$\frac{\langle\Psi_{q_{1}}|\hat{H}|\Psi_{q_{2}}\rangle}{\langle\Psi_{q_{1}}|\Psi_{q_{2}}\rangle}\approx
E_{0}-\sum_{ph}c^{2}_{ph}(\varepsilon_{p}-\varepsilon_{h})$ (10)

where the particle-hole amplitudes have been determined with respect to a
reference configuration at $q=(q_{1}+q_{2})/2$. It is remarkable that all
dependence on the residual interaction has disappeared.

The formula can be easily applied to calculate the inertial parameter
associated with the collective coordinate $q$. For small $q$ the amplitudes
may be expanded as a power series with the leading term
$c_{ph}=\frac{1}{2}(q_{1}-q_{2})\,\,dc_{ph}/dq$ to give

$\frac{\langle\Psi_{q_{1}}|\hat{H}|\Psi_{q_{2}}\rangle}{\langle\Psi_{q_{1}}|\Psi_{q_{2}}\rangle}\approx
E_{0}-\frac{1}{4}\sum_{ph}\left(\frac{dc_{ph}}{dq}\right)^{2}(q_{1}-q_{2})^{2}(\varepsilon_{p}-\varepsilon_{h}).$
(11)

Finally, one can calculate the collective inertial parameter $B$ by applying
the formula to the plane-wave state

${|k\rangle=\int dq|\Psi_{q}\rangle e^{ikq}.}$ (12)

Writing the derivative as $dc_{ph}/dq=c^{\prime}_{ph}$, one has

${\langle\Psi_{q}|\Psi_{q^{\prime}}\rangle\sim
e^{-\sum_{ph}(c^{\prime}_{ph})^{2}(q-q^{\prime})^{2}/2}}$ (13)

for small $(q-q^{\prime})$. One thus obtains

${\frac{\langle k|\hat{H}|\Psi_{q}\rangle}{\langle k|\Psi_{q}\rangle}\approx
E_{0}}-\frac{\sum_{ph}(c^{\prime}_{ph})^{2}(\varepsilon_{p}-\varepsilon_{h})}{4\sum_{ph}(c^{\prime}_{ph})^{2}}+\frac{{\sum_{ph}(c^{\prime}_{ph}})^{2}(\varepsilon_{p}-\varepsilon_{h})}{4(\sum_{ph}(c^{\prime}_{ph})^{2})^{2}}k^{2}.$
(14)

The third term is the collective kinetic energy $k^{2}/2B$ with $B$ given by

$B=\frac{2(\sum_{ph}(c^{\prime}_{ph})^{2})^{2}}{{\sum_{ph}(c^{\prime}_{ph}})^{2}(\varepsilon_{p}-\varepsilon_{h})}.$
(15)

This agrees with the inertia derived from GCM/GOA using sum rules of a
constraining field (fiolhais83, , Eq.(22)).

##  III Model 1: Translational motion

Here we consider a Hamiltonian of $N$ distinguishable particles interacting
through a translationally invariant quadratic potential in one dimension,

$H=T+v_{2}\\\
=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2m}+\frac{m\omega_{v}^{2}}{2}\sum_{i<j}(x_{i}-x_{j})^{2}.$
(16)

The GCM algebra below generalizes the discussion of the harmonic oscillator
model in Ref. (RS, , Sec. 10.5).

The first task is to find a mean-field ground-state wave function $\Psi_{0}$.
Since the particles are distinguishable, all the orbitals $\phi_{i}$ are the
same and the many-body ground-state configuration has the form

$\Psi_{0}=\prod_{i=1}^{N}\phi(x_{i}).$ (17)

We may assume that the orbitals are Gaussian functions centered at $x=0$,

$\phi(x)=\left(\frac{\alpha}{\pi}\right)^{1/4}e^{-\alpha x^{2}/2},$ (18)

with a common size parameter $\alpha$. It can be verified later that the
Gaussian form allows a self-consistent solution of the mean-field equations.
The expectation value of the Hamiltonian is

$\langle\Psi_{0}|H|\Psi_{0}\rangle=N\frac{\alpha{\hbar^{2}}}{4m}+N(N-1)\frac{m\omega_{v}^{2}}{{4}\alpha}.$
(19)

The first term is from the kinetic energy operator $T$ and the second is from
the interaction $v_{2}$. Minimizing this expression with respect to $\alpha$
yields

$\alpha=(N-1)^{1/2}\omega_{v}m{/\hbar}$ (20)

giving a total energy

$E_{0}=\langle\Psi_{0}|H|\Psi_{0}\rangle=N\frac{(N-1)^{1/2}}{2}{\hbar}\omega_{v}.$
(21)

The next GCM task is to define a generator coordinate; we take it to be a
displacement of the ground-state configuration by an amount $z$,

$\Psi_{z}=\left(\frac{\alpha}{\pi}\right)^{N/4}\prod_{i=1}^{N}e^{-\alpha(x_{i}-z)^{2}/2}{\equiv\prod_{i=1}^{N}\phi_{z}(x_{i}).}$
(22)

Applications of the GCM require the overlap matrix elements
$\langle\Psi_{z}|\Psi_{z^{\prime}}\rangle$ and the Hamiltonian matrix elements
$\langle\Psi_{z}|H|\Psi_{z^{\prime}}\rangle$. The overlaps are given by

$\langle\Psi_{z}|\Psi_{z^{\prime}}\rangle=\exp\left(-N\frac{\alpha}{4}(z-z^{\prime})^{2}\right).$
(23)

Decomposing the matrix elements of the Hamiltonian into the kinetic term and
the interaction term, we have

${\langle\Psi_{z}|T|\Psi_{z^{\prime}}\rangle=}\left(N\frac{\alpha{\hbar^{2}}}{4m}-N\frac{\alpha^{2}{\hbar^{2}}}{8m}(z-z^{\prime})^{2}\right){\langle\Psi_{z}|\Psi_{z^{\prime}}\rangle}$
(24)

for the kinetic term and

$\langle\Psi_{z}|v_{2}|\Psi_{z^{\prime}}\rangle=N(N-1)\frac{m\omega_{v}^{2}}{4\alpha}\langle\Psi_{z}|\Psi_{z^{\prime}}\rangle,$
(25)

for the interaction term. Note that it is the same as the expectation value in
the ground state except for the overlap factor. Combining the two terms, the
Hamiltonian matrix element can be expressed

$\langle\Psi_{z}|H|\Psi_{z^{\prime}}\rangle=\langle\Psi_{z}|\Psi_{z^{\prime}}\rangle\left(E_{0}-h_{2}(z-z^{\prime})^{2}\right)$
(26)

where

$h_{2}=N\frac{\alpha^{2}{\hbar^{2}}}{8m}.$ (27)

This is identical to (RS, , Eq. 10.51) except for the $N$ dependence.

One can derive the effective mass associated with the collective coordinate
from Eq. (23) and (26). It will be shown below that it is just the value for
translational motion, namely $Nm$.

Before that, we verify the orbital energy formula, Eq. (10), choosing end
configurations at $\pm z$ and the reference configuration at $z=0$. The
$\Psi_{z}$ particle-hole amplitudes in the basis of $\Psi_{0}$ are calculated
by taking orbital overlaps. We only need orbitals up to the second excited
state for our purposes. The amplitudes are

${\phi_{z}(x)=\psi_{0}(x-z)=e^{-\alpha
z^{2}/4}\left(\phi_{0}(x)+z\left(\frac{\alpha}{2}\right)^{1/2}\phi_{1}(x)+z^{2}\frac{\alpha}{8^{1/2}}\phi_{2}(x)+...\right).}$
(28)

Next we need the orbital excitation energies, obtained from the spectrum as
$\varepsilon_{n}=n\alpha{\hbar^{2}}/m$. Then the orbital excitation formula
gives

${\sum
c^{+}_{10}c^{-}_{10}(\varepsilon_{1}-\varepsilon_{0})=N\frac{\alpha}{2}\frac{\alpha}{m}\hbar^{2}\cdot\frac{1}{4}(z-z^{\prime})^{2}},$
(29)

in agreement with Eq. (26) with Eq. (27).

Let us now discuss the associated mass term in the collective Hamiltonian. We
can construct a plane-wave state from the GCM configurations according to Eq.
(12). The overlap with $\Psi_{0}$ is

${\langle\Psi_{0}|k\rangle}=\left(\frac{4\pi}{\alpha
N}\right)^{1/2}e^{-k^{2}/\alpha N}$ (30)

and the Hamiltonian matrix element reduces to

${{\langle\Psi_{0}|H|k\rangle}=\langle\Psi_{0}|k\rangle}\left(\frac{k^{2}}{2Nm}-\frac{\alpha}{4m}\right).$
(31)

The first term in parenthesis $k^{2}/2Nm$ is just what we expected for a free
particle of mass $m$. The other term is the zero point kinetic energy
associated with the center-of-mass coordinate in the original GCM
configuration composed of Gaussian orbitals.

It is not generally recognized that the GCM/GOA methodology can deliver exact
composite-particle masses for translational motion. It was noticed as an
empirical finding in Ref. be19  that the translational masses of nuclei
calculated with a Gogny functional were close to $Am$, the number of nucleons
$A$ times the nucleon mass $m$.

##  IV Model 2: Dynamics in a deformation coordinate

We next consider a simplified Nilsson model to illustrate the use of the self-
consistent separable interaction Eq. (6) to generate a two-particle
interaction from a diabatic treatment of the mean field. The single-particle
potential is given by a deformed harmonic oscillator potential,

$V(\mbox{\boldmath$r$})=\frac{1}{2}m\omega_{z}^{2}z^{2}+\frac{1}{2}m\omega_{\perp}^{2}(x^{2}+y^{2}),$
(32)

where $\omega_{z}^{2}=\omega_{0}^{2}(1-4\delta/3)$ and
$\omega_{\perp}^{2}=\omega_{0}^{2}(1+2\delta/3)$, following the notation of
Ref. (bohr-mottelson, , Eq 5-5). This potential can also be expressed as

$V(\mbox{\boldmath$r$})=\frac{1}{2}m\omega_{0}^{2}r^{2}-\frac{{2}}{3}\delta\,m\omega_{0}^{2}\left(z^{2}-(x^{2}+y^{2})/2\right).$
(33)

We take the parameter $\delta$ as the generator coordinate333Another parameter
in common use is $\beta_{2}\approx\frac{2}{3}\sqrt{\frac{4\pi}{5}}\,\delta$..

The eigenstates of the Hamiltonian are specified by a set of quantum numbers
$(n_{x},n_{y},n_{z})$ together with the $z$-components of spin and isospin;
their spatial distributions are determined by $\omega_{0}$ and $\delta$,
combined to give oscillator length parameters $b_{i}=\sqrt{\hbar/m\omega_{i}}$
for $i={z,\perp,0}$.

For simplicity, we consider a harmonic oscillator model of the nucleus 16O.
Its 8 neutrons and 8 protons occupy the “1s” state with
$N_{osc}=n_{x}+n_{y}+n_{z}=0$ and the three “1p” states with $N_{osc}=1$. The
density distribution is given by

$\displaystyle\rho(\mbox{\boldmath$r$})$ $\displaystyle=$
$\displaystyle{4}\sum_{k}|\phi_{k}(\mbox{\boldmath$r$})|^{2}$ (35)
$\displaystyle=$ $\displaystyle\frac{{4}}{\pi b_{\perp}^{2}\sqrt{\pi
b_{z}^{2}}}\left[1+2\left(\frac{x^{2}+y^{2}}{b_{\perp}^{2}}+\frac{z^{2}}{b_{z}^{2}}\right)\right]$
$\displaystyle\times e^{-z^{2}/b_{z}^{2}}e^{-(x^{2}+y^{2})/b_{\perp}^{2}},$

where $k\equiv(n_{x},n_{y},n_{z})=(0,0,0),~{}(0,0,1),~{}(0,1,0)$, and (1,0,0).

To the linear order of $\delta$, the density distribution (35) is expanded to

$\rho(\mbox{\boldmath$r$})=\rho_{0}(r)+\rho_{2}(r)P_{2}({\cos\theta})$ (36)

with

$\rho_{0}(r)=\frac{{4}}{(\pi
b_{0}^{2})^{3/2}}\,\left(1+\frac{2r^{2}}{b_{0}^{2}}\right)e^{-r^{2}/b_{0}^{2}}$
(37)

and

$\rho_{2}(r)={\frac{\delta}{3(\pi
b_{0}^{2})^{3/2}}}\left(\frac{4r^{4}}{b_{0}^{4}}-\frac{2r^{2}}{b_{0}^{2}}\right)\,e^{-r^{2}/b_{0}^{2}}.$
(38)

Given the form of the potential (33), the many-body Hamiltonian with a
separable interaction reads

$H=\sum_{i}\left(\frac{\mbox{\boldmath$p$}_{i}^{2}}{2m}+\frac{1}{2}m\omega_{0}^{2}r_{i}^{2}\right)+\frac{\kappa}{2}\left(\sum_{i}r_{i}^{2}P_{2}({\cos\theta_{i}})\right)^{2}.$
(39)

The strength $\kappa$ evaluated with Eq. (8) is

$\kappa=-{\frac{2}{3}}\,\delta
m\omega_{0}^{2}\left({\frac{4\pi}{5}}\int^{\infty}_{0}r^{2}dr\,r^{2}\rho_{2}(r)\right)^{-1}=-\frac{m^{2}\omega_{0}^{3}}{{18}}.$
(40)

$n$ | $\langle nb|nb^{\prime}\rangle$ | $\langle nb|x^{2}|nb^{\prime}\rangle$ | $\langle nb|p^{2}|nb^{\prime}\rangle$
---|---|---|---
0 | $\sqrt{\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}}$ | $\frac{bb^{\prime}}{2}\,\left(\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}\right)^{3/2}$ | $\frac{\hbar^{2}}{b^{2}+b^{\prime 2}}\,\left(\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}\right)^{1/2}$
1 | $\left(\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}\right)^{3/2}$ | $\frac{3bb^{\prime}}{2}\,\left(\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}\right)^{5/2}$ | $\frac{3\hbar^{2}}{b^{2}+b^{\prime 2}}\,\left(\frac{2bb^{\prime}}{b^{2}+b^{\prime 2}}\right)^{3/2}$
Table 1:  Overlaps $\langle
nb|nb^{\prime}\rangle{\equiv}\langle\phi_{n}(b)|\phi_{n}(b^{\prime})\rangle$
and needed matrix elements for one-dimensional harmonic oscillator wave
functions with oscillator lengths $b$ and $b^{\prime}$. The oscillator quanta
are denoted by $n$.

Using orbital overlaps in Table I, it is easy to evaluate the many-body
overlap and Hamiltonian matrix element,

$\displaystyle\langle\Psi_{{\delta}}|\Psi_{{\delta^{\prime}}}\rangle$
$\displaystyle=$
$\displaystyle\left(\frac{2b_{\perp}b_{\perp}^{\prime}}{b_{\perp}^{2}+b_{\perp}^{\prime
2}}\right)^{{24}}\left(\frac{2b_{z}b_{z}^{\prime}}{b_{z}^{2}+b_{z}^{\prime
2}}\right)^{{12}},$ (41)
$\displaystyle\frac{\langle\Psi_{{\delta}}|H|\Psi_{{\delta^{\prime}}}\rangle}{\langle\Psi_{{\delta}}|\Psi_{{\delta^{\prime}}}\rangle}$
$\displaystyle=$
$\displaystyle\frac{{24}\hbar^{2}}{2m}\,\left(\frac{2}{b_{\perp}^{2}+b_{\perp}^{\prime
2}}+\frac{1}{b_{z}^{2}+b_{z}^{\prime 2}}\right)$
$\displaystyle+\frac{1}{2}m\omega_{0}^{2}\left(\frac{{48}b_{\perp}^{2}b_{\perp}^{\prime
2}}{b_{\perp}^{2}+b_{\perp}^{\prime 2}}+\frac{{24}b_{z}^{2}b_{z}^{\prime
2}}{b_{z}^{2}+b_{z}^{\prime 2}}\right)$
$\displaystyle+\frac{\kappa}{2}\left[24\left(\frac{b_{\perp}^{2}b_{\perp}^{\prime
2}}{b_{\perp}^{2}+b_{\perp}^{\prime 2}}-\frac{b_{z}^{2}b_{z}^{\prime
2}}{b_{z}^{2}+b_{z}^{\prime 2}}\right)\right]^{2}.$

The result for the overlap can also be obtained as a special case of the
general formula for many-particle harmonic oscillator configurations given in
Ref. (arima1959, , Eq.3.5). 444The variable $\varepsilon$ in the notation of
Ref. arima1959  is identical $(\delta-\delta^{\prime})/6$.

Finally we expand Eq. (LABEL:eq:Hov) to the quadratic order of $\delta$ and
$\delta^{\prime}$ to obtain

$\frac{\langle\Psi_{{\delta}}|H|\Psi_{{\delta^{\prime}}}\rangle}{\langle\Psi_{{\delta}}|\Psi_{{\delta^{\prime}}}\rangle}=18\hbar\omega_{0}-\hbar\omega_{0}(\delta-\delta^{\prime})^{2}+...$
(43)

One thus finds

$h_{2}=\hbar\omega_{0}.$ (44)

$\psi_{000}$ | $-{8}\hbar\omega_{0}[2(a_{2}^{(\perp)})^{2}+(a_{2}^{(z)})^{2}]$
---|---
$\psi_{001}$ | $-{8}\hbar\omega_{0}[2(a_{2}^{(\perp)})^{2}+(a_{3}^{(z)})^{2}]$
$\psi_{100}$ | $-{8}\hbar\omega_{0}[(a_{3}^{(\perp)})^{2}+(a_{2}^{(\perp)})^{2}+(a_{2}^{(z)})^{2}]$
$\psi_{010}$ | $-{8}\hbar\omega_{0}[(a_{3}^{(\perp)})^{2}+(a_{2}^{(\perp)})^{2}+(a_{2}^{(z)})^{2}]$
Table 2:  The contribution of each single-particle wave function $\psi_{k}$ to
the second order term in the energy overlap,
$-\sum_{ph}|C_{ph}|^{2}(\epsilon_{p}-\epsilon_{h})$. Here, the particle-hole
energy is $\epsilon_{p}-\epsilon_{h}=2\hbar\omega_{0}$, and the spin
degeneracy is also taken into account.

We now evaluate $h_{2}$ in the orbital energy formula Eq. (11) to compare with
Eq. (44). We apply the formula to the matrix elements between
$\Psi_{{\delta}}$ and $\Psi_{{-\delta}}$ taking $\Psi_{{0}}$ as the reference
configuration. The main task is to determine the coefficients $c_{ph}$ of the
deformed orbits in terms of the orbitals in the spherical basis. Since we are
interested in the second order behavior of the Hamiltonian overlap with
respect to $\delta$, it is sufficient to expand the wave functions up to $n=2$
and $n=3$ for $\phi_{0}(b,x)$ and $\phi_{1}(b,x)$, respectively. This leads to

$\displaystyle\phi_{0}(b,x)$ $\displaystyle\propto$
$\displaystyle\phi_{0}(b_{0},x)+a_{2}\phi_{2}(b_{0},x)+\cdots$ (45)
$\displaystyle\phi_{1}(b,x)$ $\displaystyle\propto$
$\displaystyle\phi_{1}(b_{0},x)+a_{3}\phi_{3}(b_{0},x)+\cdots$ (46)

with the coefficients of

$a_{2}=\frac{1}{\sqrt{2}}\,\frac{b^{2}-b_{0}^{2}}{b^{2}+b_{0}^{2}},$ (47)

and

$a_{3}=\sqrt{\frac{3}{2}}\,\frac{b^{2}-b_{0}^{2}}{b^{2}+b_{0}^{2}}$ (48)

for the one-dimensional Gaussians in an orbital. The three-dimensional
orbitals $\phi_{n_{x},n_{y},n_{z}}$ are expanded as

$\displaystyle\phi_{000}$ $\displaystyle\propto$
$\displaystyle\phi_{000}^{(0)}+a_{2}^{(\perp)}(\phi_{200}^{(0)}+\phi_{020}^{(0)})+a_{2}^{(z)}\phi_{002}^{(0)}+\cdots$
$\displaystyle\phi_{001}$ $\displaystyle\propto$
$\displaystyle\phi_{001}^{(0)}+a_{2}^{(\perp)}(\phi_{201}^{(0)}+\phi_{021}^{(0)})+a_{3}^{(z)}\phi_{003}^{(0)}+\cdots$
$\displaystyle\phi_{100}$ $\displaystyle\propto$
$\displaystyle\phi_{100}^{(0)}+a_{3}^{(\perp)}\phi_{300}^{(0)}+a_{2}^{(\perp)}\phi_{120}^{(0)}+a_{2}^{(z)}\phi_{102}^{(0)}+\cdots$

where $\phi_{k}^{(0)}$ is a spherical orbital. The contribution of each
single-particle state to the second order term in the energy overlap,
$-\sum_{ph}|C_{ph}|^{2}(\epsilon_{p}-\epsilon_{h})$, is shown in Table 2. From
this table, one has

$\displaystyle-\sum_{ph}|C_{ph}|^{2}(\epsilon_{p}-\epsilon_{h})$ (53)
$\displaystyle=$
$\displaystyle-{8}\hbar\omega_{0}[6(a_{2}^{(\perp)})^{2}+2(a_{3}^{(\perp)})^{2}+3(a_{2}^{(z)})^{2}+(a_{3}^{(z)})^{2}]$
$\displaystyle=$
$\displaystyle-{8}\hbar\omega_{0}\left\\{6\left(\frac{b_{\perp}^{2}-b_{0}^{2}}{b_{\perp}^{2}+b_{0}^{2}}\right)^{2}+3\left(\frac{b_{z}^{2}-b_{0}^{2}}{b_{z}^{2}+b_{0}^{2}}\right)^{2}\right\\}.$

Expanding this quantity up to the second order of $\delta$, one finds

$-\sum_{ph}|C_{ph}|^{2}(\epsilon_{p}-\epsilon_{h})\sim-{4}\hbar\omega_{0}\delta^{2}.$
(54)

This coincides with the second order term in Eq. (43) with
$\delta^{\prime}=-\delta$. It is easy to confirm that the same relation holds
also with a deformed reference configuration, with $b_{0}=b(\delta_{\rm ref})$
and $b=b(\delta_{\rm ref}+\delta)$.

##  V Summary

We have derived a simple formula for the Hamiltonian matrix elements by the
Generator Coordinate Method. The formula is based on a residual interaction of
separable form, determined by a diabatic treatment of the generator
coordinate. The formula was shown to be exact for the leading dependence on
$q_{1}-q_{2}$ for two solvable models.

Although it is approximate, the formula has several attractive features.
First, it does not require full details of the many-body Hamiltonian. Thus one
can carry out GCM calculations even when only the mean-field potential is
known. This is the case when using a phenomenological mean-field potential
such as a Woods-Saxon potential. Second, the formula is much simpler than the
original multi-step procedure based on the generalized Wick theorem or the
Löwdin formula. Third, by using a separable interaction, one can avoid the
well-known difficulties of treating energy functionals as Hamiltonians.

We plan to report numerical calculations of more realistic models in a
separate publication.

## Acknowledgements

We thank P.-G. Reinhard for comments on the manuscript. This work was
supported in part by JSPS KAKENHI Grant Numbers JP19K03861 and JP21H00120.

## References

  * (1) P. Ring and P. Schuck, The Nuclear Many Body Problem (Springer-Verlag, New York, 1980).
  * (2) P.-G. Reinhard and K. Goeke, Rep. Prog. in Phys. 50, 1 (1987).
  * (3) L.M. Robledo, T.R. Rodriguez, and R.R. Rodriguez-Guzman, J. of Phys. G46, 013001 (2019).
  * (4) R. Balian and E. Brézin, Nuovo Cimento 64B, 34 (1969).
  * (5) P. Bonche et al., Nucl. Phys. A510, 466 (1990).
  * (6) M Baldo, L.M. Robledo, P. Schuck and X. Vinas, Phys. Rev. C 87, 064305 (2013).
  * (7) R. Roth, S. Binder, K. Vobig, et al., Phys. Rev. Lett. 109, 052501 (2012).
  * (8) D.J. Thouless, Nucl. Phys. 21, 225 (1960).
  * (9) A. Bohr and B.R. Mottelson, Nuclear Structure, Vol. II (Benjamin, New York, 1975).
  * (10) G.F. Bertsch, Prog. Theo. Phys. Suppl. 74&75, 115 (1983).
  * (11) T. Suzuki and H. Sagawa, Prog. Theo. Phys. 65, 565 (1981).
  * (12) S. Tagami, Y.R. Shimizu, and J. Dudek, Phys. Rev. C87, 054306 (2013).
  * (13) Y. Alhassid, L. Fang, and H. Nakada, Phys. Rev. Lett. 101, 082501 (2008).
  * (14) P.O. Löwdin, Phys. Rev. 97, 1490 (1955).
  * (15) D.M. Brink, Proc. International School of Physics Enrico Fermi’ course XXXVI (1965), p. 247.
  * (16) C. Fiolhais and R.M. Dreizler, Nucl. Phys. A393, 205 (1983).
  * (17) G.F. Bertsch and W. Younes, Ann. Phys. 403, 68 (2019).
  * (18) A. Arima and S. Yoshida, Nucl. Phys. 12, 139 (1959).


©2024 IEEE. Published in 2024 International Conference on Acoustics, Speech
and Signal Processing (ICASSP), scheduled for 14-19 April 2024 in Seoul,
Korea. Personal use of this material is permitted. However, permission to
reprint/republish this material for advertising or promotional purposes or for
creating new collective works for resale or redistribution to servers or
lists, or to reuse any copyrighted component of this work in other works, must
be obtained from the IEEE. DOI: tba

# A Guided Upsampling Network for Short Wave Infrared Images Using Graph
Regularization

###### Abstract

Exploiting the infrared area of the spectrum for classification problems is
getting increasingly popular, because many materials have characteristic
absorption bands in this area. However, sensors in the short wave infrared
(SWIR) area and even higher wavelengths have a very low spatial resolution in
comparison to classical cameras that operate in the visible wavelength area.
Thus, in this paper an upsampling method for SWIR images guided by a visible
image is presented. For that, the proposed guided upsampling network (GUNet)
uses a graph-regularized optimization problem based on learned affinities is
presented. The evaluation is based on a novel synthetic near-field visible-
SWIR stereo database. Different guided upsampling methods are evaluated, which
shows an improvement of nearly 1 dB on this database for the proposed
upsampling method in comparison to the second best guided upsampling network.
Furthermore, a visual example of an upsampled SWIR image of a real-world scene
is depicted for showing real-world applicability.

Index Terms—  Image Processing, Deep Learning, Short Wave Infrared Imaging,
Guided Upsampling

##  1 Introduction

The infrared area of the spectrum is particularly interesting for a lot of
real world classification problems, since different materials often have
unique spectral fingerprints in this wavelength range. For example, it can be
used in agriculture [1] to retrieve the plant health, in medicine to determine
the degree of burn [2], in the area of recycling to discriminate between
different types of plastic [3], or in difficult imaging scenarios like haze,
since longer wavelengths are advantageous in case of Rayleigh scattering [4].
Due to the different spectral area from roughly from 1000 nm to 2000 nm that
is recorded by SWIR cameras, they typically have indium gallium arsenide
(InGaAs) sensors. Caused by the worse thermal conductivity of InGaAs sensors,
the dark current noise is much stronger in comparison to classical visual
range of spectrum cameras [5]. Therefore, these cameras typically have huge
pixel sizes and hence a much lower spatial resolution. However, a high spatial
resolution is highly desirable in many applications to reveal details and help
classification and tracking algorithms. Hence, the goal of the paper is to
upsample the SWIR image by exploiting the structure of a corresponding high
resolution visible image, which typically shows the spectral area from 400 nm
to 700 nm.

(a) Setup(b) SWIR image(c) Visible image Fig. 1: The built stereo setup (a)
producing an SWIR image (b) a visible image (c).
VISDeepFeatureExtractorImageFeatures ${\bm{f}}$AffinitygraphGraph-
basedOptimiz.SWIR LRSWIR HR Fig. 2: The pipeline of the proposed guided
upsampling network (GUNet).

To test the real-world guided upsampling capability, a visible-SWIR stereo
camera setup was built as shown in Fig. 1. The SWIR camera has a resolution of
320 × 256 pixels with a pixel size of 30 µm × 30 µm, while the visible camera
has a resolution of 2448 × 2048 using pixels of size 3.45 µm × 3.45 µm. The
lenses of both cameras have a focal length of 16 mm. Due to the slightly
bigger sensor size, the SWIR camera is able to capture a wider field of view.

The visible-SWIR stereo setup shown in Fig. 1 operates in near-field, which
means that the objects are relatively close in comparison to, e.g., an
airborne device. In near-field imaging, the objects in the scene have a depth-
dependent offset, also called disparity, from the perspective of spatially
distributed cameras. Moreover, due to the different viewing angles on objects,
the different cameras also have a different perspective on objects and see
different parts of the background behind an object. Hence, a registration and
reconstruction process is necessary. For this, an image processing pipeline
like introduced by Genser et al. [6] can be used. The assumption for the
proposed method is that such a pipeline was already deployed.

The task of guided upsampling is popular for increasing the resolution of
sparse depth maps based on a high resolution image. There, the goal is to
upsample a sparse depth map using a high resolution RGB image as guide [7, 8,
9]. In this paper, the idea of using affinities [9] is examined, since SWIR
images typically have less texture than its corresponding visible image. This
is visible in Fig. 1, where the logo of the cup is not visible in the SWIR
image. Moreover, an affinity-based upsampling is more robust towards
registration and reconstruction errors made by the image processing pipeline
for overlaying the SWIR and visible image.

##  2 Guided Upsampling

$\begin{pmatrix}.\overline{4}&.\overline{2}&0&.\overline{2}&.\overline{1}&0&0&0&0\\\
0&.\overline{2}&.\overline{4}&0&.\overline{1}&.\overline{2}&0&0&0\\\
0&0&0&.\overline{2}&.\overline{1}&0&.\overline{4}&.\overline{2}&0\\\
0&0&0&0&.\overline{1}&.\overline{2}&0&.\overline{2}&.\overline{4}\end{pmatrix}$${\bm{D}}$
Fig. 3: Illustration of the downsampling matrix ${\bm{D}}$. A high resolution
$3\times 3$ grid is downscaled to a $2\times 2$ image. The low resolution
pixels are shown in different colors.

The proposed guided upsampling network (GUNet) is based on a work for guided
depth upsampling by Lutio et al. [9]. In this paper, this network is improved
and modified to upsample SWIR image data using the visible image as guide. For
that, a more general affinity function is introduced and the optimal working
point is examined. Moreover, arbitrary scaling factors are made possible by
introducing a more general downsampling operator. Finally, it is proven that
the optimization problem can be backpropagated, and thus the network is end-
to-end trainable. It is assumed that the images are already registered and
reconstructed using e.g. [6].

###  2.1 Network Architecture

As depicted in Fig. 2, the first step is to extract pixel-wise deep features
${\bm{f}}$ of length ${M}$ from the visible image using a UNet [10]
architecture, namely ResNet50 [11]. In this paper, the length ${M}$ of the
features is set to 64. With the help of these features, the affinity matrix
${\bm{A}}$, which contains the affinity for each pixel to its four direct
neighbors, can be calculated by the similarity function

${}A_{ij}=\mathrm{e}^{-\frac{{d}({f}_{i},{f}_{j})}{{\mu}}},$ (1)

where ${\mu}$ is a learnable scaling parameter and ${d}({f}_{i},{f}_{j})$ is a
distance function between features. This distance function is discussed in
Section 2.2. ${\bm{A}}$ can be interpreted as adjacency matrix of an affinity
graph. Afterwards, the Laplacian matrix ${\bm{L}}={\bm{U}}-{\bm{A}}$ is
calculated, where the degree matrix ${\bm{U}}$ is a diagonal matrix with the
entries being ${U}_{ii}=\sum_{j}{}A_{ij}$. Finally, the optimization problem

$\operatorname*{\arg\\!\min}_{{\bm{y}}}\mathcal{L}_{\text{rec}}({\bm{y}},{\bm{L}})=\operatorname*{\arg\\!\min}_{{\bm{y}}}||{\bm{D}}{\bm{y}}-{\bm{s}}||_{2}^{2}+\lambda{\bm{y}}^{\text{T}}{\bm{L}}{\bm{y}}$
(2)

is solved, where the first term is the data term and contains the downsampling
matrix ${\bm{D}}$, the vectorized low resolution SWIR image ${\bm{s}}$ and the
vectorized high resolution SWIR image ${\bm{y}}$ to determine. The second term
is the regularizer, which implicitly incorporates the affinity matrix and thus
the structure of the high resolution visible image. The form of this term
results from any smoothness regularizer. In this case, the estimated signal
should be smooth on the affinity graph. $\lambda$ is a learnable trade-off
parameter steering how close to stay to the original image. Since the scale
factor between warped SWIR image and the visible image is non-integer with a
very high probability, ${\bm{D}}$ is able to split a high resolution pixel to
influence several low resolution pixels. For example, for a scale factor of
$1.5$, the first low resolution pixel will contain the content of the first
high resolution pixel, half of the pixels to the right and to the bottom, and
a quarter of the diagonal pixel. This example is depicted in Fig. 3. The
solution to the optimization problem is calculated by solving the linear
systems of equations

$\left({\bm{D}}^{\text{T}}{\bm{D}}+\lambda{\bm{L}}\right){\bm{y}}={\bm{D}}^{\text{T}}{\bm{s}}.$
(3)

The optimization problem shown in (2) can be backpropagated. The gradients
through the optimization problem can be calculated using the implicit function
theorem [12]. In general, the optimization problem for training contains the
optimization problem shown in (2)

${\bm{\theta}}^{*}=\operatorname*{\arg\\!\min}_{{\bm{\theta}}}\mathcal{L}_{\text{train}}\left(\operatorname*{\arg\\!\min}_{{\bm{y}}}\mathcal{L}_{\text{rec}}({\bm{y}},{\bm{L}}({\bm{\theta}}))\right),$
(4)

where it is being made explicit, that only the Laplacian matrix depends on the
network parameters ${\bm{L}}$. The gradients for this type of loss can be
calculated using the implicit function theorem. Suppose ${\bm{L}}$ and
${\bm{y}}$ are related to each other through function $g$

$g({\bm{y}},{\bm{L}})=0.$ (5)

Then, by assuming $g({\bm{y}},{\bm{L}})$ to be smooth, if ${\bm{y}}$ is
slightly changed by $\Delta_{\bm{y}}$, ${\bm{L}}$ is also slightly changed by
$\Delta_{\bm{L}}$ to preserve the constraint

$g(\Delta_{\bm{y}},\Delta_{\bm{L}})=0.$ (6)

Taking the first-order expansion of this leads to

$g({\bm{y}},{\bm{L}})+\Delta_{\bm{y}}\frac{\delta
g}{\delta{\bm{y}}}+\Delta_{\bm{L}}\frac{\delta g}{\delta{\bm{L}}}=0.$ (7)

Since (5) holds, the relationship

$\Delta_{\bm{y}}\frac{\delta g}{\delta{\bm{y}}}=-\Delta_{\bm{L}}\frac{\delta
g}{\delta{\bm{L}}}$ (8)

is established. Finally,

$\frac{\Delta_{\bm{y}}}{\Delta_{\bm{L}}}=-\left(\frac{\delta
g}{\delta{\bm{y}}}\right)^{-1}\frac{\delta g}{\delta{\bm{L}}}$ (9)

holds. Here, the function
$g({\bm{y}},{\bm{L}})=\frac{\delta\mathcal{L}_{\text{rec}}({\bm{y}},{\bm{L}})}{\delta{\bm{y}}}$
is fulfilled, which leads to

$\displaystyle\frac{\delta{\bm{y}}}{\delta{\bm{L}}}$
$\displaystyle=-\left(\frac{\delta^{2}\mathcal{L}_{\text{rec}}({\bm{y}},{\bm{L}})}{\delta{\bm{y}}\delta{\bm{y}}^{\text{T}}}\right)^{-1}\frac{\delta^{2}\mathcal{L}_{\text{rec}}({\bm{y}},{\bm{L}})}{\delta{\bm{y}}\delta{\bm{L}}^{\text{T}}}$
(10)
$\displaystyle=-\left({\bm{D}}^{\text{T}}{\bm{D}}+\lambda{\bm{L}}\right)^{-1}\lambda{\bm{y}}$

Instead of calculating the inverse directly, a linear system of equations
$\left({\bm{D}}^{\text{T}}{\bm{D}}+\lambda{\bm{L}}\right)\bm{x}=\lambda{\bm{y}}$
is solved very similar to (3). Now, this gradient can be used in the chain
rule to get the gradient of the training loss function with respect to the
network parameters

$\frac{\delta\mathcal{L}_{\text{train}}}{\delta{\bm{L}}}=\frac{\delta\mathcal{L}_{\text{train}}}{\delta{\bm{y}}}\frac{\delta{\bm{y}}}{\delta{\bm{L}}}.$
(11)

This shows that the whole network is end-to-end trainable.

The problem is that no large scale visible-SWIR training data is available.
Therefore, the RGB database Places [13] is used and the cross spectral data is
augmented. For this, the same data augmentation is used as in [14]. This data
augmentation exploits the HSV color space to assign random grayscale values to
a couple of colors. In between these colors, the remaining colors are linearly
interpolated using the random grayscale values. Patches of size $256\times
256$ were used as well as a downscale factor of 8. Adam optimizer with
parameters $\beta_{1}=0.5$ and $\beta_{2}=0.999$ and a learning rate of
$0.0001$ was used. The best model on training data after one epoch is
deployed.

###  2.2 Affinity Function

Table 1: Average PSNR and SSIM of different distance functions $d_{o}({f}_{i},{f}_{j})$ on the hyperspectral database CAVE. | $o=1$ | $o=1.5$ | $o=2$ | $o=4$ | $o=10$
---|---|---|---|---|---
PSNR | 32.13 | 32.32 | 32.02 | 32.25 | 31.88
SSIM | 0.890 | 0.892 | 0.889 | 0.889 | 0.880

One of the very few choices to make in this architecture is the affinity and
thus distance function. For this, the network was trained using several
distance functions. The distance functions have the general form of

$d_{o}({f}_{i},{f}_{j})=\frac{1}{{M}}\sum_{m=1}^{{M}}|{f}_{i,m}-{f}_{j,m}|^{{o}}$
(12)

and were used at points ${o}\in\\{1,1.5,2,4,10\\}$.

To optimize this parameter, the individual trained networks were evaluated on
the hyperspectral database CAVE [15]. For this an image at wavelength 500 nm
served as guide while the image at wavelength 650 nm was scaled down by factor
8 and used as low resolution image. Note that the reconstruction and
registration process necessary for visible-SWIR stereo imaging is not part of
the evaluation here. Only the upsampling performance is considered for finding
the best distance function.

The results are summarized in Table 1. The different distance functions
$d_{o}({f}_{i},{f}_{j})$ all work well on the CAVE database. However,
parameters ${o}=1.5$ and ${o}=4$ work better than the others. Due to the
slightly better performance, ${o}=1.5$ is chosen in the upcoming evaluation.

##  3 Evaluation

The evaluation is split into two parts. First, a quantitative evaluation based
on a novel visible-SWIR database is shown. Afterwards, a qualitative
evaluation is performed using a record from a real-world visible-SWIR stereo
setup.

###  3.1 Quantitative Evaluation

Since it is impossible to record ground-truth data using the presented setup,
a synthetic visible-SWIR database is created based on the database HyViD [16],
which contains wavelengths from 400 nm to 700 nm. This database was created
using the 3D modeling software Blender. Textures were extracted from a real-
world database containing scenes from different environments. Light sources
were emulated using Planck’s law. The database contains seven moving scenes,
each with 30 frames. For this evaluation, a visible camera and a SWIR camera
in a stereo setup was synthetically created and inserted into the scenes. The
visible camera was rendered at 500 nm, while the SWIR camera was rendered at
650 nm, since no SWIR textures are available. Apart from the focal length (6
mm), all other aspects stayed as close to the real setup as possible. An
example frame from both cameras of the scene city is shown in Fig. 4.

Note that the cameras are perfectly aligned in Blender, and thus the
quantitative evaluation skips the calibration process. Instead, since the SWIR
camera has a slightly bigger field of view, the parts of the sensor, which are
not visible to the visible camera, are cropped away. These areas were
identified by comparing physical sensor sizes.

Fig. 4: The novel synthetic visible-SWIR stereo database. The image on the left was rendered from the low-resolution SWIR camera, while the right image depicts the high resolution visible image. Table 2: Average PSNR in dB and SSIM of different upsampling methods on the novel synthetic database. | BIC | VDSR | HAN | SwinIR | GAD | FDKN | GUNet
---|---|---|---|---|---|---|---
|  | [17] | [18] | [19] | [7] | [8] |
PSNR | 25.15 | 25.14 | 24.51 | 25.03 | 23.83 | 26.40 | 27.31
SSIM | 0.755 | 0.754 | 0.629 | 0.754 | 0.745 | 0.820 | 0.834

GUNet is compared against single image super resolution methods bicubic
interpolation (BIC), Very Deep Super Resolution [17] (VDSR), Holistic
Attention Network [18] (HAN) and Swin Image Restoration [19] (SwinIR), as well
as guided upsampling methods Guided Anisotropic Diffusion [7] (GAD), and Fast
Deformable Kernel Networks [8] (FDKN). The guided upsampling networks GAD and
FDKN were retrained using the same procedure as GUNet, since they all
originate from guided depth upsampling.

The results in terms of average PSNR and SSIM are summarized in Table 2. The
single image super resolution methods are able to upsample the SWIR image, but
cannot keep up with the performance of the guided methods according to PSNR
and SSIM. Surprisingly, the bicubic upsampled version performs best in the
group of single-image super resolution methods. This may originate from the
fact, that these methods were trained on only downsampled versions of the
images to upsample. Hence, registration and reconstruction errors are
interpreted as details that needs to be sharpened and thus overexaggerated,
while BIC rather smoothens these errors. A similar effect is happening with
the guided diffusion-based network GAD, which is performing much worse than
just a simple bicubic upscaling. On the other hand, FDKN and GUNet are able to
properly outperform the bicubic upscaled version. The proposed GUNet performs
best according to PSNR and SSIM.

BICVisibleSwinIRGADFDKNGUNet Fig. 5: The qualitative evaluation on a real-
world scene.

###  3.2 Qualitative Evaluation

For the qualitative evaluation, a real-world scene was recorded using the
presented setup in Fig. 1. In Fig. 5, the results of the bicubic upsampled
image, the visible guide image, various upsampling methods are depicted. Due
to the affinity-based optimization problem, GUNet provides sharper edges than
all other methods. This is well visible for the tree trunk, the color
calibration chart and the plant. Note that GUNet can even better conceal some
of the reconstruction errors made by the cross spectral reconstruction, e.g.,
to the left of the car.

##  4 Conclusion

In this paper, a novel near-field visible-SWIR stereo camera setup was
introduced. Due to unfavorable thermal properties, the SWIR sensor has a much
lower resolution than the visible sensor. Therefore, a guided upsampling
network using a graph-regularized optimization problem was presented
exploiting the high resolution visible image. In the evaluation, the proposed
affinity-based network outperformed its learned competitors trained with the
same procedure by 1 dB and gains 2 dB over bicubic upscaling. Moreover, this
guided upsampling method also provides satisfying results on a real-world
visible-SWIR record.

## References

  * [1] Matheus Cardim Ferreira Lima, Anne Krus, Constantino Valero, Antonio Barrientos, Jaime del Cerro, and Juan Jesús Roldán-Gómez,  “Monitoring plant status and fertilization strategy through multispectral images,”  Sensors, vol. 20, no. 2, 2020.
  * [2] Michael G. Sowa, Lorenzo Leonardi, Jeri R. Payette, K. M. Cross, Manuel Gomez, and Joel Fish,  “Classification of burn injuries using near-infrared spectroscopy,”  Journal of Biomedical Optics, vol. 11, no. 5, pp. 054002, 2006.
  * [3] Monica Moroni, Alessandro Mei, Alessandra Leonardi, Emanuela Lupo, and Floriana La Marca,  “PET and PVC Separation with Hyperspectral Imagery,”  Sensors, vol. 15, no. 1, pp. 2205–2227, Jan. 2015.
  * [4] Marc P. Hansen and Douglas S. Malchow,  “Overview of SWIR detectors, cameras, and applications,”  in Thermosense XXX. International Society for Optics and Photonics, 2008, vol. 6939, p. 69390I, SPIE.
  * [5] Michael MacDougal, Jon Geske, Chad Wang, Shirong Liao, Jonathan Getty, and Alan Holmes,  “Low dark current InGaAs detector arrays for night vision and astronomy,”  in Infrared Technology and Applications XXXV, Bjørn F. Andresen, Gabor F. Fulop, and Paul R. Norton, Eds. International Society for Optics and Photonics, 2009, vol. 7298, p. 72983F, SPIE.
  * [6] Nils Genser, Jürgen Seiler, and André Kaup,  “Camera array for multi-spectral imaging,”  IEEE Transactions on Image Processing, vol. 29, pp. 9234–9249, 2020\.
  * [7] Nando Metzger, Rodrigo Caye Daudt, and Konrad Schindler,  “Guided depth super-resolution by deep anisotropic diffusion,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2023, pp. 18237–18246.
  * [8] Beomjun Kim, Jean Ponce, and Bumsub Ham,  “Deformable kernel networks for joint image filtering,”  in International Journal of Computer Vision volume, 2021, vol. 129, pp. 579–600.
  * [9] Riccardo de Lutio, Alexander Becker, Stefano D’Aronco, Stefania Russo, Jan D. Wegner, and Konrad Schindler,  “Learning graph regularisation for guided super-resolution,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2022, pp. 1979–1988.
  * [10] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, Eds., Cham, 2015, pp. 234–241, Springer International Publishing.
  * [11] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016, pp. 770–778.
  * [12] Krisorn Jittorntrum,  “An implicit function theorem,”  Journal of Optimization Theory and Applications volume, vol. 25, no. 4, pp. 575–577, 1978.
  * [13] Bolei Zhou, Agata Lapedriza, Aditya Khosla, Aude Oliva, and Antonio Torralba,  “Places: A 10 million image database for scene recognition,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, no. 6, pp. 1452–1464, 2018.
  * [14] Frank Sippel, Jürgen Seiler, and André Kaup,  “Cross spectral image reconstruction using a deep guided neural network,”  in 2023 IEEE International Conference on Image Processing (ICIP), 2023, pp. 226–230.
  * [15] Fumihito Yasuma, Tomoo Mitsunaga, Daisuke Iso, and Shree K. Nayar,  “Generalized assorted pixel camera: Postcapture control of resolution, dynamic range, and spectrum,”  IEEE Transactions on Image Processing, vol. 19, no. 9, pp. 2241–2253, 2010.
  * [16] Frank Sippel, Jürgen Seiler, and André Kaup,  “Synthetic hyperspectral array video database with applications to cross-spectral reconstruction and hyperspectral video coding,”  J. Opt. Soc. Am. A, vol. 40, no. 3, pp. 479–491, Mar. 2023.
  * [17] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee,  “Accurate image super-resolution using very deep convolutional networks,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016.
  * [18] Ben Niu, Weilei Wen, Wenqi Ren, Xiangde Zhang, Lianping Yang, Shuzhen Wang, Kaihao Zhang, Xiaochun Cao, and Haifeng Shen,  “Single image super-resolution via a holistic attention network,”  in Computer Vision – ECCV 2020, Cham, 2020, pp. 191–207, Springer International Publishing.
  * [19] Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte,  “SwinIR: Image restoration using swin transformer,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV) Workshops, Oct. 2021, pp. 1833–1844.