# Predictions for feed-down enhancements at the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ thresholds via the triangle and box singularities

Ming-Xiao Duan1<EMAIL_ADDRESS>Lin Qiu1,2<EMAIL_ADDRESS>Xi-Zhe Ling1
<EMAIL_ADDRESS>Qiang Zhao1,2<EMAIL_ADDRESS>1 Institute of High Energy
Physics, Chinese Academy of Sciences, Beijing 100049, P.R. China and
2 University of Chinese Academy of Sciences, Beijing 100049, P.R. China

###### Abstract

We demonstrate that triangle singularity (TS) and box singularity (BS)
mechanisms can produce unique narrow enhancements at the $\Lambda_{c}\bar{D}$
and $\Lambda_{c}\bar{D}^{*}$ thresholds in the invariant mass spectra of
$J/\psi p$ and $J/\psi p\pi$, respectively. Taking into account that such
mechanisms only depend on the initial $\Sigma_{c}^{(*)}\bar{D}^{(*)}$
interactions near threshold within the TS or BS kinematic regimes, the
$\Lambda_{c}\bar{D}$ and $\Lambda_{c}\bar{D}^{*}$ threshold enhancements can
be regarded as a feed-down phenomenon originated from both the heavier
pentaquark decays and the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ scatterings from the
continuum. A search for these structures in the $J/\psi p$ and $J/\psi p\pi$
spectra in both exclusive and semi-inclusive processes will provide a smoking-
gun evidence for the hadronic molecule nature of those observed pentaquarks
and clarify the role played by the TS and BS in the near-threshold dynamics.

###### pacs:

11.55.Fv, 12.40.Yx ,14.40.Gx

Introduction The observation of hidden charm pentaquarks by the LHCb
Collaboration in 2015 LHCb:2015yax  provided a strong evidence for QCD exotics
in the baryon sector and has initiated dramatic activities in both experiment
and theory to understand their nature (see review articles Zhao:2016akg ;
Guo:2017jvc ; Chen:2016qju ; Olsen:2017bmm ; Esposito:2016noz ; Lebed:2016hpi
; Brambilla:2019esw ; Liu:2019zoy ; Ali:2017jda ; Karliner:2017qhf ;
Chen:2022asf  for details). Based on the Run-1 data sample, two states,
$P_{c}(4380)$ and $P_{c}(4450)$, were observed in the $J/\psi p$ invariant
mass spectrum in $\Lambda_{b}^{0}\to J/\psi pK^{-}$ LHCb:2015yax . Since their
masses were close to the threshold of $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ channel,
they were interpreted as hadronic molecules formed by the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ near-threshold interactions Chen:2015loa ;
Chen:2015moa ; Karliner:2015ina ; Roca:2015dva ; He:2015cea ; Lu:2016nnt ;
Shen:2016tzq ; Yamaguchi:2016ote ; Lin:2017mtz ; Yamaguchi:2017zmn ;
Wang:2018waa . With the Run-2 data samples, the reanalysis of
$\Lambda_{b}^{0}\to J/\psi pK^{-}$ exposes more detailed structures in the
$J/\psi p$ invariant mass spectrum LHCb:2019kea . $P_{c}(4450)$ splits into
two states, $P_{c}(4440)$ and $P_{c}(4457)$, which fit well the
$J^{P}=1/2^{-}$ and $3/2^{-}$ doublet formed by $\Sigma_{c}\bar{D}^{*}$, and a
new state $P_{c}(4312)$, which fits the $\Sigma_{c}\bar{D}$ molecular state,
is identified Chen:2019bip ; Chen:2019asm ; Liu:2019tjn ; He:2019ify ;
Guo:2019kdc ; Xiao:2019mvs ; Xiao:2019aya ; Meng:2019ilv ; Yamaguchi:2019seo ;
Du:2019pij ; Wang:2019ato . Meanwhile, signals for $P_{c}(4380)$, formerly
assigned as the $\Sigma_{c}^{*}\bar{D}$ molecule in the Run-1 analysis, become
obscure in the Run-2 analysis. It turns out to be broad and can be even
absorbed by the background LHCb:2019kea . It should be noted that the quantum
numbers determined by the Run-2 analysis are more consistent with the
molecular scenario based on the relative $S$-wave interactions and also
consistent with the early theoretical predictions Wu:2010vk ; Yang:2011wz ;
Wang:2011rga ; Wu:2012md .

A key issue clarified by the Run-2 analysis is that these hidden charmonium
pentaquark candidates couple to the nearby $\Sigma_{c}\bar{D}^{(*)}$ in a
relative $S$-wave. This makes the hadronic molecule scenario a natural
interpretation for their internal structures and provide important insights
into the threshold dynamics. Note that these observed $P_{c}$ states are
correlated with the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ thresholds instead of
$\Lambda_{c}^{(*)}\bar{D}^{(*)}$. It is also interesting to observe that this
feature seems to be consistent with the one boson exchange (OBE) dynamics
between the charmed meson and baryon. Namely, it manifests to be attractive
between $\Sigma_{c}^{(*)}$ and $\bar{D}^{(*)}$, but appears to be repulsive
between $\Lambda_{c}^{(*)}$ and $\bar{D}^{(*)}$ Wang:2011rga ; Yang:2011wz ;
Chen:2019bip ; Chen:2019asm ; Liu:2019tjn ; He:2019ify ; Guo:2019kdc ;
Xiao:2019mvs ; Xiao:2019aya ; Meng:2019ilv ; Yamaguchi:2019seo ; Du:2019pij ;
Wang:2019ato . This can explain that no obvious enhancements have been seen at
the thresholds of $\Lambda_{c}\bar{D}$ and $\Lambda_{c}\bar{D}^{*}$ in
exclusive decays of $\Lambda_{b}^{0}\to J/\psi pK^{-}$ LHCb:2019kea .

Another important threshold dynamics to be exposed by the pentaquark
production in the molecular picture is that the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ configurations in the $\Lambda_{b}$ decays
would be either produced via the color-suppressed process, or via their
subleading couplings to $\Lambda_{c}^{(*)}\bar{D}^{(*)}$ through the triangle
singularity (TS) Landau:1959fi ; Cutkosky:1960sp  enhanced transitions as
pointed out in Refs. Liu:2015fea ; Guo:2015umn . Notice that the
$\Lambda_{c}\bar{D}^{(*)}$ configurations can be accessed by both the direct
emission processes and TS enhanced transitions. The absence of enhancements
around the thresholds of $[\Lambda_{c}\bar{D}](4150)$ and
$[\Lambda_{c}\bar{D}^{*}](4290)$ seems to confirm that those observed hadronic
molecules states are genuine states driven by the
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ interactions, while the repulsive nature of
$\Lambda_{c}^{(*)}\bar{D}^{(*)}$ would not produce apparent threshold
structures at leading order in the exclusive decays of $\Lambda_{b}$. In Ref.
Wang:2015jsa  it was proposed to search for these heavy pentaquark states in
$J/\psi$ photoproduction where the cross sections at large scattering angles
should serve as direct probe for the production of the pentaquarks as genuine
states. Although these hidden charm pentaquark states are searched at GlueX
based on their 12 GeV data, the present statistics do not allow a conclusion
for their existence at this moment GlueX:2019mkq .

Apparently, the confirmation or refusal of $P_{c}(4380)$ will have a strong
impact on our understanding of the nature of these hidden charm pentaquarks.
In this work, based on the hadronic molecule picture, we identify the crucial
role played by the TS and BS in the pentaquark three-body and four-body decays
into $J/\psi p\pi$ and $J/\psi p\pi\pi$, respectively. It will strongly
enhance the width of $P_{c}(4380)$ due to the special kinematic in
$P_{c}(4380)\to J/\psi p\pi$, and naturally accommodate it into the hadronic
molecule picture. As a smoking-gun evidence for this mechanism, we predict two
feed-down enhancements at the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ thresholds in the $J/\psi p$ and $J/\psi p\pi$
invariant mass spectra, respectively. We will discuss that high-statistics
data for exclusive processes or semi-inclusive channels of $J/\psi\pi$ and
$J/\psi\pi\pi$ will be ideal for the search for these two threshold
structures.

Formalism Given the productions of these $P_{c}$ states via various processes,
their decays into $J/\psi p$ can be via a direct two-body decays into $J/\psi
p$, three-body decays into $J/\psi p\pi$, or four-body decays into $J/\psi
p\pi\pi$. An illustration of these transitions is presented in Fig. 1. Note
that in both Fig. 1 (b) and (c) the transition satisfy the TS condition, i.e.
all the internal particles can approach their on-shell condition
simultaneously. As the leading singularity in the loop functions, it may
provide non-negligible partial widths for these $P_{c}$ states via their
$J/\psi p\pi$ three-body decays. In comparison with the two-body decay
amplitude which counts an order of $v$, where $v$ is the typical non-
relativistic on-shell velocity of the intermediate heavy hadrons, the triangle
loop of Fig. 1 (a) and (b) count $(v^{5}/v^{6})v\sim{\cal O}(1)$.

The transitions of Fig. 1 (d) fulfill the so-called “box singularity” (BS),
where all the internal particles can also approach their on-shell conditions
simultaneously. The amplitudes of Fig. 1 (d) count $(v^{5}/v^{8})v^{2}\sim
1/v$ which appears to be more leading than (b) and (c). However, taking into
account the limited phase space, contributions from $P_{c}(4440)$ and
$P_{c}(4457)\to J/\psi p\pi\pi$ will be suppressed compared to those three-
body decays.

In the feed-down scenario the semi-inclusive spectra of $J/\psi p$ or $J/\psi
p\pi$ will be measured. It means that all the transitions in Fig. 1 will
contribute and sum from the $J/\psi p$ or $J/\psi p\pi$ threshold to higher
energies. Due to the triangle singularity, we would expect the threshold
enhancement at the mass of $\Lambda_{c}\bar{D}$ in the $J/\psi p$ invariant
mass spectrum, and due to the box singularity, we would expect the threshold
enhancements at the thresholds of $\Lambda_{c}\bar{D}^{*}$ and
$\Sigma_{c}\bar{D}$ in the $J/\psi p\pi$ invariant mass spectrum.

Figure 1: Feynman diagrams illustrating the $P_{c}$ states decay into $J/\psi
p$, $J/\psi p\pi$, and $J/\psi p\pi\pi$. Figures (b)-(d) manifest the TS and
BS mechanisms that can feed down to the semi-inclusive $J/\psi p$ and $J/\psi
p\pi$ spectra, respectively.

To proceed, we employ an effective Lagrangian approach to evaluate the feed-
down mechanisms. The $S$-wave interactions between $P_{c}$ and
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$ can be expressed as Zou:2002yy

$\begin{split}\mathcal{L}_{P_{c}(4312)\Sigma_{c}\bar{D}}&=g_{P_{c}(4312)}\bar{D}^{\dagger}\bar{\Sigma}_{c}P_{c},\\\
\mathcal{L}_{P_{c}(4380)\Sigma_{c}^{*}\bar{D}}&=g_{P_{c}(4380)}\bar{D}^{\dagger}\bar{\Sigma}^{*}_{c}\cdot
P_{c},\\\
\mathcal{L}_{P_{c}(4440)\Sigma_{c}\bar{D}^{*}}&=g_{P_{c}(4440)}\bar{D}^{*\mu{\dagger}}\bar{\Sigma}_{c}\gamma_{5}\tilde{\gamma}_{\mu}(p_{A})P_{c},\\\
\mathcal{L}_{P_{c}(4457)\Sigma_{c}\bar{D}^{*}}&=g_{P_{c}(4457)}\bar{\Sigma}_{c}\bar{D}^{*{\dagger}}\cdot
P_{c},\\\ \end{split}$ (1)

where $\bar{\Sigma}_{c}^{(*)}$ and $\bar{D}^{(*){\dagger}}$ denote the field
operators for $\Sigma_{c}^{*}$ and $\bar{D}^{(*)}$, respectively, and $P_{c}$
denotes the spinor field of $P_{c}$ states; $g_{P_{c}}$ are the coupling
constants between $P_{c}$ and corresponding $\Sigma_{c}^{(*)}\bar{D}^{(*)}$
channels. They are determined in the molecular scenario by reproducing the
mass and width of the $P_{c}$ states. We adopt the coupling constants
$g_{P_{c}(4312)}=2.16$, $g_{P_{c}(4380)}=2.19$, $g_{P_{c}(4440)}=2.6$, and
$g_{P_{c}(4457)}=1.7$ determined from Ref. Xie:2022hhv .
$\tilde{\gamma}_{\mu}(p_{A})$ is a momentum dependant function, which is
defined as
$\tilde{\gamma}^{\mu}(p_{A})=(-g^{\mu\nu}+p_{A}^{\mu}p_{A}^{\nu}/p_{A}^{2})\gamma_{\nu}$
Zou:2002yy . The projection operator for the spin-3/2 particle can be written
as

$\displaystyle\sum u^{\nu}\bar{u}^{\alpha}$ $\displaystyle=$
$\displaystyle(\not{p}+m)[-g^{\nu\alpha}+\frac{1}{3}\gamma^{\nu}\gamma^{\alpha}$
(2)
$\displaystyle+\frac{1}{3m}(\gamma^{\nu}p^{\alpha}-\gamma^{\alpha}p^{\nu})+\frac{2}{3m^{2}}p^{\nu}p^{\alpha}],$

with which the propagator of $\Sigma_{c}^{*}$ in Fig. 1 (b) can be obtained.
Lagrangians for other vertices in Fig. 1 are collected as follows:

$\displaystyle\mathcal{L}_{\Lambda_{c}\bar{D}J/\psi p}$ $\displaystyle=$
$\displaystyle
ig_{x}\psi^{\mu}\bar{N}\gamma_{5}\gamma_{\mu}\Lambda_{c}\bar{D},$
$\displaystyle\mathcal{L}_{\Sigma_{c}\Lambda_{c}\pi}$ $\displaystyle=$
$\displaystyle
g_{\Sigma_{c}\Lambda_{c}\pi}\bar{\Lambda}_{c}\partial^{\mu}\pi\gamma_{5}\gamma_{\mu}\Sigma_{c},$
$\displaystyle\mathcal{L}_{\Sigma^{*}_{c}\Lambda_{c}\pi}$ $\displaystyle=$
$\displaystyle
g_{\Sigma_{c}^{*}\Lambda_{c}\pi}\bar{\Lambda}_{c}\partial^{\mu}\pi\Sigma_{c\mu}^{*},$
$\displaystyle\mathcal{L}_{D^{*}D\pi}$ $\displaystyle=$ $\displaystyle
ig_{D^{*}D\pi}D_{\mu}^{*}\bar{D}\partial^{\mu}\pi\ ,$ (3)

where the coupling constants $g_{\Sigma_{c}\Lambda_{c}\pi}$,
$g_{\Sigma_{c}^{*}\Lambda_{c}\pi}$, and $g_{D^{*}D\pi}$ are determined by the
partial decay widths of the corresponding processes ParticleDataGroup:2022pth
. Coupling $g_{x}$ describes the contact interaction for the
$\Lambda_{c}\bar{D}J/\psi p$ vertex. It can be estimated by a $t$-channel $D$
meson exchange near threshold as shown in Fig. 2.

Figure 2: The near-threshold interaction between $\Lambda_{c}\bar{D}$ and
$J/\psi p$ is described by a $t$-channel $D$ meson exchange from which the
effective contact coupling $g_{x}$ is extracted.

The Lagrangians involved in the determination of $g_{x}$ in Fig. 2 (a) are

$\begin{split}\mathcal{L}_{\Lambda_{c}Dp}&=ig_{DN\Lambda_{c}}D^{\dagger}\bar{N}\gamma_{5}\Lambda_{c},\\\
\mathcal{L}_{J/\psi D\bar{D}}&=g_{J/\psi
D\bar{D}}(D\partial_{\mu}\bar{D}-\partial_{\mu}D\bar{D})\psi^{\mu},\\\
\end{split}$ (4)

where the coupling constant $g_{DN\Lambda_{c}}=13.2$ is estimated the same as
the coupling of $KN\Lambda$ by employing the $SU(4)$ symmetry Rijken:1998yy ;
Stoks:1999bz ; Oh:2006hm ; Liu:2010um ; Coupling $g_{J/\psi
D\bar{D}}=2g_{2}\sqrt{m_{J/\psi}m_{D}m_{\bar{D}}}$ with
$g_{2}=\sqrt{m_{\psi}}/(2m_{D}f_{J/\psi})$ and $f_{J/\psi}=0.405~{}{\rm GeV}$
Colangelo:2003sa  is adopted from the heavy quark symmetry relation. As
illustrated in Fig. 2, the relationship
$|\mathcal{M}^{(a)}|^{2}=|\mathcal{M}^{(b)}|^{2}$ is employed, which has a
concrete form like:

$\displaystyle g_{x}^{2}$ $\displaystyle=$
$\displaystyle\frac{g_{DN\Lambda_{c}}^{2}g_{J/\psi
D\bar{D}}^{2}|\bar{u}_{p}\gamma_{5}u_{\Lambda_{c}}(p_{\bar{D}}-q_{D})_{\mu}\frac{1}{q_{D}^{2}-m_{D}^{2}}\varepsilon^{*\mu}_{J/\psi}|^{2}}{|\bar{u}_{p}\gamma_{5}\gamma_{\mu}u_{\Lambda_{c}}|^{2}},$
(5)

which yields $|g_{x}|=3.0~{}{\rm GeV}^{-1}$ at the threshold of $\Lambda_{c}$
and $\bar{D}$.

With the above Lagrangians, the triangle and box loop amplitudes in Fig. 1 can
be expressed as

$\begin{split}&\mathcal{M}^{P_{c}(4312)\to J/\psi p\pi}\\\
&=\int\frac{dp^{4}_{\Lambda_{c}}}{(2\pi)^{4}}\bar{u}_{p}\gamma_{5}\gamma_{\mu}\varepsilon_{J/\psi}^{*\mu}(\not{p}_{\Lambda_{c}}+m_{\Lambda_{c}})\gamma_{5}\gamma_{\nu}p_{\pi}^{\mu}(\not{p}_{\Sigma_{c}}+m_{\Sigma_{c}})u_{P_{c}}\\\
&\quad\times\frac{g_{P_{c}(4312)}g_{\Sigma_{c}\Lambda_{c}\pi}g_{x}}{(p_{\Sigma_{c}}^{2}-m_{\Sigma_{c}}^{2})(p_{D}^{2}-m_{D}^{2})(p_{\Lambda_{c}}^{2}-m_{\Lambda_{c}}^{2})},\\\
&\mathcal{M}^{P_{c}(4380)\to J/\psi p\pi}\\\
&=\int\frac{d^{4}p_{\Lambda_{c}}}{(2\pi)^{4}}\bar{u}_{p}\gamma_{5}\gamma_{\mu}\varepsilon_{J/\psi}^{*\mu}(\not{p}_{\Lambda_{c}}+m_{\Lambda_{c}})p_{\pi\nu}(\not{p}_{\Sigma_{c}}+m_{\Sigma_{c}})\\\
&\quad\times[-g^{\nu\alpha}+\frac{1}{3}\gamma^{\nu}\gamma^{\alpha}+\frac{1}{3m_{\Sigma_{c}}}(\gamma^{\nu}p^{\alpha}-\gamma^{\alpha}p^{\nu})+\frac{2}{3p_{\Sigma_{c}}^{2}}p_{\Sigma_{c}}^{\nu}p_{\Sigma_{c}}^{\alpha}]\\\
&\quad\times
u_{P_{c}\alpha}\frac{g_{P_{c}(4380)}g_{\Sigma_{c}^{*}\Lambda_{c}\pi}g_{x}}{(p_{\Sigma_{c}^{*}}^{2}-m_{\Sigma_{c}^{*}}^{2})(p_{D}^{2}-m_{D}^{2})(p_{\Lambda_{c}}^{2}-m_{\Lambda_{c}}^{2})},\\\
&\mathcal{M}^{P_{c}(4440)\to J/\psi p\pi\pi}\\\
&=\int\frac{d^{4}p_{D}}{(2\pi)^{4}}(-g^{\alpha\rho}+\frac{p_{D^{*}}^{\alpha}p_{D^{*}}^{\rho}}{m_{D^{*}}^{2}})p_{\pi\rho}\\\
&\quad\times\bar{u}_{p}\gamma_{5}\gamma_{\mu}\varepsilon_{J/\psi}^{*\mu}(\not{p}_{\Lambda_{c}}+m_{\Lambda_{c}})\gamma_{5}\gamma_{\nu}p_{\pi}^{\nu}(\not{p}_{\Sigma_{c}}+m_{\Sigma_{c}})\gamma_{5}\tilde{\gamma}_{\alpha}u_{P_{c}}\\\
&\quad\times\frac{g_{P_{c}(4440)}g_{\Sigma_{c}^{*}\Lambda_{c}\pi}g_{x}g_{D^{*}D\pi}}{(p_{\Sigma_{c}}^{2}-m_{\Sigma_{c}}^{2})(p_{D^{*}}^{2}-m_{D^{*}}^{2})(p_{\Lambda_{c}}^{2}-m_{\Lambda_{c}}^{2})(p_{D}^{2}-m_{D}^{2})},\\\
&\mathcal{M}^{P_{c}(4457)\to J/\psi p\pi\pi}\\\
&=\int\frac{d^{4}p_{D}}{(2\pi)^{4}}\bar{u}_{p}\gamma_{5}\gamma_{\mu}\varepsilon_{J/\psi}^{*\mu}(\not{p}_{\Lambda_{c}}+m_{\Lambda_{c}})\gamma_{5}\gamma_{\nu}p_{\pi}^{\nu}(\not{p}_{\Sigma_{c}}+m_{\Sigma_{c}})u_{P_{c}\alpha}\\\
&\quad\times(-g^{\alpha\beta}+\frac{p_{D^{*}}^{\alpha}p_{D^{*}}^{\beta}}{m_{D^{*}}^{2}})p_{\pi\beta}\\\
&\quad\times\frac{g_{P_{c}(4457)}g_{\Sigma_{c}^{*}\Lambda_{c}\pi}g_{x}g_{D^{*}D\pi}}{(p_{\Sigma_{c}}^{2}-m_{\Sigma_{c}}^{2})(p_{D^{*}}^{2}-m_{D^{*}}^{2})(p_{\Lambda_{c}}^{2}-m_{\Lambda_{c}}^{2})(p_{D}^{2}-m_{D}^{2})}.\\\
\end{split}$ (6)

In the kinematic regions of the TS and BS, the internal particles are in the
vicinity of their on-shell conditions. We further simplify the integrands by
adopting the following eigenfunction relations for the spin-$1/2$, spin-$1$,
and spin-$3/2$ particles can be simplified as $(\not{p}+m)\to 2m$,
$(-g^{\mu\nu}+p^{\mu}p^{\nu}/p^{2})\to\delta^{ij}$, and
$(\not{p}+m)[-g^{\mu\nu}+\frac{1}{3}\gamma^{\mu}\gamma^{\nu}+\frac{1}{3m}(\gamma^{\mu}p^{\nu}-\gamma^{\nu}p^{\mu})+\frac{2}{3m^{2}}p^{\mu}p^{\nu}]\to
2m[-g^{\mu\nu}+\frac{1}{3}\gamma^{\mu}\gamma^{\nu}+\frac{1}{3}(\gamma^{\mu}g^{\nu
0}-\gamma^{\nu}g^{\mu 0})+\frac{2}{3}g^{\mu 0}g^{\nu 0}]$. As a result, the
numerators in the integrands do not contain the inner momentum and the loop
integrals are actually convergent.

For the relatively broad states $\Sigma_{c}^{(*)}$, complex masses with the
widths are adopted, while for $D^{*}$ we have ignored its width effects since
they do not cause significant changes to the lineshapes of the $J/\psi p$ or
$J/\psi p\pi$ invariant mass spectra.

Results and discussions Proceeding to the numerical calculations we first list
in Table. 1 the masses and widths adopted for these observed $P_{c}$ states
LHCb:2019kea ; ParticleDataGroup:2022pth . We stress that it is the TS and BS
mechanisms driving the enhancements expected at the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ thresholds. These two threshold structures are
insensitive to the uncertainties of the masses and widths of the $P_{c}$
states, and no need for us to re-derive these quantities here.

Table 1: Masses (unit: MeV) adopted for the $P_{c}$ states.  $m_{P_{c}(4312)}=4312$ | $m_{P_{c}(4380)}=4380$ | $m_{P_{c}(4440)}=4440$
---|---|---
$m_{P_{c}(4457)}=4457$ | $m_{\Sigma_{c}}=2453-i1.89/2$ | $m_{\Sigma_{c}^{*}}=2518-i15/2$
$m_{\Lambda_{c}}=2286$ | $m_{p}=938$ | $m_{J/\psi}=3097$
$m_{\pi}=137$ | $m_{D}=1867$ | $m_{D^{*}}=2009$
Figure 3: The $J/\psi p$ invariant mass spectra (red solid lines) from
different $P_{c}$ state decays. The dashed vertical lines indicate the
$\Lambda_{c}\bar{D}$ threshold, and the purple bands denote the range of
$4.13{\rm GeV}<m_{J/\psi p}<4.17{\rm GeV}$ the invariant mass spectra, which
are collected in Fig. 5.

For the feed-down processes, it is not necessary to have the initial states to
be the pentaquark states. In principle, as far as the initial
$\Sigma_{c}^{(*)}$ and $\bar{D}^{(*)}$ scatterings occur in an $S$-wave near
threshold, they will feed down to contribute to the production of the
threshold enhancements at the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ thresholds in the $J/\psi p$ and $J/\psi p\pi$
invariant mass spectra, respectively, via the TS or BS mechanisms. If the
initial states are genuine states, the enhancements will be strongly enhanced
by the propagators when they are on-shell. Thus, even though we do the
calculations with the $P_{c}$ states as the initial states, the feed-down
mechanisms actually involve both possible genuine states and continuum
contributions. This unique feature makes the semi-inclusive measurements of
the $J/\psi p$ and $J/\psi p\pi$ spectra interesting and useful for clarifying
not only the nature of the $\Lambda_{c}\bar{D}^{(*)}$ interactions, but also
the role played by the TS and BS mechanisms.

In Fig. 3 we plot the $J/\psi p$ invariant mass spectra for four observed
$P_{c}$ states in their three-body or four-body decays into either $J/\psi
p\pi$ or $J/\psi p\pi\pi$. Clear and predominant peaks appear at the
$\Lambda_{c}\bar{D}$ threshold due to the TS and BS mechanisms. Meanwhile, it
shows that $P_{c}(4380)$ contributes the largest to the $J/\psi p$ spectrum.
It is due to the relatively large phase space for
$\Sigma_{c}^{*}\to\Lambda_{c}\pi$ within the TS kinematic region that strongly
enhances the decay rate of $P_{c}(4380)\to J/\psi p\pi$. Such an effect
actually will add additional partial decay widths to the $P_{c}(4380)$ decay
111We note that the TS also occurs in another situation in the three-body and
four-body decays into $\Lambda_{c}\bar{D}\pi$ and $\Lambda_{c}\bar{D}\pi\pi$,
respectively. Different from Fig. 1, here the significant contributions from
the $S$-wave coupling of $\bar{D}\pi$ into $\bar{D}_{0}(2400)$. These
processes will also contribute to the total width of $P_{c}(4380)$ although it
cannot be as large as 100 MeV suggested by the LHCb fits LHCb:2019kea . More
detailed studies of such processes will be reported elsewhere..

Figure 4: The $J/\psi p\pi$ invariant mass spectra (red solid lines) from
$P_{c}(4440)$ (left) and $P_{c}(4457)$ (right) decays. The
$\Lambda_{c}\bar{D}^{*}$ threshold is indicated by the dashed vertical lines.

The semi-inclusive $J/\psi p$ spectrum should contain contributions from the
exclusive decay channels such as the $P_{c}$ discovery channel $\Lambda_{b}\to
K^{-}J/\psi p$ LHCb:2015yax ; LHCb:2019kea . As shown in Ref. Liu:2015fea ,
the threshold of $\Lambda_{c}\bar{D}$ will appear as a two-body CUSP because
the mass of $\Lambda_{b}$ is far above the thresholds of
$\Lambda_{c}\bar{D}_{s}^{*}$ and $\Lambda_{c}\bar{D}_{s1}$ and the TS
condition cannot be satisfied. This situation actually helps us clarify the
role played by the TS and BS mechanisms since a narrow peak at the threshold
of $\Lambda_{c}\bar{D}$ in the semi-inclusive $J/\psi p$ spectrum will
indicate the leading TS and/or BS mechanisms instead of a CUSP effect
Guo:2014iya .

The feed-down process via the BS mechanism also gives access to the
$\Lambda_{c}\bar{D}^{*}$ threshold which can manifest as a narrow threshold
enhancement as shown by the invariant mass spectrum of $J/\psi p\pi$ in Fig.
4. We stress that although such a threshold enhancement is illustrated by the
four-body decays of $P_{c}(4440)$ and $P_{c}(4457)$, any processes involving
the near-threshold $\Sigma_{c}^{(*)}\bar{D}^{*}$ $S$-wave scatterings may feed
down to the $J/\psi p\pi\pi$ channel and produce the $\Lambda_{c}\bar{D}^{*}$
threshold enhancement. This phenomenon can also serve as a unique probe for
the role played by the BS mechanism and can be searched in experiment.

In Fig. 1 (c) in addition to the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ in the invariant mass spectra of $J/\psi p$ and
$J/\psi p\pi$, respectively, the $\Sigma_{c}^{*}\bar{D}$ threshold can also be
accessed due to the BS mechanism. However, since the $\Sigma_{c}^{*}\bar{D}$
and $J/\psi p\pi$ thresholds are close to each other the
$\Sigma_{c}^{*}\bar{D}$ threshold structure will be strongly suppressed.

It should be clarified that, although it is possible to observe the
$\Lambda_{c}\bar{D}$ and $\Lambda_{c}\bar{D}^{*}$ threshold structures in
exclusive processes, e.g. in $\Lambda_{b}\to K^{-}J/\psi p$ and
$\Lambda_{b}\to K^{-}J/\psi p\pi$, the signal significance should be less than
those of the $P_{c}$ states. The reason is that the relative $S$-wave
interactions between $\Lambda_{c}$ and $\bar{D}^{(*)}$ are repulsive in the
hadronic molecular picture. Without the pole structures, the signals
accumulated from the $P_{c}$ decays due to the feed-down mechanism will be
limited.

In order to demonstrate the feed-down effects we sum up the feed-down
contributions from the $P_{c}$ decays and plot the $J/\psi p$ and $J/\psi
p\pi$ spectra in Fig. 5 (a) and (b), respectively. We neglect the possible
interferences at this moment. They may occur if the feed-down contributions
are from the same exclusive process. One can see that as the lowest quark-
rearranging threshold above the $J/\psi p$ channel, the TS and BS kinematics
turn out to be unique. The signals originated from the TS and BS mechanisms
will drop smoothly in the higher energy end. It means that these feed-down
signals will stand predominantly on the incoherent background in the semi-
inclusive measurement.

Figure 5: The accumulated $J/\psi p$ (left) and $J/\psi p\pi$ (right)
invariant mass spectra from all different $P_{c}$ states are shown by the
black solid lines. On the left panel the red dashed and blue dot dashed lines
represent the results from $P_{c}(4380)$ and the sum of all the others,
respectively. On the right panel the red dashed and blue dot dashed lines
represent the results from $P_{c}(4457)$ and $P_{c}(4440)$, respectively.

Summary We study the feed-down transitions of the observed $P_{c}$ states into
$J/\psi p$ and $J/\psi p\pi$ via rescatterings, and predict that predominant
narrow structures can be produced at the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ thresholds due to the TS and BS mechanisms. Since the
$S$-wave interactions between $\Lambda_{c}$ and $\bar{D}^{(*)}$ are repulsive
in the hadronic molecule picture, a search for the $\Lambda_{c}\bar{D}$ and
$\Lambda_{c}\bar{D}^{*}$ threshold enhancements can provide supportive
evidences for the hadronic molecule interpretation for these $P_{c}$ states
including $P_{c}(4380)$. Meanwhile, it will also provide a unique way for
pinning down the role played by the TS and BS mechanisms. Since such feed-down
transitions, in principle, can occur given the production of near-threshold
$\Sigma_{c}^{(*)}\bar{D}^{(*)}$, we anticipate that the $\Lambda_{c}\bar{D}$
and $\Lambda_{c}\bar{D}^{*}$ threshold enhancements can be accessed in the
semi-inclusive $J/\psi p$ and $J/\psi p\pi$ spectra in various processes, such
as hadron-hadron collisions at LHCb, $e^{+}e^{-}$ annihilations at Belle-II,
and $J/\psi$ photoproduction reactions at GlueX. Future experimental search
for them are strongly recommended.

Acknowledgement This work is supported, in part, by the National Natural
Science Foundation of China (Grant No. 12235018), DFG and NSFC funds to the
Sino-German CRC 110 “Symmetries and the Emergence of Structure in QCD” (NSFC
Grant No. 12070131001, DFG Project-ID 196253076), National Key Basic Research
Program of China under Contract No. 2020YFA0406300, and Strategic Priority
Research Program of Chinese Academy of Sciences (Grant No. XDB34030302).

## References

  * (1) R. Aaij et al. [LHCb], Observation of $J/\psi p$ Resonances Consistent with Pentaquark States in $\Lambda_{b}^{0}\to J/\psi K^{-}p$ Decays, Phys. Rev. Lett. 115, 072001 (2015).
  * (2) Q. Zhao, The Excitement and Challenging Issues Concerning the Discovery of Heavy Pentaquark Candidates with H, AAPPS Bull. 26, no.4, 8-14 (2016).
  * (3) F. K. Guo, C. Hanhart, U. G. Meißner, Q. Wang, Q. Zhao and B. S. Zou, Hadronic molecules, Rev. Mod. Phys. 90, no.1, 015004 (2018), [erratum: Rev. Mod. Phys. 94, no.2, 029901 (2022)].
  * (4) H. X. Chen, W. Chen, X. Liu and S. L. Zhu, The hidden-charm pentaquark and tetraquark states, Phys. Rept. 639, 1-121 (2016).
  * (5) S. L. Olsen, T. Skwarnicki and D. Zieminska, Rev. Mod. Phys. 90, no.1, 015003 (2018) doi:10.1103/RevModPhys.90.015003 [arXiv:1708.04012 [hep-ph]].
  * (6) A. Esposito, A. Pilloni and A. D. Polosa, Phys. Rept. 668, 1-97 (2017) doi:10.1016/j.physrep.2016.11.002 [arXiv:1611.07920 [hep-ph]].
  * (7) R. F. Lebed, R. E. Mitchell and E. S. Swanson, Prog. Part. Nucl. Phys. 93, 143-194 (2017) doi:10.1016/j.ppnp.2016.11.003 [arXiv:1610.04528 [hep-ph]].
  * (8) N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C. P. Shen, C. E. Thomas, A. Vairo and C. Z. Yuan, Phys. Rept. 873, 1-154 (2020) doi:10.1016/j.physrep.2020.05.001 [arXiv:1907.07583 [hep-ex]].
  * (9) Y. R. Liu, H. X. Chen, W. Chen, X. Liu and S. L. Zhu, Prog. Part. Nucl. Phys. 107, 237-320 (2019) doi:10.1016/j.ppnp.2019.04.003 [arXiv:1903.11976 [hep-ph]].
  * (10) A. Ali, J. S. Lange and S. Stone, Prog. Part. Nucl. Phys. 97, 123-198 (2017) doi:10.1016/j.ppnp.2017.08.003 [arXiv:1706.00610 [hep-ph]].
  * (11) M. Karliner, J. L. Rosner and T. Skwarnicki, Ann. Rev. Nucl. Part. Sci. 68, 17-44 (2018) doi:10.1146/annurev-nucl-101917-020902 [arXiv:1711.10626 [hep-ph]].
  * (12) H. X. Chen, W. Chen, X. Liu, Y. R. Liu and S. L. Zhu, Rept. Prog. Phys. 86, no.2, 026201 (2023) doi:10.1088/1361-6633/aca3b6 [arXiv:2204.02649 [hep-ph]].
  * (13) R. Chen, X. Liu, X. Q. Li and S. L. Zhu, Identifying exotic hidden-charm pentaquarks, Phys. Rev. Lett. 115, no.13, 132002 (2015).
  * (14) H. X. Chen, W. Chen, X. Liu, T. G. Steele and S. L. Zhu, Towards exotic hidden-charm pentaquarks in QCD, Phys. Rev. Lett. 115, no.17, 172001 (2015).
  * (15) M. Karliner and J. L. Rosner, New Exotic Meson and Baryon Resonances from Doubly-Heavy Hadronic Molecules, Phys. Rev. Lett. 115, no.12, 122001 (2015).
  * (16) L. Roca, J. Nieves and E. Oset, LHCb pentaquark as a $\bar{D}^{*}\Sigma_{c}-\bar{D}^{*}\Sigma_{c}^{*}$ molecular state, Phys. Rev. D 92, no.9, 094003 (2015).
  * (17) J. He, $\bar{D}\Sigma^{*}_{c}$ and $\bar{D}^{*}\Sigma_{c}$ interactions and the LHCb hidden-charmed pentaquarks, Phys. Lett. B 753, 547-551 (2016).
  * (18) Q. F. Lü and Y. B. Dong, Strong decay mode $J/\psi p$ of hidden charm pentaquark states $P_{c}^{+}(4380)$ and $P_{c}^{+}(4450)$ in $\Sigma_{c}\bar{D}^{*}$ molecular scenario, Phys. Rev. D 93, no.7, 074020 (2016).
  * (19) C. W. Shen, F. K. Guo, J. J. Xie and B. S. Zou, Disentangling the hadronic molecule nature of the $P_{c}(4380)$ pentaquark-like structure, Nucl. Phys. A 954, 393-405 (2016).
  * (20) Y. Yamaguchi and E. Santopinto, Hidden-charm pentaquarks as a meson-baryon molecule with coupled channels for $\bar{D}^{(\ast)}\Lambda_{\rm c}$ and $\bar{D}^{(\ast)}\Sigma^{(\ast)}_{\rm c}$, Phys. Rev. D 96, no.1, 014018 (2017).
  * (21) Y. H. Lin, C. W. Shen, F. K. Guo and B. S. Zou, Decay behaviors of the $P_{c}$ hadronic molecules, Phys. Rev. D 95, no.11, 114017 (2017).
  * (22) Y. Yamaguchi, A. Giachino, A. Hosaka, E. Santopinto, S. Takeuchi and M. Takizawa, Hidden-charm and bottom meson-baryon molecules coupled with five-quark states, Phys. Rev. D 96, no.11, 114031 (2017).
  * (23) Z. G. Wang, Analysis of the $\bar{D}\Sigma_{c}$, $\bar{D}\Sigma_{c}^{*}$, $\bar{D}^{*}\Sigma_{c}$ and $\bar{D}^{*}\Sigma_{c}^{*}$ pentaquark molecular states with QCD sum rules, Int. J. Mod. Phys. A 34, no.19, 1950097 (2019).
  * (24) R. Aaij et al. [LHCb], Observation of a narrow pentaquark state, $P_{c}(4312)^{+}$, and of two-peak structure of the $P_{c}(4450)^{+}$, Phys. Rev. Lett. 122, no.22, 222001 (2019).
  * (25) H. X. Chen, W. Chen and S. L. Zhu, Possible interpretations of the $P_{c}(4312)$, $P_{c}(4440)$, and $P_{c}(4457)$, Phys. Rev. D 100, no.5, 051501 (2019).
  * (26) R. Chen, Z. F. Sun, X. Liu and S. L. Zhu, Strong LHCb evidence supporting the existence of the hidden-charm molecular pentaquarks, Phys. Rev. D 100, no.1, 011502 (2019).
  * (27) M. Z. Liu, Y. W. Pan, F. Z. Peng, M. Sánchez Sánchez, L. S. Geng, A. Hosaka and M. Pavon Valderrama, Emergence of a complete heavy-quark spin symmetry multiplet: seven molecular pentaquarks in light of the latest LHCb analysis, Phys. Rev. Lett. 122, no.24, 242001 (2019).
  * (28) J. He, Study of $P_{c}(4457)$, $P_{c}(4440)$, and $P_{c}(4312)$ in a quasipotential Bethe-Salpeter equation approach, Eur. Phys. J. C 79, no.5, 393 (2019).
  * (29) Z. H. Guo and J. A. Oller, Anatomy of the newly observed hidden-charm pentaquark states: $P_{c}(4312)$, $P_{c}(4440)$ and $P_{c}(4457)$, Phys. Lett. B 793, 144-149 (2019).
  * (30) C. J. Xiao, Y. Huang, Y. B. Dong, L. S. Geng and D. Y. Chen, Exploring the molecular scenario of Pc(4312) , Pc(4440) , and Pc(4457), Phys. Rev. D 100, no.1, 014022 (2019).
  * (31) C. W. Xiao, J. Nieves and E. Oset, Heavy quark spin symmetric molecular states from ${\bar{D}}^{(*)}\Sigma_{c}^{(*)}$ and other coupled channels in the light of the recent LHCb pentaquarks, Phys. Rev. D 100, no.1, 014021 (2019).
  * (32) L. Meng, B. Wang, G. J. Wang and S. L. Zhu, The hidden charm pentaquark states and $\Sigma_{c}\bar{D}^{(*)}$ interaction in chiral perturbation theory, Phys. Rev. D 100, no.1, 014031 (2019).
  * (33) Y. Yamaguchi, H. García-Tecocoatzi, A. Giachino, A. Hosaka, E. Santopinto, S. Takeuchi and M. Takizawa, $P_{c}$ pentaquarks with chiral tensor and quark dynamics, Phys. Rev. D 101, no.9, 091502 (2020).
  * (34) M. L. Du, V. Baru, F. K. Guo, C. Hanhart, U. G. Meißner, J. A. Oller and Q. Wang, Interpretation of the LHCb $P_{c}$ States as Hadronic Molecules and Hints of a Narrow $P_{c}(4380)$, Phys. Rev. Lett. 124, no.7, 072001 (2020).
  * (35) B. Wang, L. Meng and S. L. Zhu, Hidden-charm and hidden-bottom molecular pentaquarks in chiral effective field theory, JHEP 11, 108 (2019).
  * (36) J. J. Wu, R. Molina, E. Oset and B. S. Zou, Dynamically generated $N^{*}$ and $\Lambda^{*}$ resonances in the hidden charm sector around 4.3 GeV, Phys. Rev. C 84, 015202 (2011).
  * (37) Z. C. Yang, Z. F. Sun, J. He, X. Liu and S. L. Zhu, The possible hidden-charm molecular baryons composed of anti-charmed meson and charmed baryon, Chin. Phys. C 36, 6-13 (2012).
  * (38) W. L. Wang, F. Huang, Z. Y. Zhang and B. S. Zou, $\Sigma_{c}\bar{D}$ and $\Lambda_{c}\bar{D}$ states in a chiral quark model, Phys. Rev. C 84, 015203 (2011).
  * (39) J. J. Wu, T. S. H. Lee and B. S. Zou, Nucleon Resonances with Hidden Charm in Coupled-Channel Models, Phys. Rev. C 85, 044002 (2012).
  * (40) X. H. Liu, Q. Wang and Q. Zhao, Understanding the newly observed heavy pentaquark candidates, Phys. Lett. B 757, 231-236 (2016).
  * (41) L. D. Landau, On analytic properties of vertex parts in quantum field theory, Nucl. Phys. 13, no.1, 181-192 (1959).
  * (42) R. E. Cutkosky, Singularities and discontinuities of Feynman amplitudes, J. Math. Phys. 1, 429-433 (1960).
  * (43) F. K. Guo, U. G. Meißner, W. Wang and Z. Yang, How to reveal the exotic nature of the Pc(4450), Phys. Rev. D 92, no.7, 071502 (2015).
  * (44) Q. Wang, X. H. Liu and Q. Zhao, Photoproduction of hidden charm pentaquark states $P_{c}^{+}(4380)$ and $P_{c}^{+}(4450)$, Phys. Rev. D 92, 034022 (2015).
  * (45) A. Ali et al. [GlueX], First Measurement of Near-Threshold J/$\psi$ Exclusive Photoproduction off the Proton, Phys. Rev. Lett. 123, no.7, 072001 (2019).
  * (46) B. S. Zou and F. Hussain, Covariant L-S scheme for the effective N*NM couplings, Phys. Rev. C 67, 015204 (2003).
  * (47) J. M. Xie, X. Z. Ling, M. Z. Liu and L. S. Geng, Search for hidden-charm pentaquark states in three-body final states, Eur. Phys. J. C 82, no.11, 1061 (2022).
  * (48) R. L. Workman et al. [Particle Data Group], Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * (49) T. A. Rijken, V. G. J. Stoks and Y. Yamamoto, Soft core hyperon - nucleon potentials, Phys. Rev. C 59, 21-40 (1999).
  * (50) V. G. J. Stoks and T. A. Rijken, Soft core baryon baryon potentials for the complete baryon octet, Phys. Rev. C 59, 3009-3020 (1999).
  * (51) Y. Oh and H. Kim, K* photoproduction off the nucleon: gamma N —$>$ K* Lambda, Phys. Rev. C 73, 065202 (2006).
  * (52) X. H. Liu and Q. Zhao, Further study of the helicity selection rule evading mechanism in $\eta_{c}$, $\chi_{c0}$ and $h_{c}$ decaying to baryon anti-baryon pairs, J. Phys. G 38, 035007 (2011).
  * (53) P. Colangelo, F. De Fazio and T. N. Pham, Nonfactorizable contributions in B decays to charmonium: The Case of B- —$>$ K- h(c), Phys. Rev. D 69, 054023 (2004).
  * (54) F. K. Guo, C. Hanhart, Q. Wang and Q. Zhao, Could the near-threshold $XYZ$ states be simply kinematic effects?, Phys. Rev. D 91, no.5, 051504 (2015).


# MS-HLMO: Multi-scale Histogram of Local Main Orientation for Remote Sensing
Image Registration

Chenzhong Gao, Wei Li, _Senior Member, IEEE_ , Ran Tao, _Senior Member, IEEE_
, Qian Du, _Fellow IEEE_ This work is supported by the National Key R&D
Program of China (Grant No. 2021YFB3900502) and National Natural Science
Foundation of China (61922013). (Corresponding author: Wei Li; e-mail:
liwei089@ieee.org)C. Gao, W. Li, R. Tao are with the School of Information and
Electronics, Beijing Institute Technology, Beijing 100081 China; Beijing Key
Laboratory of Fractional Signals and Systems, Beijing Institute of Technology,
Beijing 100081, China (e-mail: gao-pingqi@qq.com).Q. Du is with the Department
of Electrical and Computer Engineering, Mississippi State University,
Mississippi State, MS 39762, USA (e-mail: du@ece.msstate.edu).

###### Abstract

Multi-source image registration is challenging due to intensity, rotation, and
scale differences among the images. Considering the characteristics and
differences of multi-source remote sensing images, a feature-based
registration algorithm named Multi-scale Histogram of Local Main Orientation
(MS-HLMO) is proposed. Harris corner detection is first adopted to generate
feature points. The HLMO feature of each Harris feature point is extracted on
a Partial Main Orientation Map (PMOM) with a Generalized Gradient Location and
Orientation Histogram-like (GGLOH) feature descriptor, which provides high
intensity, rotation, and scale invariance. The feature points are matched
through a multi-scale matching strategy. Comprehensive experiments on 17
multi-source remote sensing scenes demonstrate that the proposed MS-HLMO and
its simplified version MS-HLMO+ outperform other competitive registration
algorithms in terms of effectiveness and generalization.

###### Index Terms:

Image registration, multi-source, remote sensing, multi-modal, multi-scale,
histogram of local main orientation (HLMO)

##  I Introduction

Image registration is a key preprocessing step in remote sensing applications.
With earth observation systems being developed rapidly in recent years, to
achieve in-depth ground analysis, the use of multi-source remote sensing
images becomes popular in achieving in-depth ground analysis [1, 2, 3, 4, 5].
The accurate alignment of the multi-source images is a prerequisite of
subsequent applications such as data fusion, change detection, joint analysis,
and other techniques.

The purpose of image registration is to spatially align multiple images for
the same scene. In general, remote sensing images are spatially corrected by
reference coordinates (such as fictitious graticule) or image control points.
However, there are inconsistent or missing spatial references. In particular,
in multi-source remote sensing, the inconsistency of spatial correction
between different sources is more obvious. This phenomenon can be seen from
the examples shown in Fig.1, which is about the correction of the multi-sensor
image solely through the spatial reference information provided by the data
source. It is obvious that the ground covers correspond inconsistently. In
addition, for some UAV or ground-based images, there is no spatial reference
information. It is very important to correct multi-source remote sensing
images automatically by registration algorithms in order to provide accurately
registered data for the same scene space.

Figure 1: Examples of alignment failure of multi-source remote sensing images
with only calibration information.

Existing automatic image registration algorithms are generally divided into
two categories, i.e., the traditional methods and the deep learning-based
methods. Among them, the deep learning-based methods [6] are relatively
advanced, but a large amount of annotated samples are needed for model
training, and the model is often highly targeted. The traditional ones are
systematically classified into area-based and feature-based methods [7]. The
area-based, also called intensity-based methods register images by
establishing the similarity measure of the intensity values or in a transform
domain [8, 9, 10, 11, 12, 13, 14, 15], and optimal geometric transformation
parameters between the image pairs are found through optimization [16]. Apart
from efficiency and other factors, for multi-source images, this type
algorithm is prone to fail when multi-modal differences are present.

Registration based on feature matching is a relatively mature technique, in
which the transformation parameters are estimated through the coordinate
correspondence of matched features. Scale-Invariant Feature Transform (SIFT)
[17, 18, 19] is one of the most classic, effective, and commonly used feature
extraction and matching methods, and became the basis of many improved
algorithms [20, 21, 22, 23, 24, 25, 26]. SIFT uses Gaussian scale-space,
supported by scale-space theory [27], for feature points extraction, which
also provides scale invariance. The feature points are then described by
statistical local gradient information for points matching. This also lays the
framework for the subsequent algorithms. In addition, many feature extraction
algorithms have been proposed and applied to image registration [28, 29, 30,
31, 32, 33, 34, 35, 36, 37].

However, most algorithms are only applicable to single-modal image
registration. Algorithms for multi-modal image registration are more critical.
Focusing on this, multi-modal image registration algorithms have been
proposed. Chen et al. [38] proposed Partial Intensity Invariant Feature
Descriptor (PIIFD) based on SIFT algorithm for multi-source retinal image
registration, which overcame the problem of intensity difference problems.
PSO-SIFT [39] presented a new gradient definition and an enhanced feature
matching method for the registration of remote sensing images with intensity
differences. Ye et al. developed Histogram of Orientated Phase Congruency
(HOPC) [40] and Local HOPC (LHOPC) [41] for multi-modal remote sensing
registration, in which Minimum Moment of Phase Congruency (MMPC)-Lap is used
for keypoints detection and an extended phase congruency model is used for
feature description. Radiation-variation insensitive feature transform (RIFT)
[42] utilized a maximum index map (MIM) as the feature map, which is invariant
to intensity diffeerence. The MIM is obtained by assigning the strongest
response of the log-Gabor filtering at several predetermined orientations as
the maximum index at each pixel. Multi-scale PIIFD (MS-PIIFD) [43] improved
PIIFD by adopting Gaussian scale-space, in which feature extraction and
matching are completed in a multi-scale strategy. It achieves scale robustness
and performs well on multi-modal images with scale differences.

The multi-modal differences of multi-source images result in multiple
inconsistency of local features in the same area of image pairs, which makes
it difficult to match correctly unless all inconsistencies are conquered. It
is needed to handle the differences and enhance the similarity of the
features. Therefore, through deep-going analysis of practical work and real
data, we refine and model the differences of multi-source remote sensing
images as

${{\bf{I}}_{2}}={{\cal F}_{{\rm{Cut}}}}({{\cal F}_{{\rm{Tran}}}}({{\cal
F}_{{\rm{Sc}}}}({{\cal F}_{{\rm{Rot}}}}({{\cal F}_{{\rm{Int}}}}({{\cal
F}_{{\rm{Chg}}}}({{\bf{I}}_{1}}))))))$ (1)

where ${{\bf{I}}_{1}}$ and ${{\bf{I}}_{2}}$ are a multi-source image pair,
${{\cal F}_{{\rm{Chg}}}}(\bullet)$ represents the changes of the real ground
object content sensed in the images, which is often caused by multi-temporal
sensing. ${{\cal F}_{\rm{Int}}}(\bullet)$ denotes intensity transformation.
${{\cal F}_{{\rm{Rot}}}}(\bullet)$, ${{\cal F}_{{\rm{Sc}}}}(\bullet)$, and
${{\cal F}_{{\rm{Tran}}}}(\bullet)$ are spatial transformation of rotation,
scale change, and translation, respectively, which are the basic operation of
similarity transformation. ${{\cal F}_{{\rm{Cut}}}}(\bullet)$ represents the
spatial cutting operation, which results in images covering different ranges
of vision. There may also be large perspective differences or severe local
geometric distortions between images, which leads to more complex image
distortion and registration modeling. Our goal is to first resolve the most
common and fundamental problem of image differences. Among the above multi-
modal differences, ${{\cal F}_{{\rm{Cut}}}}(\bullet)$, ${{\cal
F}_{{\rm{Tran}}}}$, and ${{\cal F}_{{\rm{Chg}}}}(\bullet)$ have been solved by
feature-based process. Thus, the primary goal is to find a feature robust to
${{\cal F}_{{\rm{Sc}}}}(\bullet)$, ${{\cal F}_{{\rm{Rot}}}}(\bullet)$, and
${{\cal F}_{\rm{Int}}}(\bullet)$ in multi-source images.

Based on the above analysis, a novel feature-based registration framework for
multi-source remote sensing images is proposed through designing effective
strategies for robust feature extraction that overcomes multiple image
differences. It is able to automatically register remote sensing images under
various conditions without human intervention. The main contributions are
summarized as follows:

  1. 1.

An image registration algorithm named Multi-scale Histogram of Local Main
Orientation (MS-HLMO) is designed to cope with various multi-source remote
sensing images with multi-modal differences. Through comprehensive
experimental verification, MS-HLMO is able to effectively deal with common
multi-source remote sensing image registration, including multi-sensor, multi-
temporal, and multi-viewpoint image pairs with different resolutions and data
size, which reflects strong robustness and generalization.

  2. 2.

A local feature named Histogram of Local Main Orientation (HLMO) is proposed,
in which a basic feature map called Partial Main Orientation Map (PMOM) is
developed for local feature extraction, which handles multi-modal differences
and provides robust orientation information. In addition, HLMO utilizes a
generalized GLOH-like (GGLOH) feature descriptor to extract local features.
The proposed HLMO is invariant to intensity and rotation, and has outstanding
ability in robust feature extraction from multi-source images.

  3. 3.

To overcome the problem of scale differences in multi-source images, a multi-
scale feature extraction and matching strategy based on Gaussian scale-space
is applied in MS-HLMO, which has a good effect on solving scale differences
and optimizing feature matching. MS-HLMO has a complete processing flow
without manual intervention and is able to directly applied to the practical
multi-source remote sensing joint analysis. To further research and assessment
the proposed method, the code will be released on https://github.com/MrPingQi.

The rest of this paper is arranged as follows. Section II describes the
overall framework of the proposed MS-HLMO and each critical process in detail.
Section III presents experiments and discussion in terms of discrete and
overall testing. Finally, this research is summarized in Section IV.

##  II Proposed Registration Method

Figure 2: The proposed multi-source remote sensing image registration
framework of MS-HLMO, including Harris feature point detection, Histogram of
Local Main Orientation feature extraction, and multi-scale registration
strategy.

The framework of the proposed MS-HLMO registration algorithm is shown in
Fig.2. The input multi-source image pair to be registered is preprocessed,
which includes data normalization and basic denoising. Then, the preprocessed
single-band images are used for feature points detection and feature
extraction. Harris corner point detection, which contains detail treatments
for multi-source images, is adopted to generate feature points between the
image pair for matching. The key process of the proposed HLMO feature
extraction is carried out in a multi-scale strategy, in which Gaussian
pyramids are built to create a scale-space of the images. The HLMO feature
descriptors of each Harris corner point are extracted on the PMOM of the
images. The feature points between the image pair are then matched according
to the descriptors, and Fast Sample Consensus (FSC) is carried out to remove
the outliers. The matching results in the scale-space are combined through a
multi-scale matching strategy. Finally, the spatial transformation between the
original image pair is determined by the coordinate relationship between
matched feature points according to a selected transformation model.

###  II-A Harris Feature Point Detection

Harris corner [29] is one of the most stable feature points, which is slightly
affected by intensity and scale difference and has high computational
efficiency [44, 43]. It has the advantage in multi-source remote sensing
images with multi-modal properties and large data size. Here, the similar
strategy [43] is used for feature points detection. The Harris corner response
of each pixel is calculated by:

$cornerness=\frac{{{\rm{det}}(\textbf{\emph{M}})}}{{{\rm{tr}}(\textbf{\emph{M}})}}$
(2)
$\textbf{\emph{M}}=\left[{\begin{array}[]{*{20}{l}}{\sum\limits_{{\bf{W}_{\sigma}}}{{{\bf{G}}_{x}}^{2}}}&{\sum\limits_{{\bf{W}_{\sigma}}}{{{\bf{G}}_{x}}{{\bf{G}}_{y}}}}\\\
{\sum\limits_{{\bf{W}_{\sigma}}}{{{\bf{G}}_{x}}{{\bf{G}}_{y}}}}&{\sum\limits_{{\bf{W}_{\sigma}}}{{{\bf{G}}_{y}}^{2}}}\end{array}}\right]$
(3)

where $\rm{det}(\textbf{\emph{M}})$ and $\rm{tr}(\textbf{\emph{M}})$ are the
determinant and trace of M, respectively, ${{\bf{G}}_{x}}$ and
${{\bf{G}}_{y}}$ are the image’s gradient along $x$ and $y$ directions,
respectively, and $\bf{W}_{\sigma}$ is a Gaussian window with variance
$\sigma$. Pixels with strong response are considered to be feature points with
distinct structure and stability between multi-source images.

An important issue in practical multi-source remote sensing image registration
is that the data size and scale relations between images to be registered are
diverse. For example, an image with high resolution covers a smaller spatial
area. Many of the existing algorithms only deal with the ideal case that the
image pair has the same scale and size. This paper focuses on solving several
key problems at the same time, that is, two uncertain factors of image scale
and size should be considered simultaneously. The proposed MS-HLMO adopts
local non-maximum suppression (LNMS) to solve this problem. Since the size and
scale difference of the image pairs are uncertain, in Harris corner detection,
it is expected that the feature points in the image pair are distributed as
uniformly as possible with the use of LNSM. Then, the ratio of the window size
in LNMS is set depending on the ratio of the data size of the image pair:

$ratio=\sqrt{\frac{{M\times N}}{{m\times{\rm{n}}}}}$ (4)

where $M,N$ and $m,n$ are the length and width of the two images,
respectively.

(a)

(b)

Figure 3: Examples of feature points detection results with LNMS. (a)
Detection results of the image pair with different scale. (b) Detection
results of the image pair with different scale and size.

Consequently, the distribution of feature points is consistent with the
images’ scale proportion when there are scale differences, as shown in Fig.3a.
When there is a size difference, the feature points are far more uniformly
distributed, and the repeatability is higher, as shown in Fig.3b.

###  II-B Histogram of Local Main Orientation

####  II-B1 Partial main orientation map

Feature-point-based registration algorithms often use a local descriptor to
extract the neighborhood information of the keypoints, and generate their
feature vectors for similarity matching. For example, SIFT [19], HOG [30],
SURF [22] and PIIFD [38] employ gradient information as the basic feature.
However, the performance of these algorithms is greatly degraded when
processing multi-source, especially multi-sensor images. So, it is critical to
extract invariant feature that is robust to ${{\cal F}_{\rm{Int}}}(\bullet)$,
${{\cal F}_{\rm{Rot}}}(\bullet)$, and ${{\cal F}_{\rm{Chg}}}(\bullet)$ for
feature points description. A partial main orientation map (PMOM) is designed
as the feature map in MS-HLMO, in which the Average Squared Gradient (ASG)[45]
is adopted.

The ASG is a gradient weighting method. The elementary gradient of the image
along $x$ and $y$ directions, i.e., ${\bf{G}}_{x}$ and ${\bf{G}}_{y}$ are
calculated as

$\left[\begin{array}[]{l}{{\bf{G}}_{x}}(x,y)\\\
{{\bf{G}}_{y}}(x,y)\end{array}\right]=\left[\begin{array}[]{l}\frac{\partial}{{\partial
x}}{\bf{I}}(x,y)\\\ \frac{\partial}{{\partial
y}}{\bf{I}}(x,y)\end{array}\right]$ (5)

where ${\bf{I}}(x,y)$ represents the single-layer gray-scale image. The
magnitude and orientation of its gradient, i.e., ${\bf{G}}_{\rho}$ and
${\bf{G}}_{\varphi}$ are

$\left[\begin{array}[]{l}{{\bf{G}}_{\rho}}\\\
{{\bf{G}}_{\varphi}}\end{array}\right]=\left[\begin{array}[]{l}\sqrt{{{\bf{G}}_{x}}^{2}+{{\bf{G}}_{y}}^{2}}\\\
\arctan\frac{{\bf{G}}_{y}}{{\bf{G}}_{x}}\end{array}\right]$ (6)

In the ASG, a locally weighted squared gradient [45] along $x$ and $y$
directions, ${\bf{G}}_{{{\bf{W}}_{\sigma}},s,x}$ and
${\bf{G}}_{{\bf{W}_{\sigma}},s,y}$ are obtained as

$\left[\begin{array}[]{l}{{\bf{G}}_{{{\bf{W}}_{\sigma}},s,x}}\\\
{{\bf{G}}_{{{\bf{W}}_{\sigma}},s,y}}\end{array}\right]=\left[\begin{array}[]{l}\sum\limits_{{\bf{W}}_{\sigma}}{{{\bf{G}}_{x}}^{2}-{{\bf{G}}_{y}}^{2}}\\\
\sum\limits_{{\bf{W}}_{\sigma}}{2{{\bf{G}}_{x}}{{\bf{G}}_{y}}}\end{array}\right]$
(7)

where ${\bf{W}}_{\sigma}$ is a Gaussian window with variance $\sigma$.
Accordingly, the orientation of this gradient is

${{\bf{G}}_{{{\bf{W}}_{\sigma}},s,\varphi}}=\angle({{\bf{G}}_{{{\bf{W}}_{\sigma}},s,x}},{{\bf{G}}_{{{\bf{W}}_{\sigma}},s,y}})$
(8)

where $\angle(X,Y)$ is defined as

$\angle(X,Y)=\left\\{\begin{array}[]{l}\arctan(\frac{Y}{X}),X\geq 0\\\
\arctan(\frac{Y}{X})+\pi,X<0,Y\geq 0\\\
\arctan(\frac{Y}{X})-\pi,X<0,Y<0\end{array}\right.$ (9)

making ${\bf{G}}_{{{\bf{W}}_{\sigma}},s,\varphi}$ within $(-\pi,\pi)$.
According to [45], this gradient is obtained by doubling the angle of the
original gradient, so the orientation of the ASG is

${{\bf{G}}_{{{\bf{W}}_{\sigma}},\varphi}}=\frac{1}{2}{{\bf{G}}_{{{\bf{W}}_{\sigma}},s,\varphi}}$
(10)

Compared with the classical gradient operator, the ASG orientation
${{\bf{G}}_{{{\bf{W}}_{\sigma}},\varphi}}\in(-\frac{\pi}{2},\frac{\pi}{2})$
reflects the weighted gradient orientation of a local region
${\bf{W}}_{\sigma}$, which is more robust and stable. In addition, the $x$
direction gradient is constant, and this characteristic meets the requirement
that not affected by the reversal of gradient in intensity difference. Note
that when $\sigma$ increases, the scale of ASG increases, which makes the
local orientation more invariant to intensity difference and noise, but the
uniqueness of local features decreases. From this, the following function is
defined:

${{\bf{G}}_{PMOM}}=\frac{1}{2}\angle(\sum\limits_{\sigma}{\sum\limits_{{\bf{W}}_{\sigma}}{{{\bf{G}}_{x}}^{2}-{{\bf{G}}_{y}}^{2}}},\sum\limits_{\sigma}{\sum\limits_{{\bf{W}}_{\sigma}}{2{{\bf{G}}_{x}}{{\bf{G}}_{y}}}})$
(11)

where a series of scale $\sigma$ are preset, the weighted responses
${\bf{G}}_{{{\bf{W}}_{\sigma}},s,x}$ and ${\bf{G}}_{{{\bf{W}}_{\sigma}},s,y}$
at each scale are added, and the ASG orientation is obtained. By filtering the
image with Eq.(11), the PMOM is obtained, where its value reflects the overall
orientation of multiple scales in each partial area of the image.

Figure 4: Comparison of feature maps of a selected scene, including a visible
and an infrared image.

A visualized comparison of PMOM with other feature maps of typical multi-
source data is shown in Fig.4. The original data is a visible-infrared image
pair, which contains obvious intensity difference. For comparison purposes,
the images have been registered manually, basically eliminating the spatial
differences of scale, rotation, and size. The magnitude and orientation of
images’ gradient are shown in Fig.4, which are obtained using Eqs.(5)-(6).
These are the basic feature information in most algorithms [19, 30, 22, 38,
43]. It is observed that, due to the multi-modal properties of the original
image pair, these two feature maps have large differences and instability,
which is the main reason for the failure of most algorithms. The MIM [42] of
RIFT shown in Fig.4 also focuses on the local orientation of the image, where
the maximum index is the main orientation among several ones. Compared with
the directional gradient, Gabor transformation has a more stable response,
which leads to RIFT robust to intensity difference. However, its value will
also mutate due to local changes in images, and the rotation invariance is
slightly poor. The image pair’s PMOMs are shown at the bottom of Fig.4.
Compared with MIM, the proposed PMOM is not only more robust and stable
between multi-modal images, but also continuous in value, which is conducive
to achieving effective rotation invariance. In HLMO, PMOM is used as the
unique feature information to extract local features of keypoints that are
invariant to multi-modal properties.

####  II-B2 Descriptor extraction

After determining the feature points and the specific feature for
discrimination, the next step is to make use of the local feature information
around each point and generate descriptors. Gradient Location and Orientation
Histogram (GLOH) has shown excellent ability through experiments [21], and has
been applied in multi-source remote sensing image registration [25, 39]. The
original GLOH descriptor is a circular region divided by three circles,
similar to that shown in Fig.6, in which the two outer circular regions are
divided into 4 parts, and the radius of the circular region divided are 5, 9,
and 11. The partition size and the number are then improved [21, 25, 39].
However, different parameters have various effects when treating multifarious
types of images. In addition, if the number of outer ring regions is too
small, the character of feature points is not significant, which makes it
difficult to match accurately. If it is too large, the features are unstable,
and the dimension of the descriptor is too high, which increases the burden of
redundant calculation. To deal with this, a generalized GLOH-like (GGLOH)
descriptor is proposed, and its structure of which is shown in Fig.5a.

(a)

(b)

Figure 5: Descriptor structure of the proposed GGLOH. (a) Subregion partition
of the local neighborhood of a feature point. (b) Angle quantification within
each subregion.

Let ${{\rm{A}}^{0}}$ denote the central circular region, and
${\rm{A}}_{j}^{i},(i=1,2,j=1,...,{N_{\rm{A}}})$ represent the sector subregion
$j$ in the outer ring region $i$. Let ${N_{\rm{A}}}$ be the number of the
subregions in each out ring region, which is even, ${\theta_{\rm{0}}}$ be the
main orientation of the feature point, and $R_{0}$, $R_{1}$, $R_{2}$ be the
radii of the central and outer regions, respectively. Note that the
orientations of pixels’ gradient in each region are counted as feature
information, therefore, fair use of information in each region is expected.
The number of pixels in each region should be roughly the same, and the weight
of the outer regions should not change due to the change of ${N_{\rm{A}}}$. So
the area of each region should be the same, that is

${N_{\rm{A}}}\cdot\pi{\rm{R}}_{\rm{0}}^{2}=\pi({\rm{R}}_{\rm{1}}^{2}-{\rm{R}}_{\rm{0}}^{2})=\pi({\rm{R}}_{\rm{2}}^{2}-{\rm{R}}_{\rm{1}}^{2})$
(12)

which also fixes the relationship between $R_{0}$, $R_{1}$, $R_{2}$ and
${N_{\rm{A}}}$. When ${N_{\rm{A}}}$ is given different values, the stability
and importance of each region’s feature remains the same. In HLMO, the GGLOH
is used to extract local features on the PMOM, where the orientation values
within $(-\frac{\pi}{2},\frac{\pi}{2})$ are uniformly quantified to
${N_{\rm{O}}}$ values, as shown in Fig.5b, where
${\phi_{k}}(k=1,2,...,{N_{\rm{O}}})$ is the angle after quantization. A
histogram with ${N_{\rm{O}}}$ bins is obtained in each region.

It is simple to achieve rotation invariance of HLMO. For each keypoint, the
PMOM value at its location is the main orientation, that is, the reference
orientation ${\theta_{0}}$ of the GGLOH. Then, all of the PMOM values within
the local area of GGLOH also take ${\theta_{0}}$ as the reference (0°), that
is, all angle values minus ${\theta_{0}}$, and those beyond
$(-\frac{\pi}{2},\frac{\pi}{2})$ are flipped to their opposite angles.

Figure 6: The problem caused by the jump of main orientation near
$-\frac{\pi}{2}$ or $\frac{\pi}{2}$.

Another key problem is that the rotation and nonlinear intensity difference
may cause the jump of the main orientations of some feature points near
$-\frac{\pi}{2}$ and $\frac{\pi}{2}$. An example is shown in Fig.6. In PIIFD
[38], a similar problem has been discovered and improvement has been made for
SIFT. Then a similar strategy is adopted to process GLOH-like descriptors,

${{\bf{D}}_{1}}=\left[{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{{\bf{H}}_{1}^{1}}&{{\bf{H}}_{2}^{1}}&\cdots&{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}}^{1}}\end{array}}\\\
{\begin{array}[]{*{20}{c}}{{\bf{H}}_{1}^{2}}&{{\bf{H}}_{2}^{2}}&\cdots&{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}}^{2}}\end{array}}\end{array}}\right]$
(13)

${{\bf{D}}_{2}}=\left[{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}+1}^{1}}&{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}+2}^{1}}&\cdots&{{\bf{H}}_{{N_{\rm{A}}}}^{1}}\end{array}}\\\
{\begin{array}[]{*{20}{c}}{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}+1}^{2}}&{{\bf{H}}_{{{{N_{\rm{A}}}}\mathord{\left/{\vphantom{{{N_{\rm{A}}}}2}}\right.\kern-1.2pt}2}+2}^{2}}&\cdots&{{\bf{H}}_{{N_{\rm{A}}}}^{2}}\end{array}}\end{array}}\right]$
(14)

${\bf{D}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{D}}_{1}}+{{\bf{D}}_{2}}}\\\
{c\left|{{{\bf{D}}_{1}}-{{\bf{D}}_{2}}}\right|}\end{array}}\right]$ (15)

where ${\bf{H}}_{j}^{i}$ is the histogram vector of gradient orientation of
region ${\rm{A}}_{j}^{i}$. In this way, no matter whether the main orientation
of feature points is reversed 180° or not, descriptor ${\bf{D}}$ is composed
of the addition and subtraction of the upper and lower parts of GGLOH
according to the main orientation axis, without changing the regions’ order.
Finally, a descriptor vector ${\bf{D}}_{P}$ is generated for the feature point
$P$, whose dimension is $(2\cdot{N_{\rm{A}}}+1)\cdot{N_{\rm{O}}}$.

###  II-C Multi-scale Registration Strategy

Scale difference ${{\cal F}_{{\rm{Sc}}}}(\bullet)$ of multi-source images has
a great influence on local features. Some algorithms have quantitative scale
judging methods, such as SIFT [19] and LHOPC. However, it is found that these
methods are invalid in images with large modal differences. The reason is that
when images do not belong to the same degradation model, it is not credible to
judge the scale quantitatively through local image feature information. Multi-
source images often have scale differences, and sometimes the scale proportion
is unknown. In order to deal with this key problem and realize scale
robustness, a multi-scale feature extraction and matching strategy is designed
in MS-HLMO.

Figure 7: Structure of the Gaussian pyramid used in MS-HLMO.

Local information of feature points is extracted in the scale-space of the
images. Based on the scale-space theory [27], the method of building image’s
Gaussian pyramids is adopted. The schematic diagram of establishing Gaussian
pyramid of the image in the proposed algorithm is shown in Fig.7. The original
image is first sampled down step by step to obtain a series of images with
different resolutions, that is, the first layer in each octave. Then in each
octave, a series of Gaussian blurs are performed:

${\bf{L}}={\bf{G}}*{\bf{I}}$ (16)
${\bf{G}}=\frac{1}{{\sqrt{2\pi{\sigma^{2}}}}}{e^{\frac{{-({x^{2}}+{y^{2}})}}{{2{\sigma^{2}}}}}}$
(17)

where ${\bf{I}}$ is the original image, ${\bf{G}}$ is a Gaussian kernel with a
standard deviation of $\sigma$, and ${\bf{L}}$ is the Gaussian blur image.

Algorithm 1 Proposed MS-HLMO Feature Extraction

Input: single-band image $\mathbf{I}$, feature point set $P_{\mathbf{I}}$,
total number of octaves ${N_{\rm{GO}}}$ and layers ${N_{\rm{GL}}}$ in Gaussian
pyramid, subregion and angle partition parameters ${N_{\rm{A}}}$,
${N_{\rm{O}}}$ in GGLOH , patch size $S$ of HLMO.

Through down-sampling and Eq.(16)(17), the Gaussian pyramid
${\bf{G}}_{\bf{I}}(O,L)$ of image $\mathbf{I}$ is built with ${N_{\rm{GO}}}$
octaves and ${N_{\rm{GL}}}$ layers in each octave.

In each layer of ${\bf{G}}_{\bf{I}}(O,L)$:

Calculate the PMOM of this layer to get ${\bf{F}}_{\bf{I}}(O,L)$ according to
Eq.(5)(11)

For each feature point $P$ in $P_{\bf{I}}$:

Calculate the corresponding position

Take the PMOM value at the position as the main orientation ${\theta_{0}}$

Taking the main orientation as the reference direction ($0^{\circ}$),
establish a GGLOH window with size (diameter) of S

Statistics the PMOM value within each region of GGLOH to obtain the basic
feature descriptor $D_{1}(P,O,L)$ and $D_{2}(P,O,L)$

Obtain the descriptor $D(P,O,L)$ of $P$ with Eq.15.

Output: feature descriptor set $D_{P_{\bf{I}}}(O,L)$

After the scale-space of the images is established, for each Harris corner
point, the HLMO descriptor is calculated by obtaining the local information at
the corresponding location of each feature point in the scale-space. The
proposed multi-scale HLMO feature extraction method is provided in Algorithm
1, where $O$ is the octave number in the Gaussian pyramid, and $L$ is the
layer number. The algorithm outputs the feature point descriptor set
$D_{P_{\bf{I}}}(O,L)$, which contains
$(2\cdot{N_{\rm{A}}}+1)\cdot{N_{\rm{O}}}$-dimensional vectors for each feature
point at each scale.

The next step is to match the feature point sets of the image pair according
to the descriptor sets. The process of the multi-scale feature matching is
provided in Algorithm 2. In the process, each scale is matched in turn. Then
the matching results are merged step by step while the outlier removal is
carried out to realize the optimization of matching points. The final matching
results are used to determine the spatial transformation between images. The
most critical is to combine all the matching results of feature points and
remove outliers, so as to maximize the correct matches of all scales. Fig.8
shows this process visually. Obviously, this is a general approach to handle
all kinds of unknown scale differences. When the scale proportion of images is
known or can be estimated, then the above process can be greatly simplified.
In this case, the proposed multi-scale strategy still has the advantages of
enhancing feature matching and maximizing the number of matching points.

Algorithm 2 Proposed MS-HLMO Feature Matching

Input: feature point set of the image pair $P_{{\bf{I}}1}$, $P_{{\bf{I}}2}$,
feature descriptor set of the image pair $D_{P_{{\bf{I}}1}}(O_{1},L_{1})$,
$D_{P_{{\bf{I}}2}}(O_{2},L_{2})$

Take each layer of $D_{P_{{\bf{I}}1}}(O_{1},L_{1})$:

Take each layer of $D_{P_{{\bf{I}}2}}(O_{2},L_{2})$:

Match $P_{{\bf{I}}1}$ and $P_{{\bf{I}}2}$ using Euclidean distance of the
descriptorss

Remove outliers, producing the matching result of a single scale
$M(O_{1},O_{2},L_{1},L_{2})$

The matching results of all layers in each octave of the scale-space are union
and then optimized with outlier removal, producing the matching result
$M_{L}(O_{1},O_{2})$

The matching results of all octaves in $M_{L}(O_{1},O_{2})$ are union and then
optimized with outlier removal, producing the final matching result
$M_{OL}(P_{{\bf{I}}1},P_{{\bf{I}}2})$

Output: feature points matching set $M_{OL}(P_{{\bf{I}}1},P_{{\bf{I}}2})$

Figure 8: Multi-scale keypoints matching strategy in MS-HLMO.

##  III Experimental Results and Discussion

To validate the registration performance of the proposed framework, a set of
typical multi-source remote sensing images are carefully selected, and six
effective registration algorithms are used for comparison, including SIFT
[19], SAR-SIFT [39], PSO-SIFT [39], PIIFD [38], MS-PIIFD [43], and RIFT [42],
all of which contain processes dealing with multi-modal properties. The
experiments are implemented using MATLAB2021a on a platform with AMD-
Ryzen-7-5800X 3.80GHz CPU, 32GB RAM.

###  III-A Experimental Data and Preprocessing

The data used in experiments are multi-source remote sensing images of 17
scenes, labeled a$\sim$q. These images include hyperspectral image (HSI),
multispectral (MSI), synthetic aperture radar (SAR), visible images, infrared
images, depth maps such as LiDAR-derived DSM (digital surface model), and even
artificial maps. They contain day-night image pair and images from different
years. The scenes cover spaceborne, airborne, and ground-based remote sensing
data. There are various intensity and geometric distortions with the images.
Table I shows the types and sizes of each set of data, as well as whether they
have significant intensity difference, rotation, or scale difference. The
specific form of the scenes is shown in Fig.13 and 14. Except for artificial
maps, these are all real data without manual manipulation.

TABLE I: Multi-source remote sensing data of 17 scenes. Scene | Type | Size | Intensity | Rotation | Scale
---|---|---|---|---|---
Difference | Difference
a | a1 | MSI | 3555$\times$4026 | $\checkmark$ |  | $\checkmark$
HSI | 1185$\times$1342 |
SAR | 3555$\times$4026 |
a2 | MSI | 1287$\times$2035 | $\checkmark$ |  | $\checkmark$
HSI | 422$\times$678 |
SAR | 1366$\times$887 |
b | MSI | 3529$\times$1756 | $\checkmark$ |  | $\checkmark$
HSI | 1175$\times$585 |
c | MSI | 2461$\times$4139 | $\checkmark$ |  | $\checkmark$
MSI | 2271$\times$3152 |
d | PAN | 1000$\times$1000 | $\checkmark$ |  | $\checkmark$
HSI | 550$\times$450 |
e | visible | 614$\times$614 | $\checkmark$ | $\checkmark$ |
infrared | 611$\times$611 |
f | visible | 617$\times$593 | $\checkmark$ | $\checkmark$ |
infrared | 617$\times$593 |
g | visible | 1920$\times$1080 | $\checkmark$ |  | $\checkmark$
infrared | 667$\times$504 |
h | visible | 4056$\times$3040 | $\checkmark$ |  | $\checkmark$
infrared | 640$\times$512 |
i | visible | 633$\times$460 | $\checkmark$ |  |
infrared | 670$\times$508 |  |
j | HSI | 349$\times$1905 | $\checkmark$ |  |
LiDAR | 349$\times$1905 |  |
k | HSI | 166$\times$600 | $\checkmark$ |  |
LiDAR | 166$\times$600 |  |
l | HSI | 325$\times$220 | $\checkmark$ |  |
LiDAR | 325$\times$220 |  |
m | visible | 500$\times$500 | $\checkmark$ |  |
depth | 500$\times$500 |  |
n | visible | 500$\times$500 | $\checkmark$ |  |
map | 500$\times$500 |  |
o | visible | 500$\times$500 | $\checkmark$ |  |
visible | 500$\times$500 |  |
p | SAR | 600$\times$500 | $\checkmark$ | $\checkmark$ |
SAR | 600$\times$500 |
q | visible | 300$\times$300 |  | $\checkmark$ |
visible | 300$\times$300 |  |

As MS-HLMO and other methods only focus on spatial rather than waveband
features, while some types of images have multiple wavebands, such as MSI,
HSI, or LiDAR, single-band gray-scale images need to be obtained from the
original multi-band ones. It is acceptable to use data dimension reduction,
such as Principal Component Analysis (PCA), to obtain single-band images.
However, in multi-source registration, the images may have multi-modal
properties. Due to its enhancing pixel-level features, data dimension
reduction could break the stable and unified spatial structures of each type
of ground covers in the image, and increase the structural differences between
images. Then, registration would be more difficult as ${{\cal
F}_{{\rm{Chg}}}}(\bullet)$ is exacerbated, which is undesired. Note that HLMO
is not subject to the changes in image’s grayscale, so it is optional to
select any one of the bands as the input, or just simply add up all the bands,
which has little effect on the registration result.

###  III-B Experiment Settings

The proposed MS-HLMO contains seven main parameters that greatly affect the
registration, which are the number of feature points, the scale of PMOM, the
$R_{2}$, ${N_{\rm{A}}}$, and ${N_{\rm{O}}}$ of GGLOH, and ${N_{\rm{GO}}}$ and
${N_{\rm{GL}}}$ of the Gaussian pyramid. In general, the lower the number of
subregions ${N_{\rm{A}}}$ and angles ${N_{\rm{O}}}$ are, the higher the
robustness of features is, but the lower the separability of the feature
points is. Conversely, the more distinct the features are, but the lower the
stability is, which also increases the computational burden. In the
experiment, both ${N_{\rm{A}}}$ and ${N_{\rm{O}}}$ are set to 12, and $R_{2}$
is set to 48, which has the best performance and stability in the majority of
multi-source data. Then, $R_{1}$ and $R_{0}$ in GGLOH are determined by Eq.12,
and the patchsize of HLMO is $2\times R_{2}=96$. Through a simple set of
tests, the threshold of the Harris corner points’ number is set at 2000, and
${N_{\rm{GO}}}$ and ${N_{\rm{GL}}}$ are set to 3 and 4, which handles the
registration effectively enough with high efficiency.

In the calculation of PMOM, 10 scales are set for fusion, of which the
radiuses are evenly spaced between $R_{0}$ and $R_{2}$, and the $\sigma$ is
1/3 of each radius. In this way, when GGLOH is applied to a feature point, the
main orientation of the patch, that is the PMOM value of the center point, is
only determined by the pixels’ values within the patch, excluding information
beyond the GGLOH region. 10 scales are sufficient to fully obtain the local
multi-scale weighted main orientation information in various multi-source
images, which balances the robustness and uniqueness with a low computational
burden.

The number of correct matches (NCM) is employed as the main evaluation metric
of feature matching. The correct matches refer to the matched points that are
basically at the same position in the real space, and the NCM is positively
correlated with the accuracy of the transformation model. Note that if the NCM
of an image pair is less than 3, the parameters of the spatial transformation
model cannot be solved, then the registration process is considered to fail.

To ensure the fairness of the comparison, except for the feature point
detection and descriptor extraction methods, all the processes are the same,
including preprocessing, feature point matching, outlier removal,
transformation, etc. The feature points used in PIIFD, MS-PIIFD, and RIFT are
the same as the proposed algorithm with the strategy in [43]. Different
algorithms use various feature matching and outlier removal methods in their
original process, which interferes with the comparison of registration
effects. In this experiment, Euclidean distance measurement is used for all
feature matching, and FSC is used for outlier removal for a fair comparison.
As FSC has random property, the result fluctuates in the experiment, which is
not very large but does have an impact. So each test is repeated 50 times and
the result with the most NCM is taken.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 9: NCMs of MS-HLMO on different types of multi-source remote sensing
scenes as the rotation angles from 0° to 359°. (a) MSI-HSI. (b) MSI-SAR. (c)
MSI-map. (d) PAN-HSI. (e) visible-infrared. (f) HSI-LiDAR. (g) visible-depth.
(h) SAR-SAR.

###  III-C Invariance Tests

The crux of multi-source image registration based on feature matching lies in
various multimodal differences, and the robustness of local features to these
differences is crucial. According to the analysis, intensity difference
${{\cal F}_{\rm{Int}}}(\bullet)$, rotation ${{\cal F}_{{\rm{Rot}}}}(\bullet)$,
and scale difference ${{\cal F}_{{\rm{Sc}}}}(\bullet)$ are the most common
problems. So we design experiments to singly test the invariance of MS-HLMO
from these three aspects, so as to test whether it deals with each multimodal
problem between multi-source images.

####  III-C1 Intensity invariance

In order to independently test the intensity invariance of MS-HLMO, the
influence of other factors is excluded. Several image pairs are selected and
manually registered in advance, providing image pairs of the same size that
basically eliminated the differences in scale and rotation as the input image
of registration. The NCMs of the registration results are listed in Table II.
The MS-HLMO+ represents the proposed algorithm without the operation of
rotation invariance, that is, the main orientation is not assigned to each
feature point, and the reference orientation ${\theta_{\rm{0}}}$ in feature
description is set to 0. In this way, the robustness to intensity difference
is evaluated more accurately, and the influence of descriptor rotation and
PMOM value changes is excluded. As RIFT is an effective feature focusing on
intensity difference, it is chosen for comparison.

TABLE II: NCMs of 8 scenes with only intensity difference by MS-HLMO and RIFT. Method | Scene
---|---
a1 | a2 | d | j | g | m | a4 | p
RIFT+ | 1032 | 385 | 393 | 176 | 104 | 189 | 33 | 92
MS-HLMO | 1417 | 620 | 1321 | 149 | 280 | 374 | 220 | 575
MS-HLMO+ | 1434 | 1017 | 1454 | 771 | 351 | 846 | 386 | 868

The 8 scenes basically cover common remote sensing image types, including
various optical images, SAR, depth map, and artificial map. From the numerical
results, each group has obtained sufficient correct matches, especially in
scene (a1) and (d), about 1400 of 2000 detected Harris corner points have been
successfully matched. In most scenes, the NCMs of the MS-HLMO are far more
than RIFT, reflecting stronger robustness to intensity difference. In scene
(j), the MS-HLMO with rotation invariance obtained fewer matches than RIFT.
Through analysis, this is because scene (j) is an HSI-LiDAR image pair with
large size and complex content, where the modal difference between them is
large. Thus, the stability of the feature points’ main orientation is not
high. The inconsistent reference direction leads to large differences in some
feature points’ descriptors, while RIFT does not have influence caused by the
feature point’s orientation. However, without rotation interference, the MS-
HLMO+ obtains a significant number of correct matches, more than four times to
RIFT. The experiment shows that MS-HLMO with PMOM as its basic feature has
strong intensity invariance and is very effective for multi-modal image
registration.

####  III-C2 Rotation invariance

To test rotation invariance, for each one of the above image pairs without
scale and rotation differences, one of them is fixed and the other is rotated.
The rotation angles are from 0° to 359° with an interval of 1°, from which a
total of 360 image pairs of each scene are obtained. Then, the registration
algorithms are performed on each image pair, and the corresponding NCMs are
plotted in Fig.9. Samples of registration results under different rotation
angles are shown in Fig.10.

Figure 10: Examples of feature matching results with image rotation.

It is observed that the NCMs of the 8 scenes vary as the rotation angles
change, but are all greater than 100. None of the image pairs fails to be
registered at any angles, indicating that the proposed MS-HLMO has rotation
invariance over the entire angle interval of 0° to 360°. In most scenes, as
the images rotate, the NCMs fluctuate with a cycle of 90°, the reason of this
phenomenon is that the image rotation involves data interpolation as a digital
image is stored a matrix. In other words, the change of the local gradient
orientation is not linear, which is inevitable in digital image processing. As
a result, the main orientation of each feature point is not completely stable
with a certain offset from the image rotation angle. Relatively, the rotation
of 0°, 90°, 180°, and 270° do not change the data value, and the main
orientation of each feature point will not shift, resulting in the four peaks
in the curve. Theoretically, the four intervals evenly divided by 90° from 0°
to 360° are completely consistent, that is, the descriptors of a feature point
are completely the same when rotated by 90°, 180°, and 270° which is confirmed
by the experimental results. MS-HLMO has a strong rotation invariance and
fully copes with image rotation from various angles. In addition, the rotation
invariance is not restricted by the types of imaging sensors.

####  III-C3 Scale invariance

To evaluate the scale invariance of MS-HLMO, the same strategy is used by
fixing one image and scaling the other. The scaling ratios are from 1 to 2
with an interval of 0.1, providing a total of 11 image pairs in each scene. An
example of the feature points matching results of MS-HLMO and MS-HLMO+ are
shown in Fig.11-12, respectively.

Figure 11: NCMs of MS-HLMO and MS-HLMO+ on scene (m) as the scaling ratios
from 1 to 2.

(a)

(b)

Figure 12: Examples of feature matching results of scene (m) with scale
differences. (a) Matching results by MS-HLMO with scale ratio of 1.5 and 2.
(b) Matching results by MS-HLMO+ with scale ratio of 1.5 and 2.

Even if there are scale differences between images with integer or non-integer
proportions, the algorithm remains effective and obtains accurate matching. By
observing the NCM curve, with the gradual change of the scale ratio from 1 to
2, the number of matching points gradually decreases, and reaches the lowest
when the scale ratio was 1.5, and then gradually increases. This is because,
in the experiment, the down-sampling ratio in building gaussian scale-space is
set as 2. The interference caused by scale difference is well resisted when
the sampling scale is integer multiple of 2. In other scale proportions, a
significant difference of the receptive field in ground covers of the feature
descriptor is caused between the image pair. So the local windows used for
feature extraction between the two images are different in the receptive
field. In practice, the sampling step in Gaussian pyramid is further refined
to achieve better scale invariance.

This experiment demonstrates that MS-HLMO features are also robust to a
certain extent when images are at scale proportion between two octaves in
scale space and The proposed MS-HLMO has good robustness to scale difference
of multi-source images, which plays an extremely positive role in the actual
multi-source remote sensing image registration.

TABLE III: NCMs of the 17 remote sensing scenes comparison by eight registration methods. Method | Scene
---|---
a1 | a2 | a3 | a4 | b | c | d | e | f | g | h | i | j | k | l | m | n | o | p | q
SIFT | 794 | 30 | 25 | —— | —— | 145 | 86 | —— | —— | —— | —— | —— | —— | —— | —— | —— | —— | 15 | 14 | 256
SAR-SIFT | 16 | 6 | —— | —— | —— | 66 | 4 | —— | —— | —— | —— | —— | —— | —— | —— | —— | —— | 9 | 36 | 262
PSO-SIFT | 3 | —— | —— | —— | —— | —— | —— | 171 | 13 | —— | —— | —— | 9 | 6 | 10 | —— | —— | 19 | 24 | 259
PIIFD | —— | —— | —— | —— | —— | —— | —— | 53 | —— | —— | —— | —— | —— | —— | —— | —— | —— | —— | —— | 493
MS-PIIFD | 1105 | —— | 481 | —— | 444 | 145 | 999 | 460 | 78 | —— | —— | 25 | 5 | 339 | 30 | 74 | 165 | —— | 322 | 860
RIFT+ | —— | 21 | —— | 34 | —— | —— | —— | —— | —— | —— | —— | 39 | 176 | 109 | 188 | 189 | 96 | 179 | —— | ——
MS-HLMO | 1126 | 112 | 865 | 119 | 603 | 394 | 1030 | 689 | 161 | 130 | 231 | 79 | 149 | 511 | 420 | 374 | 161 | 97 | 432 | 883
MS-HLMO+ | 1194 | 321 | 996 | 149 | 984 | 564 | 1163 | —— | 693 | 162 | 337 | 145 | 761 | 789 | 861 | 846 | 273 | 451 | —— | ——

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

Figure 13: Feature points matching results of the 17 remote sensing scenes by
the proposed MS-HLMO. (a) a2:SAR-MSI, a1:MSI-HSI, a4:map-MSI, and a3:MSI-HSI
(multi-temporal). (b) MSI-HSI. (c) MSI-MSI. (d) PAN-HSI. (e) visible-infrared.
(f) visible-infrared. (g) visible-infrared. (h) visible-infrared. (i) visible-
infrared. (j) HSI-LiDAR. (k) HSI-LiDAR. (l) HSI-LiDAR. (m) visible-depth. (n)
visible-map. (o) day-night (visible). (p) multi-source SAR. (q) single-source
visible.

###  III-D Registration Performance

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

(m)

(n)

(o)

(p)

(q)

Figure 14: Registration results of the 17 remote sensing scenes by the
proposed MS-HLMO, shown in fusion form and alternation form. Feature points
matching results of the 17 remote sensing scenes by the proposed MS-HLMO. (a)
a2:SAR-MSI, a1:MSI-HSI, a3:MSI-HSI, and a4:map-MSI. (b) MSI-HSI. (c) MSI-MSI.
(d) PAN-HSI. (e) visible-infrared. (f) visible-infrared. (g) visible-infrared.
(h) visible-infrared. (i) visible-infrared. (j) HSI-LiDAR. (k) HSI-LiDAR. (l)
HSI-LiDAR. (m) visible-depth. (n) visible-map. (o) day-night (visible). (p)
multi-source SAR. (q) single-source visible.

In this section, the registration is performed directly on each original data
of the 17 scenes, to evaluate the overall registration performance. The NCMs
of the eight registration methods with all 17 scenes are listed in Table III.
Fig.13 shows the corresponding visualization of feature point matching results
of each scene. Due to space limitation, the one with the better result between
MS-HLMO and MS-HLMO+ is selected to display in each group.

In most scenes, MS-HLMO obtains the most NCM, and basically far more than
other algorithms except for scene (j) and (o), in which RIFT+ has more NCM.
The cases in the two scenes have already been discussed in the former section.
However, the proposed MS-HLMO still obtains a considerable number of correct
matches that are sufficient for accurate registration, and without the
interference of orientation, MS-HLMO+ obtains more matches.

SIFT, SAR-SIFT, PSO-SIFT, and PIIFD fail in registration in most multi-source
scenes, which cannot well cope with multi-source remote sensing data images in
various situations. Note that some of the test data selected in this paper are
also used in MS-PIIFD, and the registration is successful when the number of
Harris corners is set to 1000. However, when the number of Harris corners is
increased to 2000 in this experiment, the effect of MS-PIIFD is greatly
reduced, and it even fails. This is because the proportion of correct matching
of feature points obtained by MS-PIIFD before going through FSC is low, and
more outliers may be introduced after increasing the number of feature points,
which makes it difficult to calculate a consistent transformation model and
leads to performance degradation. It indicates that although MS-PIIFD is able
to handle multi-scale and multi-mode image registration, its feature
robustness and stability are limited. Under the same conditions, MS-HLMO
obtains a large number of correct matches, which indicates that its feature
robustness and stability are relatively high. The RIFT+ works well in scenes
with only intensity difference, but the NCMs are relatively low. RIFT+ cannot
cope with the images containing rotation and scale differences.

The last scene is a single-source image with only rotation. This scene proves
that MS-HLMO can not only deal with multi-source images effectively but also
has a good effect on single-source images, which still performs better than
other algorithms. Among all the competitive methods, MS-HLMO is the only
algorithm that does not fail in all scenes and obtains sufficient feature
matches. In the scenes without obvious image rotation difference, MS-HLMO+
achieves the most feature point matching, far outperforming other algorithms.
Through the above experiments, the ability of MS-HLMO is verified from a
detailed and comprehensive perspective, and it can effectively handle the
registration task of multi-source remote sensing images. Compared with the
results of the existing effective algorithms, the proposed MS-HLMO has much
more significant effect than other algorithms in most cases. MS-HLMO+ better
deals with the situation without obvious rotation difference, which is also
the majority of practical multi-source image registration tasks. MS-HLMO and
its simplified version MS-HLMO+ have obvious advantages in robust feature
extraction and matching.

Fig.14 illustrates the registration results of the 17 scenes, in which the
transformation parameters is obtained based on the feature matching shown in
Fig.13. For each image pair, one with the higher resolution is taken as the
reference, and the other is aligned with it through spatial affine
transformation and data interpolation. Image pairs with large spatial offset
are shown in a fusion form, and the areas with prominent structures in the
overlapping parts are selected and shown in an alternation form. All the
operations are performed autonomously without human intervention.

It is observed that most of the image pairs are well registered, the outlines
and textures of the ground covers stay continuous and consistent. The
deviations are mostly within one pixel. The only obvious deviations are found
in some areas of scene (g) and (i) shown in Fig.14. This phenomenon does not
result from inaccurate feature extraction or matching but is caused by the
characteristics of the data. Scene (g) and (i) are ground-based images with
relatively close imaging distance, in which there are obvious viewpoint
differences. In addition, nonrigid geometric distortions exist in the infrared
images, which may be caused by the imaging capability of the sensors. Due to
the above reasons, more complex spatial differences exist between the image
pairs, which is difficult to be solved with only an affine model globally
transforming the images. In practice, this problem can be solved by using more
sophisticate nonlinear transformation models, or by dividing the image into
blocks for adaptive local registration. Apart from that, most areas is well
registered, where the objects are accurately aligned.

##  IV Conclusions

In this paper, an image registration algorithm has been proposed for multi-
source remote sensing image registration. The proposed MS-HLMO utilizes
Partial Main Orientation Map (PMOM) based on Average Squared Gradient (ASG) as
the basic feature map, a generalized GLOH-like (GGLOH) descriptor structure,
and multi-scale feature extraction and matching strategy, which has the
benefit of high robustness to multimodal properties and effectively deals with
the key problems of intensity distortion, rotation, and scale differences in
registration of multi-source images. Comprehensive experiments on 17 datasets
of multi-source remote sensing scenes reveal the effectiveness and
generalization of the proposed MS-HLMO and its simplified version MS-HLMO+
compared with other state-of-the-art registration methods.

## References

  * [1] M. Zhang, W. Li, Q. Du, L. Gao, and B. Zhang, “Feature extraction for classification of hyperspectral and lidar data using patch-to-patch cnn,” _IEEE Transactions on Cybernetics_ , vol. 50, no. 1, pp. 100–111, 2020.
  * [2] X. Zhao, R. Tao, W. Li, H.-C. Li, Q. Du, W. Liao, and W. Philips, “Joint classification of hyperspectral and lidar data using hierarchical random walk and deep cnn architecture,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 58, no. 10, pp. 7355–7370, 2020.
  * [3] C. Liu, R. Tao, W. Li, M. Zhang, W. Sun, and Q. Du, “Joint classification of hyperspectral and multispectral images for mapping coastal wetlands,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 982–996, 2020.
  * [4] M. Zhang, W. Li, R. Tao, and S. Wang, “Transfer learning for optical and sar data correspondence identification with limited training labels,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 1545–1557, 2020.
  * [5] Y. Gao, W. Li, M. Zhang, J. Wang, W. Sun, R. Tao, and Q. Du, “Hyperspectral and multispectral classification for coastal wetland using depthwise feature interaction network,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, 2021.
  * [6] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical Image Analysis_ , vol. 42, pp. 60–88, 2017.
  * [7] B. Zitova and J. Flusser, “Image registration methods: a survey,” _Image and Vision Computing_ , vol. 21, no. 11, pp. 977–1000, 2003.
  * [8] W. K. Pratt, “Correlation techniques of image registration,” _IEEE Transactions on Aerospace and Electronic Systems_ , no. 3, pp. 353–358, 1974.
  * [9] E. De Castro and C. Morandi, “Registration of translated and rotated images using finite fourier transforms,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , no. 5, pp. 700–703, 1987.
  * [10] P. Viola and W. M. Wells III, “Alignment by maximization of mutual information,” _International Journal of Computer Vision_ , vol. 24, no. 2, pp. 137–154, 1997.
  * [11] J. Le Moigne, W. J. Campbell, and R. F. Cromp, “An automated parallel image registration technique based on the correlation of wavelet features,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 40, no. 8, pp. 1849–1864, 2002.
  * [12] J. P. Pluim, J. A. Maintz, and M. A. Viergever, “Mutual-information-based registration of medical images: a survey,” _IEEE Transactions on Medical Imaging_ , vol. 22, no. 8, pp. 986–1004, 2003.
  * [13] A. Mahmood and S. Khan, “Correlation-coefficient-based fast template matching through partial elimination,” _IEEE Transactions on Image Processing_ , vol. 21, no. 4, pp. 2099–2108, 2011.
  * [14] F. P. Oliveira and J. M. R. Tavares, “Medical image registration: a review,” _Computer Methods in Biomechanics and Biomedical Engineering_ , vol. 17, no. 2, pp. 73–93, 2014.
  * [15] X. Tong, Z. Ye, Y. Xu, S. Liu, L. Li, H. Xie, and T. Li, “A novel subpixel phase correlation method using singular value decomposition and unified random sample consensus,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 8, pp. 4143–4156, 2015.
  * [16] X. Zhang, C. Leng, Y. Hong, Z. Pei, I. Cheng, and A. Basu, “Multimodal remote sensing image registration methods and advancements: A survey,” _Remote Sensing_ , vol. 13, no. 24, p. 5128, 2021.
  * [17] D. G. Lowe, “Object recognition from local scale-invariant features,” in _Proceedings of the Seventh IEEE International Conference on Computer Vision_ , vol. 2. IEEE, 1999, pp. 1150–1157.
  * [18] D. G. Lowe, “Local feature view clustering for 3d object recognition,” in _Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001_ , vol. 1. IEEE, 2001, pp. I–I.
  * [19] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International Journal of Computer Vision_ , vol. 60, no. 2, pp. 91–110, 2004\.
  * [20] Y. Ke and R. Sukthankar, “Pca-sift: A more distinctive representation for local image descriptors,” in _Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004._ , vol. 2. IEEE, 2004, pp. II–II.
  * [21] K. Mikolajczyk and C. Schmid, “A performance evaluation of local descriptors,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 27, no. 10, pp. 1615–1630, 2005.
  * [22] H. Bay, T. Tuytelaars, and L. Van Gool, “Surf: Speeded up robust features,” in _European Conference on Computer Vision_. Springer, 2006, pp. 404–417.
  * [23] J.-M. Morel and G. Yu, “Asift: A new framework for fully affine invariant image comparison,” _SIAM Journal on Imaging Sciences_ , vol. 2, no. 2, pp. 438–469, 2009.
  * [24] A. Sedaghat, M. Mokhtarzade, and H. Ebadi, “Uniform robust scale-invariant feature matching for optical remote sensing images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 49, no. 11, pp. 4516–4527, 2011.
  * [25] F. Dellinger, J. Delon, Y. Gousseau, J. Michel, and F. Tupin, “Sar-sift: a sift-like algorithm for sar images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 1, pp. 453–466, 2014.
  * [26] A. Sedaghat and H. Ebadi, “Remote sensing image matching based on adaptive binning sift descriptor,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 53, no. 10, pp. 5283–5293, 2015.
  * [27] T. Lindeberg, “Scale-space theory: A basic tool for analyzing structures at different scales,” _Journal of Applied Statistics_ , vol. 21, no. 1-2, pp. 225–270, 1994.
  * [28] H. P. Moravec, “Obstacle avoidance and navigation in the real world by a seeing robot rover,” Ph.D. dissertation, Stanford University, 1980.
  * [29] C. G. Harris, M. Stephens, _et al._ , “A combined corner and edge detector.” in _Alvey Vision Conference_ , vol. 15, no. 50. Citeseer, 1988, pp. 10–5244.
  * [30] N. Dalal and B. Triggs, “Histograms of oriented gradients for human detection,” in _2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)_ , vol. 1. IEEE, 2005, pp. 886–893.
  * [31] E. Rosten, R. Porter, and T. Drummond, “Faster and better: A machine learning approach to corner detection,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 32, no. 1, pp. 105–119, 2008.
  * [32] M. Calonder, V. Lepetit, M. Ozuysal, T. Trzcinski, C. Strecha, and P. Fua, “Brief: Computing a local binary descriptor very fast,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 34, no. 7, pp. 1281–1298, 2011.
  * [33] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 2564–2571.
  * [34] S. Leutenegger, M. Chli, and R. Y. Siegwart, “Brisk: Binary robust invariant scalable keypoints,” in _2011 International Conference on Computer Vision_. IEEE, 2011, pp. 2548–2555.
  * [35] A. Alahi, R. Ortiz, and P. Vandergheynst, “Freak: Fast retina keypoint,” in _2012 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2012, pp. 510–517.
  * [36] P. F. Alcantarilla, A. Bartoli, and A. J. Davison, “Kaze features,” in _European Conference on Computer Vision_. Springer, 2012, pp. 214–227.
  * [37] Y. Ye, L. Bruzzone, J. Shan, F. Bovolo, and Q. Zhu, “Fast and robust matching for multimodal remote sensing image registration,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 57, no. 11, pp. 9059–9070, 2019.
  * [38] J. Chen, J. Tian, N. Lee, J. Zheng, R. T. Smith, and A. F. Laine, “A partial intensity invariant feature descriptor for multimodal retinal image registration,” _IEEE Transactions on Biomedical Engineering_ , vol. 57, no. 7, pp. 1707–1718, 2010.
  * [39] W. Ma, Z. Wen, Y. Wu, L. Jiao, M. Gong, Y. Zheng, and L. Liu, “Remote sensing image registration with modified sift and enhanced feature matching,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 1, pp. 3–7, 2016\.
  * [40] Y. Ye, J. Shan, L. Bruzzone, and L. Shen, “Robust registration of multimodal remote sensing images based on structural similarity,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 55, no. 5, pp. 2941–2958, 2017.
  * [41] Y. Ye, J. Shan, S. Hao, L. Bruzzone, and Y. Qin, “A local phase based invariant feature for remote sensing image matching,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 142, pp. 205–221, 2018.
  * [42] J. Li, Q. Hu, and M. Ai, “Rift: Multi-modal image matching based on radiation-variation insensitive feature transform,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 3296–3310, 2019.
  * [43] C. Gao and W. Li, “Multi-scale piifd for registration of multi-source remote sensing images,” _Journal of Beijing Institute of Technology_ , vol. 30, no. 2, pp. 113–124, 2021.
  * [44] C. Gao and W. Li, “Multi-scale harris-piifd features for registration of visible and infrared images,” in _2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS_. IEEE, 2021, pp. 5437–5440.
  * [45] M. Kass and A. Witkin, “Analyzing oriented patterns,” _Computer Vision, Graphics, and Image Processing_ , vol. 37, no. 3, pp. 362–385, 1987.