# Separability in Morse local-to-global groups

Lawk Mineh and Davide Spriano

###### Abstract.

We show that in a Morse local-to-global group where stable subgroups are
separable, the product of any stable subgroups is separable. As an
application, we show that the product of stable subgroups in virtually special
groups is separable.

##  1\. Introduction

Given a group $G$, we can equip it with the _profinite topology_ , whose basic
open subsets are cosets of finite index subgroups of $G$. A subset of $G$ is
said to be _separable_ if it is closed in the profinite topology on $G$. The
group $G$ is called _residually finite_ if the trivial subgroup is separable
in $G$.

Knowing that particular subsets of groups are separable often gives useful
information about the group. For example, in a finitely presented group,
separability of a finitely generated subgroup gives a solution to the
membership problem in that subgroup. In a geometric setting, separability
properties of the fundamental group of a space correspond to desirable lifting
properties of that space: immersed subcomplexes of a complex $X$ may be
promoted to embedded ones in a finite sheeted cover of $X$, provided that
their corresponding subgroups are separable in $\pi_{1}X$. For an example
involving subsets rather than subgroups, it has been proven that if $X$ is a
nonpositively curved cube complex in which every double coset of hyperplane
stabilisers is separable in $\pi_{1}X$, then $X$ has a finite sheeted special
cover [HW08].

It is a difficult problem to show that a given subset of a group is separable,
especially when one is only given some geometric data about the group. For
instance, even the question of whether hyperbolic groups are residually finite
is a long-standing open problem. It is known that all hyperbolic groups are
residually finite if and only if every quasiconvex subgroup is separable in
every hyperbolic group [AGM09]. Minasyan showed that if $G$ is a hyperbolic
group in which every quasiconvex subgroup is separable, the setwise product of
any finite number of quasiconvex subgroups is also separable in $G$ [Min06],
extending a result of Ribes and Zalesskii, who proved the same result in the
case $G$ is free [RZ93]. Recently, the first author and Minasyan provided
generalisation of this result in the setting of relatively hyperbolic groups
[MM22]. In this paper, we will provide another natural generalisation of this
product separability result to the class of groups with the _Morse local-to-
global_ (MLTG) property.

Introduced in [RST22], the MLTG property roughly speaking requires that
quasigeodesics with hyperbolic-like properties behave similarly to
quasigeodesics in hyperbolic spaces. Consider the following two perspectives
on hyperbolic spaces. The first involves Morse geodesics: we say that a
quasigeodesic is _Morse_ if any other quasigeodesic with the same endpoints
stays uniformly finite Hausdorff-distance from it (see Definition 2.3). It is
a well-known fact that every quasigeodesic in a hyperbolic space satisfies the
Morse property, and moreover that a space is hyperbolic if and only if all of
its geodesics are Morse [Bon96]. This motivates the study of Morse
quasigeodesics in spaces that are not hyperbolic, an approach that has been
successful in understanding the properties of spaces up to quasi-isometry
[CH17, IM21, GKLS21, CRSZ22, QR22]. On the other hand, Gromov showed that a
space is hyperbolic if and only if all of it local-quasigeodesics, i.e. paths
that are quasigeodesics on every subpath of a certain length (see Definition
2.2) are globally quasigeodesics. The Morse local-to-global property puts
these two perspectives together and prescribes that all paths that are locally
Morse quasigeodesics are globally Morse quasigeodesics.

In groups with the MLTG property, elements acting on Morse geodesics behave
“as they should”. For instance, it is appealing to think that given two
independent infinite order elements with Morse axes, then it is possible to
use a ping-pong argument to generate a free subgroup using these elements. In
general finitely generated groups this is not true, and we require the MLTG
property in order for it to be satisfied. As the above suggests, the failure
of the MLTG property seems to imply some pathological behaviour; the only
known examples of groups without the MLTG property are not finitely
presentable. On the other hand, many well-behaved classes of groups, such as
3-manifold groups, CAT(0) groups, and mapping class groups are known to
satisfy the MLTG property.

Our main theorem is concerned with separability of products of _stable
subgroups_. Stable subgroups were introduced in [DT15] where they showed that
a the convex cocompact subgroups of the mapping class groups are precisely the
stable ones. For infinite cyclic subgroups the notion of stability and fixing
a Morse quasigeodesic agree, and in general stable subgroups present many
properties akin to quasiconvex subgroups of hyperbolic groups.

######  Theorem 1.1.

Let $G$ be a finitely generated group with the Morse local-to-global property,
and suppose that any stable subgroup of $G$ is separable. Then the product of
any stable subgroups of $G$ is separable.

Recall that a group is _LERF_ (locally extended residually finite) if every
finitely generated subgroup is separable. The following statement may be of
more general interest, for instance as a criterion for showing that a given
group is not LERF. As stable subgroups are always finitely generated (see
Lemma 2.9), the hypotheses are stronger than the above theorem.

######  Corollary 1.2.

Let $G$ be a finitely generated LERF group with the Morse local-to-global
property. Then the product of any stable subgroups of $G$ is separable.

A group is _virtually special_ if it is has a finite index subgroup that is
the fundamental group of a special cube complex. Triple cosets of convex
subgroups in virtually special groups are known to be separable [She23]. We
extend this result to arbitrary products of stable subgroups, which are
quasiconvex.

######  Corollary 1.3.

Let $G$ be a virtually special group. Then the product of any stable subgroups
of $G$ is separable.

_Strongly quasiconvex_ subgroups, also known as _Morse_ subgroups, were
introduced independently by Genevois and Tran [Tra19, Gen19], and Tran showed
that a subgroup of a finitely generated group is stable if and only if it is
strongly quasiconvex and hyperbolic [Tra19, Proposition 4.3]. In the case of
right-angled Artin groups, which contain many stable subgroups [KMT17], we can
use [RST23, Corollary 7.4] to deduce the following.

######  Corollary 1.4.

Suppose that $\Gamma$ is a finite connected graph, and let $A_{\Gamma}$ be the
associated right-angled Artin group. Then the product of any strongly
quasiconvex subgroups of $A_{\Gamma}$ is separable.

###  1.1. Acknowledgments

The authors are grateful to the support of the CRM in Montréal where part of
this work was completed. We would like to thank Sam Shepherd for interesting
conversations on the topic.

##  2\. Preliminaries

In this paper, we will restrict our attention to Cayley graphs of groups.

######  Definition 2.1.

Let $G$ be a group and $S$ a generating set for $G$. The _Cayley graph_ of $G$
with respect to $S$ is the graph $\mathrm{Cay}(G,S)$ with vertex set $G$ and
elements $g,h$ connected by an edge if either $gh^{-1}\in S$ or $hg^{-1}\in
S$.

We equip the set of vertices of a graph with the metric induced by declaring
all of its edges to have length one. Given a path $\gamma$ of a graph, we will
denote its length (i.e. number of edges) by $\ell(\gamma)$. Given metric
spaces $(X,d_{X})$ and $(Y,d_{Y})$, a _$(\lambda,c)$ –quasi-isometric
embedding_ is a map $f\colon X\to Y$ such that the following holds for any
pair $x,y\in X$.

$\frac{1}{\lambda}d_{Y}(f(x),f(y))-c\leq d_{X}(x,y)\leq\lambda
d_{Y}(f(x),f(y))+c.$

A _$(\lambda,c)$ –quasigeodesic_ is a $(\lambda,c)$–quasi isometric embedding
of an interval $I\subset\mathbb{R}$.

The main geometric definition of the paper is the Morse local-to-global
property. To define it, we need to define the Morse property and what is means
for a given property of a path to be local.

######  Definition 2.2 (Local property).

A path $\gamma\colon I\to X$ is said to _$L$ –locally_ satisfy a property $P$
if each subpath of the form $\gamma|_{[t_{1},t_{2}]}$ with $t_{2}-t_{1}\leq L$
satisfies $P$. When a path $\gamma$ is $L$–locally a
$(\lambda,c)$–quasigeodesic, we say that $\gamma$ is a $(L;\lambda,c)$–local
quasigeodesic.

######  Definition 2.3 (Morse quasigeodesic).

Let $M\colon\mathbb{R}_{\geq 1}\times\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq
0}$ be a non-decreasing function. A quasigeodesic $\gamma\colon I\to X$ is
_$M$ –Morse_ if for any $(\lambda,c)$–quasigeodesic segment
$\eta\colon[a,b]\to X$ such that $\eta(a)=\gamma(t_{1}),\eta(b)=\gamma(t_{2})$
we have

$d_{H}(\gamma|_{[t_{1},t_{2}]},\eta)\leq M(\lambda,c),$

where $d_{H}$ denotes the Hausdorff-distance. We say that $\gamma$ is a
$(M;\lambda,c)$–Morse quasigeodesic.

Morse geodesics satisfy a thin-triangles condition, similar to geodesics is in
a hyperbolic metric space.

######  Lemma 2.4 ([Liu21, Lemma 3.6]).

Let $X$ be a geodesic metric space and suppose that $p$ and $q$ are $M$–Morse
geodesics with $p_{-}=q_{-}$. There is a constant $\delta=\delta(M)\geq 0$
such that for any geodesic $r$ with endpoints $r_{-}=p_{+}$ and $r_{+}=q_{+}$,
the geodesic triangle $p\cup q\cup r$ is $\delta$-thin.

######  Definition 2.5 (Local Morse quasigeodesic).

We say that a path is an _$(L;M;\lambda,c)$ –local Morse quasigeodesic_ if it
is $L$–locally an $M$–Morse $(\lambda,c)$–quasigeodesic.

######  Definition 2.6 (MLTG property).

We say that a metric space $X$ satisfies the _Morse local-to-global property_
, for short MLTG property, if for any choice of Morse gauge $M$ and constants
$\lambda\geq 1,c\geq 0$ there exist a scale $L>0$, a Morse gauge $M^{\prime}$
and constants $\lambda^{\prime}\geq 1,c^{\prime}\geq 0$ such that every
$(L;M;\lambda,c)$–local Morse quasigeodesic is a
$(M^{\prime};\lambda^{\prime},c^{\prime})$–Morse quasigeodesic.

The strength of the MLTG property is that it allows us to draw global
conclusions from local conditions, as the next lemma shows.

######  Lemma 2.7.

Let $p=p_{1},\dots,p_{n}$ be a concatenation of $M$–Morse geodesics in a MLTG
space $X$ and let $a_{i}$ and $a_{i+1}$ be the ordered endpoints of $p_{i}$.
For each $\varepsilon>0$ there are constants $B\geq 0,\lambda\geq 1,c\geq 0$,
and a gauge $N$ (all depending only on $M$ and $\varepsilon$) such that if we
have $\ell(p_{i})>B$ for all $i=2,\dots,n-1$ and
$\langle(a_{i-1}),(a_{i+1})\rangle_{a_{i}}\leq\varepsilon$ for all
$i=2,\dots,n$, then $p$ is a $(N;\lambda,c)$–Morse quasigeodesic.

###### Proof.

We will show that there exists $M^{\prime}$ depending only on $M$ and
$\varepsilon$ such that $p$ is locally a $(M^{\prime};1,2\epsilon)$–Morse
quasigeodesic, and then we will choose an appropriate $B$ to use the MLTG
property. Given the existence of such $M^{\prime}$, the MLTG property gives us
a Morse gauge $N$ and constants $\lambda\geq 1,c\geq 0,L\geq 0$ such that any
$(L;M^{\prime};1,2\varepsilon)$–local Morse quasigeodesic is also a
$(N;\lambda,c)$–Morse quasigeodesic.

We start with the quasigeodesic claim. Take $B\geq L+\varepsilon$ and observe
that if we consider two points $x,y$ at parameterised distance at most $L$,
they either lie on the same segment $p_{i}$ (which is geodesic), or on two
consecutive segments $p_{i-1}$ and $p_{i}$. In the latter case, we have

$\langle x,y\rangle_{a_{i}}\leq\langle
a_{i-1},a_{i+1}\rangle_{a_{i}}\leq\varepsilon$

which means

$d(x,y)+2\epsilon\geq d(a_{i},x)+d(a_{i},y)=\ell(p_{i-1}|_{[x,a_{i}]}\ast
p_{i}|_{[a_{i},y]}).$

Thus, $p$ is an $(L;1,2\epsilon)$–local quasigeodesic. A similar computation
to that above (with $x=a_{i-1},y=a_{i+1}$) shows that

$d(a_{i-1},a_{i+1})\geq 2B-2\varepsilon>L$

where the last inequality comes from the choice of $B$. Now, applying [RST22,
Lemma 2.15] to each concatenation $p_{i}\ast p_{i+1}$, we obtain some
$M^{\prime}$ (depending only on $M$ and $\varepsilon$) such that $p$ is an
$(L;M^{\prime};1,2\epsilon)$–local Morse quasigeodesic. Now applying the MLTG
property shows that $p$ is $(N;\lambda,c)$–Morse quasigeodesic. ∎

The notion of stability – in some sense – generalises the notion of having the
Morse property from quasigeodesics to arbitrary subgroups.

######  Definition 2.8 (Stable subgroup).

Let $G$ be a group with finite generating set $S$. A subgroup $H\leqslant G$
is called _stable_ if there is a Morse gauge $M$ and constant $\mu\geq 0$ such
that any geodesic in $\mathrm{Cay}(G,S)$ with endpoints in $H$ is $M$–Morse
and lies in the $\mu$-neighbourhood of $H$.

An immediate consequence of the definition is that a stable subgroup of a
finitely generated group is undistorted. We note that while the gauge $M$ and
constant $\mu$ in the above depend on the choice of generating set $S$, the
property of being stable does not (see, for example, [DT15, Lemma 3.4]). We
recall some basic properties of stability, which follow from [DT15].

######  Lemma 2.9.

Let $G$ be a finitely generated group and suppose $H\leqslant G$ is stable.
Then the following are true:

  1. (1)

if $K\leqslant H$ has finite index in $H$, then $K$ is stable;

  2. (2)

if $g\in G$, then $gHg^{-1}$ is stable;

  3. (3)

$H$ is finitely generated and undistorted in $G$;

  4. (4)

$H$ is hyperbolic.

######  Remark 2.10.

In Item 1 we require $K$ has finite index in $H$. However, note that for _any_
subgroup $K\leqslant H$ the following weak form of stability holds, namely
that any geodesic with endpoints on $K$ is $M$–Morse. That is to say, a
subgroup of an $(M,\mu)$–stable group if $(M,\mu^{\prime})$–stable, if we
allow $\mu^{\prime}=\infty$. The finite index assumption guarantees
$\mu^{\prime}<\infty$.

The following lemma tells us that stable subgroups cannot be “too parallel”
away from their intersection. More precisely, that there is a uniform upper
bound on the Gromov product of elements from one subgroup when taken with
minimal length coset representatives of the other.

######  Lemma 2.11.

Let $G$ be a group with finite generating set $S$, and suppose that $H$ and
$K$ are $(M,\mu)$-stable subgroups of $G$. There is a constant
${\rho}={\rho}(M,\mu,S)\geq 0$ such that if $h\in H$ is a shortest (with
respect to $S$) representative of its right coset $(H\cap K)h$, then for any
$k\in K$, we have $\langle h,k\rangle_{1}\leq{\rho}$.

###### Proof.

Suppose for a contradiction that we can find elements $h\in H,k\in K$ such
that $h$ is a shortest coset representative of $h(H\cap K)$ and $\langle
h,k\rangle_{1}$ is arbitrarily large. Since $H$ and $K$ are stable, any choice
of geodesics $p=[1,h]$ and $q=[1,k]$ are $M$–Morse and lie in a
$\mu$-neighbourhood of $H$ and $K$ respectively.

Let $a_{1},\dots,a_{n}$ be the vertices of $p$ with $d_{S}(1,a_{i})\leq\langle
h,k\rangle_{1}$. The assumption that $\langle h,k\rangle_{1}$ can be taken to
be arbitrarily large means that $n$ can be taken to be arbitrarily large.
Corresponding to each vertex $a_{i}$, there is $v_{i}\in H$ such that
$d_{S}(a_{i},v_{i})\leq\mu$ by stability of $H$. Moreover, by Lemma 2.4 there
is $\delta=\delta(M)\geq 0$ such that $d_{S}(a_{i},q)\leq\delta$ for each
$i=1,\dots,n$. By stability of $K$ and the triangle inequality, therefore, we
obtain that $d_{S}(v_{i},K)\leq 2\mu+\delta$ for each $i=1,\dots,n$. Note that
since $p$ is geodesic

(1) $d_{S}(1,v_{i})\leq i-\mu\qquad\text{ and }\qquad d_{S}(v_{i},h)\leq
d_{S}(1,h)-i-\mu$

for each $i=1,\dots,n$.

For each $i=1,\dots,n$, let $g_{i}$ be the shortest element of $G$ with
respect to $S$ such that $v_{i}g_{i}\in K$, so $\left\lvert
g_{i}\right\rvert_{S}\leq 2\mu+\delta$. Let $N$ be the number of elements in
the ball of radius $2\mu+\delta$ about the identity in $\mathrm{Cay}(G,S)$.
Taking $n$ to be sufficiently large with respect to $N$ and $\mu$, there must
be some pair $(i,j)$ with $g_{i}=g_{j}$ satisfying $j-i>2\mu$. Then equation
(1) gives

(2) $d_{S}(1,v_{i})+d_{S}(v_{j},h)<d_{S}(1,h).$

But then $v_{j}g_{j}(v_{i}g_{i})^{-1}=v_{j}v_{i}^{-1}\in H\cap K$, as
$v_{i},v_{j}\in H$ and $v_{i}g_{i},v_{j}g_{j}\in K$. Moreover

$\displaystyle d_{S}(v_{j}v_{i}^{-1},h)$ $\displaystyle\leq
d_{S}(v_{j}v_{i}^{-1},v_{j})+d_{S}(v_{j},h)$
$\displaystyle=d_{S}(1,v_{i})+d_{S}(v_{j},h)<d_{S}(1,h)$

where the last inequality is an application of (2). It follows that
$v_{i}v_{j}^{-1}h\in(H\cap K)h$ and

$\left\lvert
v_{i}v_{j}^{-1}h\right\rvert_{S}=d_{S}(v_{j}v_{i}^{-1},h)<d_{S}(h,1)=\left\lvert
h\right\rvert_{S}$

which contradicts the fact that $h$ is a minimal length representative of its
$(H\cap K)$-coset. Thus, there must be an upper bound on the Gromov product. ∎

We finish this section by recalling the key property that we need for our
proof, namely that the MLTG property allows a ping-pong type argument for
stable subgroups.

######  Proposition 2.12 ([RST22, Theorem 3.1]).

Let $G$ be a group with finite generating set $S$ and suppose $G$ has the
Morse local-to-global property. Let $Q,R\leqslant G$ be stable subgroups of
$G$ with corresponding Morse gauge $M$ and constant $\mu\geq 0$. There is a
constant $C=C(M,\mu,S)\geq 0$ such that the following is true.

Let $Q^{\prime}\leqslant Q$ and $R^{\prime}\leqslant R$ be subgroups such that
$Q^{\prime}\cap R^{\prime}=Q\cap R$ and $\left\lvert g\right\rvert_{S}\geq C$
for each $g\in(Q^{\prime}\cup R^{\prime})\setminus(Q^{\prime}\cap
R^{\prime})$. Then $\langle Q^{\prime},R^{\prime}\rangle\cong
Q^{\prime}\ast_{Q^{\prime}\cap R^{\prime}}R^{\prime}$. Moreover, if
$Q^{\prime}$ and $R^{\prime}$ are finitely generated and undistorted in $G$,
then $\langle Q^{\prime},R^{\prime}\rangle$ is stable.

##  3\. Separability of products

In this section we will prove the main theorem. We start by recalling some
basic properties of separable subgroups.

######  Remark 3.1.

If $U\subseteq G$ is a separable subset of $G$, then $U^{-1},gU$, and $Ug$ are
separable for any $g\in G$.

######  Remark 3.2.

Let $H_{1},\dots,H_{n}\leqslant G$ be subgroups of $G$ and let
$a_{0},\dots,a_{n}\in G$. Observe that

$a_{0}H_{1}a_{1}\dots a_{n-1}H_{n}a_{n}=H_{1}^{a_{0}}H_{2}^{a_{0}a_{1}}\dots
H_{n}^{a_{0}\dots a_{n-1}}a_{0}\dots a_{n},$

which is a translate of a product of conjugates of the subgroups
$H_{1},\dots,H_{n}$. The set $a_{0}H_{1}a_{1}\dots a_{n-1}H_{n}a_{n}$ is thus
separable if and only if the product of subgroups $H_{1}^{a_{0}}\dots
H_{n}^{a_{0}\dots a_{n-1}}$ is separable.

In particular, suppose there is $n\in\mathbb{N}$ such that any product of $n$
stable subgroups are separable in $G$, and suppose $H_{1},\dots,H_{n}$ are
stable subgroups of $G$. Lemma 2.9(2) gives that $H_{i}^{a_{0}\dots a_{i-1}}$
is stable for each $1\leq i\leq n$, so that the set $H_{1}^{a_{0}}\dots
H_{n}^{a_{0}\dots a_{n-1}}$ is a product of $n$ stable subgroups. By the
observation above we may conclude that the set $a_{0}H_{1}a_{1}\dots
a_{n-1}H_{n}a_{n}$ is separable in this situation.

In order to exploit the geometric edge given by the MLTG property, it is
useful to choose coset representative that are geometrically meaningful, which
we can do by the following remark.

######  Remark 3.3.

Suppose that $S$ is a generating set for group $G$ and let
$H_{1},\dots,H_{n}\leqslant G$ be subgroups of $G$. Given elements $x_{1}\in
H_{1},\dots,x_{n}\in H_{n}$, there are elements $y_{1}\in H_{1},\dots,y_{n}\in
H_{n}$ such that $x_{1}\dots x_{n}=y_{1}\dots y_{n}$ and $\left\lvert
y_{i}\right\rvert_{S}$ is minimal among elements of the coset $(H_{i-1}\cap
H_{i})y_{i}$ for each $1<i\leq n$.

Indeed, there is $y_{n}\in H_{n}$ and $z_{n}\in H_{n}\cap H_{n-1}$ such that
$x_{n}=z_{n}y_{n}$ and $\left\lvert y_{n}\right\rvert_{S}$ is minimal among
elements of $(H_{n-1}\cap H_{n})x_{n}=(H_{n-1}\cap H_{n})y_{n}$. Similarly
there is $y_{n-1}\in H_{n_{1}}$ and $z_{n-1}\in H_{n-1}\cap H_{n-2}$ such that
$x_{n-1}z_{n}=z_{n-1}y_{n-1}$ and $y_{n-1}$ is a shortest representative of
$(H_{n-2}\cap H_{n-1})x_{n-1}z_{n}=(H_{n-2}\cap H_{n-1})y_{n-1}$. We can
proceed by finite induction to find elements $y_{2}\in H_{2},\dots,y_{n2}\in
H_{n}$ and $z_{2}\in H_{2}\cap H_{3},\dots,z_{n}\in H_{n-1}\cap H_{n}$ with
the properties described above. Setting $y_{1}=x_{1}z_{2}\in H_{1}$ completes
the observation.

We conclude the section by proving the main theorem of the paper and the
related corollaries.

###### Proof of Theorem 1.1.

We proceed by induction on the number $n$ of stable subgroups. The case $n=1$
is exactly the hypothesis that stable subgroups of $G$ are separable, so let
$H_{1},\dots,H_{n}$ be stable subgroups of $G$ with $n>1$ and suppose that the
product of any $n-1$ stable subgroups of $G$ is separable.

Fix a finite generating set $S$ of $G$. By taking maxima of gauges and
constants, we may assume without loss of generality that $H_{1},\dots,H_{n}$
are $(M,\mu)$-stable. Let $\rho=\rho(M,\mu,S)$ be the constant obtained from
Lemma 2.11, and let $C=C(M,\mu,S)$ be the constant of Proposition 2.12. Let
$B,\lambda,c\geq 0$ be the constants and $N$ the Morse gauge obtained from
applying Lemma 2.7 with gauge $M$ and constant $\rho$.

Suppose for a contradiction that the product $H_{1}\dots H_{n}$ is not
separable, so that there is some $g\notin H_{1}\dots H_{n}$ belonging to the
profinite closure of $H_{1}\dots H_{n}$. This means that $g$ is contained in
every separable subset containing $H_{1}\dots H_{n}$. For ease of reading, we
will write $Q=H_{1},R=H_{2},$ and $T_{i}=H_{i+2}$ whenever $1\leq i\leq
s=n-2$. By hypothesis $Q$ and $R$ are separable, and thus their intersection
$I=Q\cap R$ is also. Let $\\{N_{i}\\}_{i\in\mathbb{N}}$ be an enumeration of
the finite index subgroups of $G$ containing $I$, and note that
$I=\bigcap_{i\in\mathbb{N}}N_{i}$ as $I$ is separable. For each $i$, we write

$N^{\prime}_{i}=\bigcap_{j=1}^{i}N_{j}$

so that $\\{N^{\prime}_{i}\\}_{i\in\mathbb{N}}$ is a sequence of nested finite
index subgroups of $G$ containing $I$ whose intersection is equal to $I$.

For each $i\in\mathbb{N}$, we define the set

$K_{i}=Q\langle Q^{\prime}_{i},R^{\prime}_{i}\rangle RT_{1}\dots T_{s}$

where $Q^{\prime}_{i}=N^{\prime}_{i}\cap Q\leqslant_{f}Q$ and
$R^{\prime}_{i}=N^{\prime}_{i}\cap R\leqslant_{f}R$. Note that $I\subseteq
N^{\prime}_{i}$ for each $i\in\mathbb{N}$, so that $Q^{\prime}_{i}\cap
R^{\prime}_{i}=I$. It is immediate from the definition that $K_{i}\supseteq
QRT_{1}\dots T_{s}$ for each $i\in\mathbb{N}$. Our aim is to show that for
sufficiently large $i$, the set $K_{i}$ is separable and excludes the element
$g$.

Let us first show that $K_{i}$ is separable when $i$ is large. Indeed, for a
given $i$, let $x_{1},\dots,x_{a}$ be left coset representatives for
$Q^{\prime}_{i}$ in $Q$ and $y_{1},\dots,y_{b}$ be right coset representatives
for $R^{\prime}_{i}$ in $R$. Then we have

$K_{i}=\bigcup^{a}_{j=1}\bigcup^{b}_{k=1}x_{j}\langle
Q^{\prime}_{i},R^{\prime}_{i}\rangle y_{k}T_{1}\dots T_{s}.$

Since $I$ is the intersection of all $N^{\prime}_{i}$, there is an index
$i_{0}\in\mathbb{N}$ such that for any $i\geq i_{0}$, any element $n$ in
$N^{\prime}_{i}$ with $\left\lvert n\right\rvert_{S}\leq C$ belongs to $I$. By
(1) and (3) of Lemma 2.9, $Q^{\prime}_{i}$ and $R^{\prime}_{i}$ are finitely
generated and undistorted, so we may apply Proposition 2.12 to obtain that
$\langle Q^{\prime}_{i},R^{\prime}_{i}\rangle$ is stable. Thus, by Remark 3.2
and the induction hypothesis, $K_{i}$ can be written as a finite union of
separable subsets, and so $K_{i}$ is separable whenever $i\geq i_{0}$.

We now show that there is $i\in\mathbb{N}$ such that $g\notin K_{i}$. As $g$
belongs to the profinite closure of $QRT_{1}\dots T_{s}$ and for each $i\geq
i_{0}$ the set $K_{i}$ is a profinitely closed subset of $G$ containing
$QRT_{1}\dots T_{s}$, $g$ belongs to $K_{i}$ for each $i\geq i_{0}$. That is,
for each $i\geq i_{0}$ we may write

(3) $g=q^{(i)}x_{1}^{(i)}\dots x_{m_{i}}^{(i)}r^{(i)}t_{1}^{(i)}\dots
t_{s}^{(i)}$

for some $m_{i}\in\mathbb{N}$ and $x_{j}^{(i)}\in Q^{\prime}_{i}\cup
R^{\prime}_{i}$ for each $1\leq j\leq m_{i}$, and where $q^{(i)}\in
Q,r^{(i)}\in R,t_{1}^{(i)}\in T_{1},\dots,t_{s}^{(i)}\in T_{s}$.

There are two cases to consider.

_Case 1:_ $\liminf_{i\to\infty}\left\lvert r^{(i)}\right\rvert_{S}<\infty$ or
$\liminf_{i\to\infty}\left\lvert t_{j}^{(i)}\right\rvert_{S}<\infty$ for some
$1\leq j<s$.

We consider only the possibility that $\liminf_{i\to\infty}\left\lvert
r^{(i)}\right\rvert_{S}<\infty$, for the other cases can be dealt with
identically. By definition, there is a subsequence of
$(r^{(i)})_{i\in\mathbb{N}}$ whose terms have length bounded by some fixed
constant. Since there are only finitely many elements of $G$ with any given
length with respect to $S$, we may pass to a further subsequence whose terms
are all equal to a single element $r\in R$. Hence we have

(4) $g\in Q\langle Q^{\prime}_{i},R^{\prime}_{i}\rangle rT_{1}\dots
T_{s}\;\;\textrm{for infinitely many }i\in\mathbb{N}.$

Now by the induction hypothesis and Remark 3.2, the set $QrT_{1}\dots T_{s}$
is separable in $G$. Since $g\notin QRT_{1}\dots T_{s}$, we have $g\notin
QrT_{1}\dots T_{s}$, and so there is $N\lhd_{f}G$ such that $g\notin
NQrT_{1}\dots T_{s}=QNrT_{1}\dots T_{s}$. The subgroup $IN\leqslant_{f}G$ is a
finite index subgroup of $G$ containing $I$, so $N^{\prime}_{i_{1}}\subseteq
IN$ for some $i_{1}\in\mathbb{N}$. Since the sequence of subgroups
$\\{N^{\prime}_{i}\\}_{i\in\mathbb{N}}$ is nested, we have thus shown that

$Q\langle Q^{\prime}_{i},R^{\prime}_{i}\rangle rT_{1}\dots T_{s}\subseteq
QN^{\prime}_{i}rT_{1}\dots T_{s}\subseteq QINrT_{1}\dots T_{s}=QNrT_{1}\dots
T_{s}$

for any $i\geq i_{1}$, where the last equality uses the fact that $QI=Q$.
However, the fact that $g\notin QNrT_{1}\dots T_{s}$ now contradicts the
inclusions of (4), so this case is impossible.

_Case 2:_ $\liminf_{i\to\infty}\left\lvert r^{(i)}\right\rvert_{S}=\infty$ and
$\liminf_{i\to\infty}\left\lvert t_{j}^{(i)}\right\rvert_{S}=\infty$ for all
$1\leq j<s$.

Define
$z_{0}=1,z_{1}=q^{(i)},z_{2}=z_{1}x_{1}^{(i)},\dots,z_{m_{i}+1}=z_{m_{i}}x_{m_{i}}^{(i)},z_{m_{i}+2}=z_{m_{i}+1}r^{(i)},z_{m_{i}+3}=z_{m_{i}+2}t_{1}^{(i)},\dots,z_{m_{i}+2+s}=z_{m_{i}+s}t_{s}^{(i)}$.
For each $0\leq j\leq m_{i}+1+s$, we let $p_{j}$ be a geodesic with
$(p_{j})_{-}=z_{j}$ and $(p_{j})_{+}=z_{j+1}$. Let $p$ be the concatenation
$p_{0}\dots p_{m_{i}+1+s}$ of these paths.

We will use Lemma 2.7 to conclude that the path $p$ is a uniform
quasigeodesic. Assuming this, the fact that $\liminf_{i\to\infty}\left\lvert
r^{(i)}\right\rvert_{S}=\infty$ and $\liminf_{i\to\infty}\left\lvert
t_{j}^{(i)}\right\rvert_{S}=\infty$ means that for sufficiently large $i$, the
distance between the endpoints of $p$ is greater than $\left\lvert
g\right\rvert_{S}$, contradicting the fact that $p$ represents $g$.

Without loss of generality, we may assume $x_{1}^{(i)}\in
R^{\prime}_{i}\setminus Q$ and $x_{m_{i}}^{(i)}\in Q^{\prime}_{i}\setminus R$,
for otherwise we may replace $q^{(i)}$ with $q_{1}^{(i)}=q^{(i)}x_{1}^{(i)}\in
Q$ and eliminate $x_{1}^{(i)}$ from the product (and likewise with $r^{(i)}$
and $x_{m_{i}}^{(i)}$). Further, we may assume by Remark 3.3 that
$x_{1}^{(i)},\dots,x_{m_{i}}^{(i)},$ and $r^{(i)}$ are shortest
representatives of their right $I$-cosets, and in particular
$x_{1}^{(i)},\dots,x_{m_{i}}^{(i)},r^{(i)}\notin I$. Similarly we take
$t_{1}^{(i)}$ to be a shortest representative of $(R\cap T_{1})t_{1}^{(i)}$
and, for $1<i\leq s$, the element $t_{j}^{(i)}$ to be a shortest
representative of $(T_{j-1}\cap T_{j})t_{j}^{(i)}$.

The above paragraph together with Lemma 2.11 shows that

(5) $\langle z_{j-1},z_{j+1}\rangle_{z_{j}}\leq{\rho}\quad\textrm{ for
}j=1,\dots,m_{i}+s.$

We now verify the hypotheses of Lemma 2.7. Each of the geodesic segments
$p_{i}$ represents an element of (a subgroup) of one of
$Q,R,T_{1},\dots,T_{s-1},$ or $T_{s}$, Therefore by Remark 2.10, we obtain
that the geodesic segments $p_{i}$ are $M$–Morse. For any given
$B^{\prime}\geq B$ (recall that $B$ is the constant of Lemma 2.7 applied with
$M$ and $\rho$) we deduce the following. Since
$\liminf_{i\to\infty}\left\lvert r^{(i)}\right\rvert_{S}=\infty$ and
$\liminf_{i\to\infty}\left\lvert t_{j}^{(i)}\right\rvert_{S}=\infty$ we have
that $\left\lvert r^{(i)}\right\rvert_{S}>B^{\prime}$ and $\left\lvert
t_{j}^{(i)}\right\rvert_{S}>B^{\prime}$ for each $j=1,\dots,s$ and
sufficiently large $i$. Moreover, since
$x_{j}^{(i)}\in\left(Q^{\prime}_{i}\cup R^{\prime}_{i}\right)\setminus
I\subseteq N_{i}^{\prime}\setminus I$ and since $\bigcap N_{i}^{\prime}=I$,
for $i$ large enough we have
$\left\lvert{x_{j}^{(i)}}\right\rvert_{S}>B^{\prime}$. Thus Lemma 2.7 implies
that $p$ is $(N;\lambda,c)$–Morse quasigeodesic. Finally, choosing
$B^{\prime}$ sufficiently large with respect to $\lambda,c,$ and $\left\lvert
g\right\rvert_{S},$ gives us that the endpoints of $p$ are a greater distance
than $\left\lvert g\right\rvert_{S}$ apart, the desired contradiction.

From the above, there is some $i\in\mathbb{N}$ such that $K_{i}$ is separable,
contains the product $QRT_{1}\dots T_{s}$, and excludes $g$. Therefore the
product $QRT_{1}\dots T_{s}=H_{1}\dots H_{n}$ is separable. ∎

###### Proof of Corollary 1.3.

Stable subgroups are quasiconvex, and quasiconvex subgroups of virtually
special groups are separable by [HW08, Corollary 7.9]. Moreover,
$\operatorname{CAT}(0)$ groups have the MLTG property by [RST22, Theorem D].
Therefore Theorem 1.1 applies to give the result. ∎

###### Proof of Corollary 1.4.

Let $H_{1},\dots,H_{n}$ be strongly quasiconvex subgroups of $A_{\Gamma}$. If
there is some $1\leq i\leq n$ such that $H_{i}$ has finite index, then the
product $H_{1}\dots H_{n}$ is a union of finitely many cosets of $H_{i}$.
Since $H_{i}$ has finite index, it is separable in $A_{\Gamma}$, whence
$H_{1}\dots H_{n}$ is separable.

Now suppose that each of the subgroups $H_{1},\dots,H_{n}$ has infinite index
in $A_{\Gamma}$. By [RST23, Corollary 7.4], they must be stable subgroups.
Noting that $A_{\Gamma}$ is $\operatorname{CAT}(0)$ and hence has the MLTG
property [RST22, Theorem D], the conclusion follows by applying Theorem 1.1. ∎

## References

  * [AGM09] Ian Agol, Daniel Groves, and Jason Fox Manning.  Residual finiteness, QCERF and fillings of hyperbolic groups.  Geom. Topol., 13(2):1043–1073, 2009.
  * [Bon96] Mario Bonk.  Quasi-geodesic segments and Gromov hyperbolic spaces.  Geom. Dedicata, 62(3):281–298, 1996.
  * [CH17] Matthew Cordes and David Hume.  Stability and the Morse boundary.  J. Lond. Math. Soc. (2), 95(3):963–988, 2017.
  * [CRSZ22] Matthew Cordes, Jacob Russell, Davide Spriano, and Abdul Zalloum.  Regularity of Morse geodesics and growth of stable subgroups.  J. Topol., 15(3):1217–1247, 2022.
  * [DT15] Matthew Gentry Durham and Samuel J. Taylor.  Convex cocompactness and stability in mapping class groups.  Algebr. Geom. Topol., 15(5):2839–2859, 2015.
  * [Gen19] Anthony Genevois.  Hyperbolicities in ${\rm CAT}(0)$ cube complexes.  Enseign. Math., 65(1-2):33–100, 2019.
  * [GKLS21] Marius Graeber, Annette Karrer, Nir Lazarovich, and Emily Stark.  Surprising circles in Morse boundaries of right-angled Coxeter groups.  Topology Appl., 294:Paper No. 107645, 3, 2021.
  * [HW08] F. Haglund and D.T. Wise.  Special cube complexes.  Geom. Funct. Anal., 17(5):1551–1620, 2008.
  * [IM21] Merlin Incerti-Medici.  Comparing topologies on the Morse boundary and quasi-isometry invariance.  Geom. Dedicata, 212:153–176, 2021.
  * [KMT17] Thomas Koberda, Johanna Mangahas, and Samuel J. Taylor.  The geometry of purely loxodromic subgroups of right-angled Artin groups.  Trans. Amer. Math. Soc., 369(11):8179–8208, 2017.
  * [Liu21] Qing Liu.  Dynamics on the Morse boundary.  Geom. Dedicata, 214:1–20, 2021.
  * [Min06] A. Minasyan.  Separable subsets of GFERF negatively curved groups.  J. Algebra, 304(2):1090–1100, 2006.
  * [MM22] A. Minasyan and L. Mineh.  Quasiconvexity of virtual joins and separability of products in relatively hyperbolic groups.  arXiv preprint, 2022.
  * [QR22] Yulan Qing and Kasra Rafi.  Sublinearly Morse boundary I: ${\rm CAT}(0)$ spaces.  Adv. Math., 404(part B):Paper No. 108442, 51, 2022.
  * [RST22] Jacob Russell, Davide Spriano, and Hung Cong Tran.  The local-to-global property for Morse quasi-geodesics.  Math. Z., 300(2):1557–1602, 2022.
  * [RST23] Jacob Russell, Davide Spriano, and Hung Cong Tran.  Convexity in hierarchically hyperbolic spaces.  Algebr. Geom. Topol., 23(3):1167–1248, 2023.
  * [RZ93] L. Ribes and P.A. Zalesskii.  On the profinite topology on a free group.  Bull. London Math. Soc., 25(1):37–43, 1993.
  * [She23] Sam Shepherd.  Imitator homomorphisms for special cube complexes.  Trans. Amer. Math. Soc., 376(1):599–641, 2023.
  * [Tra19] Hung Cong Tran.  On strongly quasiconvex subgroups.  Geom. Topol., 23(3):1173–1235, 2019.


# Feature matching for multi-epoch historical aerial images

Lulin Zhang1, Ewelina Rupnik1, Marc Pierrot-Deseilligny1 1LASTIG, Université
Gustave Eiffel, ENSG, IGN, F-94160 Saint-Mandé, France - (Lulin.Zhang,
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Historical imagery is characterized by high spatial resolution and
stereoscopic acquisitions, providing a valuable resource for recovering 3D
land-cover information. Accurate geo-referencing of diachronic historical
images by means of self-calibration remains a bottleneck because of the
difficulty to find sufficient amount of feature correspondences under evolving
landscapes. In this research, we present a fully automatic approach to
detecting feature correspondences between historical images taken at different
times (i.e., inter-epoch), without auxiliary data required. Based on relative
orientations computed within the same epoch (i.e., intra-epoch), we obtain
DSMs (Digital Surface Model) and incorporate them in a rough-to-precise
matching. The method consists of: (1) an inter-epoch DSMs matching to roughly
co-register the orientations and DSMs (i.e, the 3D Helmert transformation),
followed by (2) a precise inter-epoch feature matching using the original RGB
images. The innate ambiguity of the latter is largely alleviated by narrowing
down the search space using the co-registered data. With the inter-epoch
features, we refine the image orientations and quantitatively evaluate the
results (1) with DoD (Difference of DSMs), (2) with ground check points, and
(3) by quantifying ground displacement due to an earthquake. We demonstrate
that our method: (1) can automatically georeference diachronic historical
images; (2) can effectively mitigate systematic errors induced by poorly
estimated camera parameters; (3) is robust to drastic scene changes. Compared
to the state-of-the-art, our method improves the image georeferencing accuracy
by a factor of 2. The proposed methods are implemented in MicMac, a free,
open-source photogrammetric software.

###### keywords:

Feature matching, Historical images, Multi-epoch, Pose estimation, Self-
calibration

††journal: Nuclear Physics B

##  1 Introduction

Historical imagery chronicles worldwide land-cover information, and as a
result enables long-term environmental monitoring as well as 3D dynamic change
detection. The images are of high spatial resolution, and are acquired in
stereoscopic configuration. They have been acquired in many countries all over
the world and can be traced back to the beginning of the 20th century [1].
Millions of historical images have been digitized and made accessible through
web services [2, 3, 4]. They are often accompanied by metadata, in most cases
including the camera focal length and the physical sensor size. Other metadata
such as flight plans, camera calibration certificates or orientations are not
commonly available. Besides, inappropriate film/glass plate preservation and
the scanning process enforce reestimating of the camera calibrations (i.e.,
the self-calibrating).
Self-calibration estimates simultaneously interior and exterior orientations.
It is generally solved with a bundle block adjustment (BBA) routine, taking
feature correspondences and Ground Control Points (GCPs) as input
observations. Extracting feature correspondences within a single epoch (also
referred to as intra-epoch) can be efficiently done with local features such
as SIFT [5]. Yet, due to drastic scene changes and heterogeneous acquisition
conditions, it is challenging to automatically find feature correspondences
across different epochs (also referred to as inter-epoch).
In this work we propose a fully automatic approach to computing dense and
robust inter-epoch feature correspondences. Our method consists of a rough co-
registration by finding feature correspondences between DSMs derived within
single epochs, and a precise feature matching on original RGB images. Our main
contributions include:

  * 1.

By matching DSMs, we are able to obtain robust rough co-registration as the 3D
landscape often stays globally stable over time.

  * 2.

Under the guidance of co-registered orientations and DSMs, we reduce the
difficulty in precise matching by: (1) narrowing down the search space; (2)
reducing the combinatorial complexity since only overlapping images are
analyzed;

  * 3.

By proposing a tiling scheme (including one-to-many tiling in rough co-
registration and one-to-one tiling in precise matching), we are opening up the
possibility to scale-up the deep learning methods for feature matching. As we
have shown in Section 4.3, using them out-of-the-box is inefficient. Large
images demand high computing resources, while deep learning feature extraction
methods are presumably trained on small images.

  * 4.

By including priors about the geometry (in form of DSMs), we can filter
candidate correspondences: each three correspondences projected to DSM serve
to compute a 3D Helmert transformation between epochs, and most importantly
provide a 2D constraint on all images’ correspondences.

In the following we briefly describe the current states of local features,
inter-epoch historical images alignment as well as robust matching. In Section
3 we introduce our methodology, and in Section 4 the experiments as well as
results are given. The source code is available from MicMac Github [6].

##  2 Related Work

###  2.1 Local features

Among hand-crafted features, SIFT [5] is undoubtedly a milestone. For decades,
SIFT and its variants such as RootSIFT [7], RootSIFT-PCA [8], and DSP-SIFT [9]
have been widely used both in academic and industrial fields. Other popular
traditional features include ASIFT [10], SURF [11], Hessian Affine keypoint
detector [12] and KAZE [13]. They are respectively powerful in affine
invariance, computational efficiency, scale and affine invariance, retaining
object boundaries.
With the rise of deep learning, learned features have gained much attention
due to their good performance on certain benchmarks. LIFT [14] is the first
end-to-end network that implemented a full pipeline including detection,
orientation estimation and feature description. From then on, a large number
of learned methods related to feature matching proliferated, including:

  1. 1.

Pure descriptor such as L2-Net [15], HardNet [16] and Contextdesc [17].

  2. 2.

Detector along with descriptor. Most of the methods fall into this category,
including LF-Net [18], DELF [19], SuperPoint [20], D2-Net [21], ASLFeat [22],
R2D2 [23], D2D [24], etc.

  3. 3.

Matching methods based on off-the-shelf features, such as SuperGlue [25].
Originally combined with SuperPoint, SuperGlue was designed for real-time
processing, and unlike the previous methods, it is capable of implicitly
encoding spatial relation between neighboring features within an image and
across images with the so-called attention mechanism.

Early learned methods (LIFT [14], L2-Net [15], HardNet [16], DELF [19],
SuperPoint [20], LF-Net [18]) use only intermediate metrics (e.g.,
repeatability, matching score, mean matching accuracy, etc.) to evaluate the
matching performance. However, good performance on benchmarks does not
necessarily imply a better matching quality [26]. Jin et al. [27] introduced a
comprehensive benchmark for local features and robust estimation algorithms,
focusing on the accuracy of the reconstructed camera pose as the primary
metric. Using the new metric, SIFT [5] and SuperGlue [25] take the lead [28].

###  2.2 Inter-epoch historical images alignment

When it comes to inter-epoch historical images, however, directly applying
SIFT or SuperGlue often results in inferior results due to large radiometric
differences. In Figure 8 we showed an example where SIFT and SuperGlue failed
on an inter-epoch image pair with drastic scene changes. It is understandable
as (1) SIFT is not sufficiently invariant over time, while (2) SuperGlue is
not invariant to rotations and it underperforms on larger images because it
was presumably trained on small images.
Therefore, many previous researches bypassed the task of extracting inter-
epoch correspondences by processing different epochs separately followed by an
inter-epoch co-registration relying on Ground Control Points(GCPs). Between 10
and 169 GCPs are required in [29], [30], [31], [32], [33]. GCPs are usually
measured with the help of photointerpretation on recent orthophotos, however,
it is still monotonous and time-consuming. Furthermore, it is difficult to
find salient points that are stable over time.
Certain attempts were made to extract inter-epoch correspondences. Giordano et
al. [34] extract feature correspondences between historical and recent images
relying on HoG descriptors [35]. The authors require flight plans as input,
which are not commonly available as mentioned in Section 1. Feurer et al.
[36], Filhol et al. [37], Cook et al. [38], Parente et al. [39] and Blanch et
al. [40] assume that a sufficient number of keypoints remain invariant across
time and employ SIFT to extract inter-epoch feature correspondences. It
remains questionable whether the method is capable of handling drastic scene
changes. Zhang et al. [41] extract inter-epoch correspondences from SIFT-
detected keypoints based on the hypothesis that points follow 2D and 3D
spatial similarity model. This method works in simple cases with few scene
changes. Additionally, a stream of research works focuses on historical
terrestrial images ( [42], [43], [44], [45]) and historical video recordings (
[45]). However, their algorithms are not suitable to the aerial case.
This work is an extention of [41]. Unlike in [41], we introduce a rough co-
registration between different epochs based on matching DSMs with SuperGlue,
and use it to guide a precise matching. Our rough co-registration is robust
under extreme scene changes because (1) SuperGlue utilizes context to enhance
feature descriptors and (2) DSMs are generally stable over time. With the
guidance of roughly co-registered orientations and DSMs, both SIFT and
SuperGlue achieved good performance, as shown in our experiments.

###  2.3 Robust matching

The goal of robust matching is to tell apart inliers (i.e., true
correspondences) from outliers (i.e., false correspondences), and eliminate
the latter from further processing. Typically, an iterative sampling strategy
based on RANSAC (Random Sample Consensus) [46] relying on some geometric
model, such as homography [47] or essential matrix [47] is carried out to
remove outliers. Substitutes of RANSAC include hand-crafted methods such as
LMedS [48] (a meaningful groundwork before RANSAC), MLESAC [49] (maximizing
inlier likelihood), PROSAC [50] (choosing samples progressively for
acceleration), DEGENSAC [51] (increasing robustness when a dominant plane is
present), GC-RANSAC [52] (adopting graph-cut algorithm), MAGSAC [53] (defining
threshold automatically with marginalization), and learned methods such as
DSAC [54] (simulating RANSAC in a differentiable way) and CNe [55] (end-to-end
deep network).
In this research, we use RANSAC to estimate the 3D Helmert transformation
between surfaces (i.e., DSMs) calculated in different epochs. Compared to the
classical essential/fundamental matrix filtering, with less data (3 versus 5
points) we impose stricter rules on the sets of points. Lastly, we eliminate
the remaining false correspondences by looking at their cross-correlation.

##  3 Methodology

(a) Full processing workflow

(b) Workflow of the rough co-registration.

(c) Four rotation hypotheses

Figure 1: (a) Our full processing workflow. Intra-epoch feature
correspondences, relative orientations and DSM are obtained within each epoch
individually. Once multi-epochs are roughly co-registered based on feature
correspondences computed on DSMs, a precise matching is carried out by
performing patch matching or guided matching, followed by a 3D-RANSAC filter
and a local-context validation, giving rise to the final inter-epoch feature
correspondences. (b) Workflow of the rough co-registration. SuperGlue is
applied on tiles of DSMs, followed by a RANSAC procedure to remove outliers.
As SuperGlue is not invariant to rotations lager than 45∘, we test four
rotation hypotheses and keep the best one. (c) Four rotation hypotheses. We
rotate the secondary DSM by 90∘ four times to match with master DSM and keep
the best one with the largest number of RANSAC inliers.

Figure 1(a) exhibits the workflow of our pipeline. It includes 3 main parts:
intra-epoch processing, inter-epoch processing, and combined processing. The
images’ focal lengths and physical sensor sizes are supposed to be known, as
we already mentioned in Section 1. The images are resampled to the geometry of
the fiducial marks prior to processing. For the sake of simplicity, we only
exhibit the processing flow of two epochs, however, it can be easily extended
to more epochs.
We use the term Patch and Guided to refer to two alternatives of our rough-to-
precise matching pipeline (patch matching or guided matching). The patch
matching uses learned features, and the guided matching uses hand-crafted
features. We explore the two approaches because of their diverse
characteristics: the Patch recovers more correspondences (cf. Figure 8 (b) and
(c)), however, due to the resampling stage, they are computationally more
intensive compared to Guided. Additionally, our collateral goal is to compare
their performances.
We adopt the following naming conventions – (1) images acquired in epoch1 are
denoted as $I^{e_{1}}$, and images in epoch2 as $I^{e_{2}}$; (2) orientations
of epoch1 and epoch2 are denoted as $O^{e_{1}}$ and $O^{e_{2}}$; (3) DSMs of
epoch1 and epoch2 are denoted as $D^{e_{1}}$ and $D^{e_{2}}$.

###  3.1 Intra-epoch Processing

Intra-epoch processing furnishes complementary 3D information to help the
latter inter-epoch processing. It is a standard photogrammetry or SfM pipeline
and can be accomplished with lots of solutions (e.g. MicMac [56], COLMAP [57],
OpenMVG [58], Theia [59], etc.). The solution used in our experiment is
MicMac. It is performed within each epochi individually as follows:

  1. 1.

Extract intra-epoch correspondences between images $I^{e_{i}}$ with SIFT [5];

  2. 2.

Compute interior and relative orientations ($O_{ini}^{e_{i}}$) with the
sequential SfM;

  3. 3.

Based on image orientations $O_{ini}^{e_{i}}$, perform semi-global dense
matching [60] between images $I^{e_{i}}$ to get DSM ($D_{ini}^{e_{i}}$) in
their arbitrary coordinate frames.

###  3.2 Inter-epoch Processing

Inter-epoch processing follows a rough-to-precise matching strategy including:

  1. 1.

Rough co-registration: match DSMs based on one-to-many tiling scheme
(described in Section 3.2.1) to roughly co-register the image orientations and
DSMs into a common reference frame.

  2. 2.

Precise matching: match original RGB images under the guidance of the roughly
co-registered data.

We choose matching DSMs for rough co-registration and RGB images for precise
matching because: (1) DSMs computed from historical images turn to be more
constant over time compared to RGB images; however, (2) due to low radiometric
quality of the images, the DSMs are noisy (cf. Figure 7(d)). As a consequence,
DSMs provide correspondences perfectly tailored for rough co-registration -
less accurate but repetitive over time. Then, by narrowing down the search
space based on co-registered orientations, matching RGB images will give us
correspondences that are both accurate and robust.

####  3.2.1 Rough co-registration

The workflow of rough co-registration is displayed in Figure 1(b). The results
(co-registered orientations and DSMs) will be used to (1) guide the precise
matching, and (2) provide initial orientations for bundle adjustment in the
combined processing. Therefore, the goal is to get a moderate number of
reliable inter-epoch feature correspondences at low cost. We choose SuperGlue
111We use the pre-trained SuperGlue model provided by the authors. Keypoint
locations as well as confidence scores are predicted with SuperPoint [20]. to
do the rough matching as it is more invariant over time than SIFT. Simply
applying SuperGlue on multi-epoch images fails because it is not rotation
invariant and it underperforms on large images. We improve the matching
robustness by adopting the following modifications:

##### (1) Match DSMs instead of original images

Assuming the numbers of the original images in epoch1 and epoch2 are P and Q
respectively, instead of applying brute force rough matching P$\times$Q times
on original images, we only need to do it once on a pair of DSM images. To
obtain the DSM images, we convert the DSMs (i.e., 2.5D rasters) from floating-
point to [0-255] range grayscale images. Since the DSMs contain outliers,
pixels with elevations larger than double the standard deviation of all
elevations are ignored in the conversion process. Matching DSM images has the
following merits:

  1. 1.

Redundancy caused by the forward and side overlapping areas is removed;

  2. 2.

It implicitly enables a follow-up search for globally consistent inliers;

  3. 3.

It decreases the combinatorial complexity caused by rotation ambiguity of the
product of P and Q images;

  4. 4.

Even under important scene changes, DSMs generally provide stable information
over time.

##### (2) One-to-many tiling scheme

As SuperGlue provides unsatisfactory result on large images, we propose a one-
to-many tiling scheme to make up for the deficiency. It is performed as
follows (cf. Figure 1(b)):

  1. 1.

Crop both the master and secondary DSM images into M and N tiles of certain
size, respectively;

  2. 2.

Apply SuperGlue on M$\times$N tile pairs individually;

  3. 3.

Merge the feature correspondences and perform RANSAC based on 2D similarity
transformation to find globally consistent inliers.

##### (3) Testing 4 rotation hypotheses

As SuperGlue is not invariant to rotations larger than 45∘, the one-to-many
tiling scheme will be performed repeatedly over 4 rotation hypotheses as
follows (cf. Figure 1(c)):

  1. 1.

Rotate the secondary DSM image by 90∘ four times;

  2. 2.

Match each rotated DSM image with the master DSM image with the one-to-many
tiling scheme;

  3. 3.

Keep the hypothesis with the largest number of RANSAC inliers.

With RANSAC inliers, we build a 3D Helmert transformation to transform the
orientations $O_{ini}^{e_{1}}$ and DSM $D_{ini}^{e_{1}}$ to the frame of
${epoch_{2}}$, leading to co-registered orientations $O_{co}^{e_{1}}$ and DSM
$D_{co}^{e_{1}}$.
The feature matching method based on DSMs will fail on perfectly flat terrain.
In this rare scenario, feature correspondences can be extracted on orthophotos
(see demonstration in Section 4.3), and their elevation coordinates can be
retrieved from the DSMs. However, because the scene radiometry change is more
pronounced than the 3D landscape change, fewer matches are found.

####  3.2.2 Precise matching

To compute precise inter-epoch feature correspondences, we perform matching on
original RGB images under the guidance of co-registered orientations and DSMs.
It consists of extracting tentative inter-epoch feature correspondences,
followed by a 3D-RANSAC filter and a local-context validation stage to remove
outliers.

##### (1) Get tentative inter-epoch correspondences

We offer two alternatives to get tentative correspondences: patch matching or
guided matching. The former has better overall performance while the latter is
more efficient in terms of the use of memory and CPU resources.
(1-1) Patch matching for learned features It is based on one-to-one tiling
scheme (not to confuse with the one-to-many tiling scheme presented in the
rough matching), as shown in Figure 2(a) and detailed below:

  1. 1.

Crop the master RGB image $I^{e_{1}}$ into M tiles ($T^{e_{1}}$) of certain
size, with a buffer zone overlapped with each other;

  2. 2.

As the orientations and DSMs are roughly co-registered, we project each master
tile $T^{e_{1}}$ onto the DSM $D_{co}^{e_{1}}$ and backproject to secondary
RGB image $I^{e_{2}}$ to find the corresponding tile $T^{e_{2}}$;

  3. 3.

Resample $T^{e_{2}}$ to $\widetilde{T}^{e_{2}}$, so that the tile pair
$P({T^{e_{1}},\widetilde{T}^{e_{2}}})$ is free from differences of rotation,
scale and extent;

  4. 4.

Apply SuperGlue on tile pair $P({T^{e_{1}},\widetilde{T}^{e_{2}}})$ to find
feature correspondences $M({\mathbf{K}^{e_{1}},\mathbf{K}^{e_{2}}})$
($\mathbf{K}^{e_{i}}$ represents keypoints in $I^{e_{i}}$);

  5. 5.

Merge the feature correspondences.

(a) patch matching

(b) Buffer zone

(c) guided matching

Figure 2: (a) and (c) illustrate toy-examples of the patch matching and guided
matching, respectively, (b) displays the feature correspondences where
$\mathbf{K}^{e_{1}}$ exceeds the original tile size (dark green area) and
therefore will be abandoned.

Our patch matching experiments are performed based on SuperGlue, however,
other learned methods can be adopted readily. Because the orientations and
DSMs are only roughly co-registered, we have to take into account the margin
of error when projecting tiles to overlapping images. This is why we add a
buffer zone around the tile $T^{e_{1}}$, as shown in Figure 2(b). To avoid
duplicate matches, we remove a posteriori the matches exceeding the original
tile.
(1-2) Guided matching for hand-crafted features. The patch matching substitute
orientated towards hand-crafted features is the guided matching, as shown in
Figure 2(c). It leverages the positions of predicted keypoints, the known
scale ratio and rotation differences to narrow down the list of the matching
candidates. In our experiments, we use the SIFT points, but the method is
suitable to any hand-crafted extractor. The strategy is as follows:

  1. 1.

Compute the scale ratio $R_{scl}$ and the rotation $D_{rot}$ between two
images by sequentially projecting the $I^{e_{1}}$ image corners to the co-
registered DSM $D_{co}^{e_{1}}$ and to image $I^{e_{2}}$;

  2. 2.

Extract keypoints $\mathbf{K}^{e_{1}}$ in image $I^{e_{1}}$ and
$\mathbf{K}^{e_{2}}$ in image $I^{e_{2}}$;

  3. 3.

Intersect the keypoints $\mathbf{K}^{e_{1}}$ with the co-registered DSM
$D_{co}^{e_{1}}$;

  4. 4.

Back-project them to image $I^{e_{2}}$, giving rise to predicted keypoints
$\widetilde{\mathbf{K}}^{e_{2}}$;

  5. 5.

Search for a subset of points in $\mathbf{K}^{e_{2}}$ located within a radius
$S$ (100 pixels in our experiments) centered at the predicted positions
$\widetilde{\mathbf{K}}^{e_{2}}$;

  6. 6.

Remove candidate matches whose scales and rotations computed by SIFT are
incoherent with $R_{scl}$ and $D_{rot}$ computed from image orientations and
the co-registered DSM (i.e., step 1);

  7. 7.

Find the best match with mutual nearest neighbor and apply the first to second
nearest neighbor ratio test [5].

##### (2) Get enhanced inter-epoch correspondences

To compute enhanced correspondences, we apply a 3D-RANSAC filter on the
previously obtained tentative matches. More precisely, we do the following:
(1) for each feature correspondence
$M({\mathbf{K}^{e_{1}},\mathbf{K}^{e_{2}}})$, the keypoints
$\mathbf{K}^{e_{1}}$ and $\mathbf{K}^{e_{2}}$ are projected onto DSM
$D_{co}^{e_{1}}$ and $D_{ini}^{e_{2}}$ individually to get 3D points
$M({\mathbf{G}^{e_{1}},\mathbf{G}^{e_{2}}})$; and (2) the correspondences
$M({\mathbf{G}^{e_{1}},\mathbf{G}^{e_{2}}})$ are iteratively sampled to
compute the 3D spatial similarity RANSAC model:

${G_{i}^{e_{2}}}=\lambda\cdot\mathbf{R}\cdot{G_{i}^{e_{1}}}+\mathbf{T},\quad
i\in[1,3]$

where $\lambda$ is the scale factor, $\mathbf{T}$ is the translation vector
and $\mathbf{R}$ is the rotation matrix. We set the number of RANSAC
iterations to 1000, and consider feature correspondences within $T_{r}$ of its
predicted position as inliers. In our experiment, $T_{r}$ was set to
10$\times$$GSD$ where $GSD$ is the mean ground sampling distance in the
coordinate frame of epoch ${e_{2}}$. This distance is computed as the ground
distance between two adjacent image pixels.

##### (3) Get final inter-epoch correspondences

In the preceding step we got rid of a substantial number of outliers, however,
we believe that not all outliers could be identified. We apply cross-
correlation for final validation. Feature correspondences with their
correlation scores below a predefined threshold (0.6 in our experiments) are
discarded. The correlation window size was set to be large enough to take into
account the context around a point (32$\times$32 pixels in our experiment).
Figure 3 shows an example of a false match (red) eliminated by cross
correlation, while the true match (blue) is kept.

Figure 3: Demonstration of the validation with cross-correlation. Considering
poor quality of historical images, the window size (blue and red rectangles)
was set to 32$\times$32 pixels. False match (red) is eliminated by cross
correlation, while true match (blue) is kept.

###  3.3 Combined Processing

Based on the intra-epoch and inter-epoch feature correspondences, a free
network BBA is performed to refine all the image orientations and camera
calibrations. If the results need to be analyzed in a metric scale, a spatial
similarity transformation will be performed to move the refined acquisitions
in an arbitrary reference frame to a metric one. If the precise orientations
for one of the epochs were known (i.e., deemed as ground truth), their
parameters will be fixed during the BBA and the subsequent spatial similarity
transformation will be skipped. We adopted the Fraser model [61] to calibrate
the cameras and allowed image-dependent affine parameters, the remaining
parameters were shared among all images.

##  4 Experiment

###  4.1 Implementation details

To reduce the image noise all historical images entering the pipeline are
downsampled by a factor of 3. To calculate the DSMs, we further downsample the
images by a factor of 4, which amounts to a total downsampling factor of 12
with respect to the input images. For example, the images in Fréjus 1970 are
downsampled from [8766, 8763] to [730, 730]. Note that the DSMs serve 3
purposes: (1) rough co-registration, (2) narrowing down the search space in
precise matching, (3) providing 3D coordinates for 3D-RANSAC filter. A low
resolution surface is good enough for these tasks, and keeps the computational
cost low.
To balance the number of the intra- and inter-epoch feature correspondences,
we perform intra-epoch feature correspondences reduction available in MicMac
[62], followed by setting the relative observation weight in the bundle
adjustment, if necessary. The feature correspondences reduction algorithm
maximizes good spatial distribution, points’ multiplicity and low reprojection
error, it also helps to speed up the bundle adjustment.
Inter-epoch feature correspondences are extracted for every possible
combination of 2 epochs and finally merged. Algorithm 1 demonstrates the
pseudo code of our pipeline.

Input: Images; focal lengths; physical sensor sizes;

Output: Interior and relative orientations;

1 for  _each epoch i_ do

2 Downsample historical images by a factor of 3;

3 Extract intra-epoch feature correspondences with SIFT, followed by a
reduction of feature correspondences, giving rise to $T_{intra}$;

4 Recover interior and relative orientations $O_{ini}^{e_{i}}$;

5 Calculate DSMs $D_{ini}^{e_{i}}$ based on images further downsampled by a
factor of 4;

6

7 end for

8Set the most recent epoch as reference epoch $E_{r}$, the other epochs as
free epoch $E_{f}$;

9 for  _every epoch $E_{f}$_ do

10 Convert DSM $D_{ini}^{e_{f}}$ and $D_{ini}^{e_{r}}$ to grayscale images,
leading to $I_{DSM}^{e_{f}}$ and $I_{DSM}^{e_{r}}$;

11 Match $I_{DSM}^{e_{f}}$ and $I_{DSM}^{e_{r}}$ with one-to-many tiling
scheme;

12 Transform orientations $O_{ini}^{e_{f}}$ and DSM $D_{ini}^{e_{f}}$ to the
frame of epoch $E_{r}$, leading to orientations $O_{co}^{e_{f}}$ and DSM
$D_{co}^{e_{f}}$;

13

14 end for

15for  _each inter-epoch image pair_ do

16 Get tentative feature correspondences with patch or guided matching, based
on orientations $O_{co}$ and DSM $D_{co}$;

17 Get enhanced feature correspondences with 3D-RANSAC filter, based on
orientations $O_{co}$ and DSM $D_{co}$;

18 Get final feature correspondences (i.e., $T_{inter}$) with cross
correlation;

19

20 end for

21Refine orientations $O_{co}$ with $T_{intra}$ and $T_{inter}$;

Algorithm 1 Our pipeline

###  4.2 Datasets

We tested our method on three datasets: Pezenas, Fréjus and Kobe. Details of
the datasets are listed in Table 1 and Figure 4.

##### Pezenas

It is a 420 $km^{2}$ rectangular area located in Pezenas in the Occitanie
region in southern France. The area is mainly covered with vegetation and
several sparsely populated urban zones. We have at our disposal three sets of
images acquired in 1971, 1981 and 2015. The epoch 2015 was acquired with the
IGN’s digital metric camera [63], and is treated as the ground truth (GT)
during our processing (in other words, the 2015 image orientations are fixed
during BBA). The area exhibits changes in scene appearance in the 44-year
period.

##### Fréjus

It is a 15 $km^{2}$ rectangular area located in Fréjus, a commune in
southeastern France. The area is mainly covered with buildings along with
scattered farmlands, except a half-moon-shaped bay located in south. We have
four sets of images acquired in 1954, 1966, 1970 and 2014. The epoch 2014 is
treated as GT. The area exhibits drastic scene changes in the 60-year period,
as can be seen in Figure 5, where the evolution of a subregion is displayed.

##### Kobe

It is a 90 $km^{2}$ area of irregular shape located in the north of Awaji
Island, Japan. The well-known Kobe earthquake happened here in January 1995.
We have two sets of images: pre-event acquired in 1991 and post-event acquired
in 1995. It is mainly covered with mountain area and narrow urban zones along
the sea. There is no GT, hence we measured 2 points on Google map to roughly
scale the result to metric units. In this dataset we are interested in
localizing the earthquake fault.

| Pezenas | Fréjus | Kobe
---|---|---|---
| E1971 | E1981 | E2015 | E1954 | E1966 | E1970 | E2014 | E1991 | E1995
F [pix] | 7589 | 7607 | 9967.5 | 9204.5 | 23350 | 10230 | 10230 | 18281 | 7662 | 7662
Wid [mm] | 230 | 230 | 47 | 50 | 300 | 180 | 180 | 99.28 | 212 | 212
Hei [mm] | 230 | 230 | 35 | 36 | 300 | 180 | 180 | 72.42 | 212 | 212
GSD [m] | 0.32 | 0.59 | 0.46 | 0.5 | 0.11 | 0.17 | 0.17 | 0.35 | 0.5 | 0.18
F. o. | 60% | 60% | 60% | 60% | 60% | 60% | 60% | 60% | 65% | 65%
S. o. | 20% | 20% | 50% | 50% | 20% | 30% | 30% | 30% | 35% | 65%
H [m] | 2400 | 4500 | 4600 | 4600 | 2500 | 1700 | 1700 | 6500 | 3800 | 1400
Nb | 57 | 27 | 308 | 74 | 19 | 15 | 19 | 33 | 15 | 83
Table 1: Dataset details of Pezenas, Fréjus and Kobe. The 2015 acquisition of
Pezenas is obtained with two sets of camera. E stands for epoch, F means focal
length, Wid and Hei are the width and height of image, GSD is the ground
sampling distance, F.o. and S.o. are forward and side overlap, H is the flying
height, Nb is the number of images.

(a) Pezenas

(b) Fréjus

(c) Kobe

Figure 4: Orthophotos of the datasets used in experiments.

(a) Fréjus 1954

(b) Fréjus 2014

Figure 5: Evolution of a subregion in Fréjus.

###  4.3 Comparison of co-registration with DSMs and orthophotos

In this section, we carry out the rough co-registration twice, with DSMs and
orthophotos. We choose Fréjus in the years 1954 and 2014 for this experiment
because it is the most difficult matching scenario, demonstrating significant
scene changes and very different area extents. As can be seen in Figure
6(a-b), the whole of the epoch 1954 overlaps by solely $\approx$10 % with the
epoch 2014. The sizes of orthophotos for Fréjus 1954 and 2014 individually are
519$\times$363 pixels and 1888$\times$1361 pixels. The DSM has the same size
as orthophoto. For both types of images (i.e., DSM rasters and orthophotos) we
apply Ours (one-to-many tiling scheme with SuperGlue, followed by RANSAC),
pure SuperGlue and SIFT. Feature correspondences are visualised in Figure 6
(c-h) and the inlier ratios are given in Table 2.
As can be seen, Ours inlier ratio reached 18.6% and 37.1% for orthophotos and
DSMs, respectively. SuperGlue and SIFT reached 0% for both types of images.
DSMs tend to have more inliers than orthophotos, because there is much less
perceived change in the DSMs than there is in the orthophotos, see Figure 6
(c) and (d). Our rough co-registration is capable of recovering robust
correspondences because one-to-many tiling scheme finds many good
correspondences in the overlapping tile pairs, and few wrong ones in the non-
overlapping pairs, which ensures the success of RANSAC.

(a) Overlapping areaOrthophoto

(b) Overlapping areaDSM

(c) OursOrthophoto

(d) OursDSM

(e) SuperGlueOrthophoto

(f) SuperGlueDSM

(g) SIFTOrthophoto

(h) SIFTDSM

Figure 6: Comparison of 3 rough co-registration methods (Ours, SuperGlue and SIFT) on 2 types of images (orthophotos and DSMs). Each subfigure displays Fréjus 1954 on the left and 2014 on the right. (a, b) The overlapping areas (red rectangles) of orthophotos and DSMs individually. (c, e, g) Matching result of Ours, SuperGlue and SIFT on orthophotos. (d, f, h) Matching result of Ours, SuperGlue and SIFT on DSMs. |  | Correspondences |  |
---|---|---|---|---
|  | number | Inlier number | Inlier ratio
orthophoto | SIFT | 2 | 0 | 0%
SuperGlue | 29 | 0 | 0%
Ours | 290 | 54 | 18.6%
DSM | SIFT | 2 | 0 | 0%
SuperGlue | 24 | 0 | 0%
Ours | 302 | 112 | 37.1%
Table 2: Feature correspondences number, inlier number and inlier ratio in
Figure 6. 3 rough co-registration methods (Ours, SuperGlue and SIFT) on 2
types of images (orthophoto and DSM) are displayed.

###  4.4 Comparison of precise matching on DSMs and original images

In order to decide which type of images (DSMs or original images) is more
suitable for executing the precise matching, we apply our pipeline Patch on
both DSMs and original images of Fréjus 1970 and 2014. The final feature
correspondences are displayed in Figure 7 (a) and (b). To asses quantitatively
the results, we created a GT depth map and calculated the accuracy (correct
matches / total matches). In Figure 7 (c) we plot the accuracy curves while
varying the reprojection error threshold from 0 to 10 pixels. It is clear that
the result using the original images is more accurate, even though the DSMs
recovered more correspondences. This is because historical DSMs at full
resolution are too noisy to guarantee high precision measurements (see the DSM
shaded image in Figure 7 (d).

(a) Feature correspondences on RGB images

(b) Feature correspondences on DSMs

(c) Accuracy of (a) and (b)

(d) shaded image of historical DSM

Figure 7: Comparison of precise matching on original RGB images and DSMs.

###  4.5 Comparison between Ours (Guided and Patch), SuperGlue and SIFT

In this section we compare the feature matching result between our rough-to-
precise pipeline (Guided and Patch), SuperGlue and SIFT. The image pair used
here is Fréjus 1954 and 1966. Given that SuperGlue failed on the original pair
(cf. Figure 8 (a)), the left image was rotated to match the upright position
of the right image prior to processing. SuperGlue is employed with the off-
the-shelf trained outdoor model provided by the authors. We can observe that
pure SIFT fails to extract correct matches, and the out-of-the-box SuperGlue
finds a lot of matches, most of which seem good, but at a closer look the
details reveal poor localization precision. In Guided and Patch, all matches
are correct, and the Patch approach detects denser correspondences, as can be
seen in Figure 8 (b-e).
In Figure 8 (f) we plot the accuracy curves as we vary the reprojection error
threshold from 0 to 10 pixels based on manually created GT. Notice that SIFT
and SuperGlue had 0% and 30% correct matches under 10 pixels, respectively. In
the meantime, both Guided and Patch reached over 90% under 6 pixels.

(a) SuperGlue on images with rotation

(b) Guided

(c) Patch

(d) SIFT

(e) SuperGlue

(f) Accuracy of (b-e)

Figure 8: Comparison between Ours (Guided and Patch), SuperGlue and SIFT. (a)
Out-of-the-box SuperGlue applied on an image pair with a rotation by
$180^{\circ}$. (b-e) Inter-epoch feature correspondences of our method (Guided
and Patch), SIFT and Out-of-the-box SuperGlue applied on the same image pair
without rotation. (f) Accuracy.

###  4.6 Evaluation Method

We refine the co-registered orientations in a BBA routine using the inter- and
intra-epoch feature correspondences (cf. Figure 1(a)). The orientations of the
latest epochs in Pezenas and Fréjus were treated as fixed during the combined
BBA since they were accurately known a-priori, while all the remaining
orientations were considered as free parameters. At first, interior
orientation parameters were shared among all images. Once stable initial
values were known, interior parameters were further refined with image-
dependent affine parameters. The affine component of the camera calibration is
expected to model, at least partially, the shear of the analog film.
For comparison, we recover the image orientations with 5 methods, which can be
divided into 2 categories of ”rough” and ”refined”:

  1. 1.

3vGCPs (”rough”): We co-register all the epochs with 3 manually measured
virtual GCPs in total. By virtual GCPs we mean points of known 3D coordinates
in one of the epochs.

  2. 2.

Co-Reg (”rough”): We match DSMs based on one-to-many tiling scheme to roughly
co-register all the epochs (cf. Figure 1(b)).

  3. 3.

PureSG (”refined”): Based on the Co-Reg result, we remove the scale, rotation
and extent differences for inter-epoch image pairs (without one-to-one tiling
scheme), and apply off-the-shell SuperGlue to get inter-epoch feature
correspondences, which are then used to refine the Co-Reg result.

  4. 4.

Guided (”refined”): Based on the Co-Reg result, inter-epoch feature
correspondences are extracted with guided matching, as presented in Section
3.2.2. Finally, a BA routine refines the initial Co-Reg orientations.

  5. 5.

Patch (”refined”): Same as in Guided, except here we adopt the alternative
Patch approach as presented in Section 3.2.2.

Our first metric to evaluate the aforementioned 5 methods is DoD (Difference
of DSMs). For Pezenas and Fréjus datasets, DoDs are calculated between
historical epochs and the available GT epochs (cf. Figure 9). For Kobe dataset
there is no GT so we calculate the DoD between 1991 and 1995 instead (cf.
Figure 10 (a-e)). Statistics of the presented DoDs are demonstrated in Table
3.
Even though DoD provides a global evaluation over the entire scene, it is
contaminated with real scene changes. Therefore, we present a supplementary
metric for each dataset:

  1. 1.

For Pezenas and Fréjus, to evaluate the orientations resulting from the 5
methods, we manually measured 15 check points that are uniformly distributed
across images. The ground coordinate of the check points are measured in the
GT epoch (with measurement accuracy of 0.17m in Frejus and 0.55m in Pezenas).
Average value $\mu$, standard deviation $\sigma$, and absolute average value
$|\mu|$ of the x-, y- and z-coordinate are displayed in Table 4.

  2. 2.

For Kobe, we: (1) calculated the DSMs; (2) orthorectified the images; and (3)
performed 2D correlation of the respective orthophotos [64] to see whether we
can observe the slip of the tectonic plate. The displacement maps (i.e.,
northeastward, denoted as $Gd$) are presented in Figure 10(f-j).

|  $\mu$ [m] |  $\sigma$ [m] |  $|\mu|$ [m] |  |  $\mu$ [m] |  $\sigma$ [m] |  $|\mu|$ [m]
---|---|---|---|---|---|---|---
DoD${}_{Frejus\\_1954}^{3vGCPs}$ | 12.89 | 12.35 | 14.12 |  DoD${}_{Pezenas\\_1971}^{3vGCPs}$ | -3.89 | 18.07 | 14.06
DoD${}_{Frejus\\_1954}^{Co-Reg}$ | 1.47 | 6.04 | 4.66 |  DoD${}_{Pezenas\\_1971}^{Co-Reg}$ | 7.33 | 17.22 | 16.09
DoD${}_{Frejus\\_1954}^{PureSG}$ | 2.15 | 9.09 | 7.37 |  DoD${}_{Pezenas\\_1971}^{PureSG}$ | -0.03 | 2.92 | 1.94
DoD${}_{Frejus\\_1954}^{Guided}$ | 0.39 | 3.48 | 2.53 |  DoD${}_{Pezenas\\_1971}^{Guided}$ | -0.48 | 2.25 | 1.38
DoD${}_{Frejus\\_1954}^{Patch}$ | -1.69 | 3.44 | 2.51 |  DoD${}_{Pezenas\\_1971}^{Patch}$ | -0.42 | 2.30 | 1.40
DoD${}_{Frejus\\_1966}^{3vGCPs}$ | 7.49 | 6.93 | 8.32 |  DoD${}_{Pezenas\\_1981}^{3vGCPs}$ | -6.80 | 8.89 | 8.21
DoD${}_{Frejus\\_1966}^{Co-Reg}$ | 2.06 | 5.41 | 4.49 |  DoD${}_{Pezenas\\_1981}^{Co-Reg}$ | 0.46 | 9.87 | 7.97
DoD${}_{Frejus\\_1966}^{PureSG}$ | 2.88 | 5.76 | 4.81 |  DoD${}_{Pezenas\\_1981}^{PureSG}$ | -0.52 | 2.45 | 1.61
DoD${}_{Frejus\\_1966}^{Guided}$ | -0.32 | 3.69 | 2.62 |  DoD${}_{Pezenas\\_1981}^{Guided}$ | -0.45 | 2.07 | 1.29
DoD${}_{Frejus\\_1966}^{Patch}$ | -0.40 | 3.62 | 2.43 |  DoD${}_{Pezenas\\_1981}^{Patch}$ | -0.34 | 2.04 | 1.25
DoD${}_{Frejus\\_1970}^{3vGCPs}$ | 3.65 | 7.12 | 6.83 |  DoD${}_{Kobe}^{3vGCPs}$ | 5.58 | 12.92 | 7.18
DoD${}_{Frejus\\_1970}^{Co-Reg}$ | 1.30 | 5.34 | 4.44 |  DoD${}_{Kobe}^{Co-Reg}$ | -11.28 | 17.29 | 15.81
DoD${}_{Frejus\\_1970}^{PureSG}$ | 1.99 | 4.98 | 4.17 |  DoD${}_{Kobe}^{PureSG}$ | 2.48 | 12.11 | 4.34
DoD${}_{Frejus\\_1970}^{Guided}$ | -0.61 | 3.76 | 2.77 |  DoD${}_{Kobe}^{Guided}$ | 1.79 | 12.21 | 4.13
DoD${}_{Frejus\\_1970}^{Patch}$ | -0.74 | 3.36 | 2.17 |  DoD${}_{Kobe}^{Patch}$ | 1.58 | 9.18 | 3.62
Table 3: Average value $\mu$, standard deviation $\sigma$, and absolute average value $|\mu|$ of all the DoDs in Figure 9 and 10(a-e). |  $\mu$ [m] |  $\sigma$ [m] |  $|\mu|$ [m]
---|---|---|---
| x | y | z | x | y | z | x | y | z
DoCheckPt${}_{Frejus}^{3vGCPs}$ | -0.38 | 2.32 | -5.80 | 12.05 | 9.75 | 6.67 | 8.85 | 7.93 | 7.34
DoCheckPt${}_{Frejus}^{Co-Reg}$ | 7.23 | 6.18 | -3.08 | 13.75 | 6.98 | 3.79 | 10.60 | 7.43 | 3.96
DoCheckPt${}_{Frejus}^{PureSG}$ | 1.22 | 0.32 | -1.93 | 2.16 | 1.58 | 4.00 | 1.37 | 1.00 | 3.27
DoCheckPt${}_{Frejus}^{Guided}$ | 1.40 | 0.71 | -0.70 | 2.08 | 1.71 | 3.04 | 1.52 | 1.06 | 2.70
DoCheckPt${}_{Frejus}^{Patch}$ | 1.22 | 0.59 | -0.24 | 2.15 | 1.67 | 1.63 | 1.39 | 0.95 | 1.17
DoCheckPt${}_{Pezenas}^{3vGCPs}$ | 5.17 | -2.85 | 8.26 | 5.24 | 3.67 | 13.78 | 6.01 | 3.47 | 10.44
DoCheckPt${}_{Pezenas}^{Co-Reg}$ | 33.15 | 22.42 | 1.40 | 18.80 | 16.65 | 14.48 | 33.15 | 23.02 | 11.94
DoCheckPt${}_{Pezenas}^{PureSG}$ | -0.44 | -0.56 | 0.55 | 2.73 | 1.60 | 3.02 | 1.95 | 1.44 | 2.15
DoCheckPt${}_{Pezenas}^{Guided}$ | 0.17 | -0.29 | 0.31 | 2.40 | 1.42 | 1.55 | 1.57 | 1.24 | 1.22
DoCheckPt${}_{Pezenas}^{Patch}$ | -0.04 | -0.29 | 0.46 | 2.35 | 1.42 | 1.46 | 1.64 | 1.22 | 1.10
Table 4: Accuracy of 5 different sets of orientation resulting from 5 methods
(3vGCPs, Co-Reg, PureSG, Guided and Patch), evaluated on 15 check points
uniformly distributed in the block for both Fréjus and Pezenas. Average value
$\mu$, standard deviation $\sigma$, and absolute average value $|\mu|$ of the
x-, y- and z-coordinate are displayed for each method.

(a) DoD${}_{Pezenas\\_1971}^{3vGCPs}$

(b) DoD${}_{Pezenas\\_1971}^{Co-Reg}$

(c) DoD${}_{Pezenas\\_1971}^{PureSG}$

(d) DoD${}_{Pezenas\\_1971}^{Guided}$

(e) DoD${}_{Pezenas\\_1971}^{Patch}$

(f) DoD${}_{Pezenas\\_1981}^{3vGCPs}$

(g) DoD${}_{Pezenas\\_1981}^{Co-Reg}$

(h) DoD${}_{Pezenas\\_1981}^{PureSG}$

(i) DoD${}_{Pezenas\\_1981}^{Guided}$

(j) DoD${}_{Pezenas\\_1981}^{Patch}$

(k) DoD${}_{Fr\acute{e}jus\\_1954}^{3vGCPs}$

(l) DoD${}_{Fr\acute{e}jus\\_1954}^{Co-Reg}$

(m) DoD${}_{Fr\acute{e}jus\\_1954}^{PureSG}$

(n) DoD${}_{Fr\acute{e}jus\\_1954}^{Guided}$

(o) DoD${}_{Fr\acute{e}jus\\_1954}^{Patch}$

(p) DoD${}_{Fr\acute{e}jus\\_1966}^{3vGCPs}$

(q) DoD${}_{Fr\acute{e}jus\\_1966}^{Co-Reg}$

(r) DoD${}_{Fr\acute{e}jus\\_1966}^{PureSG}$

(s) DoD${}_{Fr\acute{e}jus\\_1966}^{Guided}$

(t) DoD${}_{Fr\acute{e}jus\\_1966}^{Patch}$

(u) DoD${}_{Fr\acute{e}jus\\_1970}^{3vGCPs}$

(v) DoD${}_{Fr\acute{e}jus\\_1970}^{Co-Reg}$

(w) DoD${}_{Fr\acute{e}jus\\_1970}^{PureSG}$

(x) DoD${}_{Fr\acute{e}jus\\_1970}^{Guided}$

(y) DoD${}_{Fr\acute{e}jus\\_1970}^{Patch}$

(z) DoD Legend

Figure 9: DoDs of Pezenas and Fréjus on 5 methods (3vGCPs, Co-Reg, PureSG,
Guided and Patch)

(a) DoD${}_{Kobe}^{3vGCPs}$

(b) DoD${}_{Kobe}^{Co-Reg}$

(c) DoD${}_{Kobe}^{PureSG}$

(d) DoD${}_{Kobe}^{Guided}$

(e) DoD${}_{Kobe}^{Patch}$

(f) Gd${}_{Kobe}^{3vGCPs}$

(g) Gd${}_{Kobe}^{Co-Reg}$

(h) Gd${}_{Kobe}^{PureSG}$

(i) Gd${}_{Kobe}^{Guided}$

(j) Gd${}_{Kobe}^{Patch}$

(k) Rupture map

Figure 10: (a-e) DoDs of Kobe on 5 methods (3vGCPs, Co-Reg, PureSG, Guided and
Patch). (f-j) Ground displacements (i.e., northeastward, denoted as $Gd$) of
Kobe on 5 methods. (k) Map of earthquake surface ruptures from independent
data sources.

###  4.7 Result

##### DoDs

Table 3 reports the average value and standard deviations computed on all the
DoDs in Figure 9 and 10(a-e). A dome artifact is present in DoD3vGCPs and
DoDCo-Reg for all the datasets. This kind of systematic error is known to
originate from poorly modeled camera internal parameters [34]. The results of
DoDPureSG for datasets with smaller time span (e.g. Kobe and Pezenas) are
relatively good, thanks to our roughly co-registered orientations and DSMs.
However, the results of Fréjus are not satisfactory – without adopting our
pipeline (i.e., patch matching, 3D-RANSAC filter and cross-correlation), the
inlier ratio and precision of the inter-epoch feature correspondences are low.
DoDPatch and DoDGuided performed best for all the datasets. Given the numerous
good inter-epoch feature correspondences, we are able to effectively mitigate
the dome effect. DoDPatch performed slightly better than DoDGuided on Fréjus
dataset, as the Patch is more robust when extreme scene changes are present.

##### Ground displacements

For ground displacements in Kobe, an up-lateral strike-slip movement along the
sea is present in Guided and Patch (cf. Figure 10(i-j)), but not the other 3
methods. The observed signal is coherent with the fault of the Kobe earthquake
known from independent data sources (cf. Figure 10(k)), according to which a
striking north $30^{\circ}$-$60^{\circ}$ east rupture occurred along the
northeast-southwest coastline across $\sim$18 km [65].

##### Check point accuracy

For the check point accuracy in Pezenas and Fréjus, the Patch performed best,
while the Guided ranked second best for both datasets. Also, the Patch
performed slightly better than Guided for Pezenas, while the gap widened for
Fréjus, which is consistent with the result of DoDs.

##### Robustness

Table 5 demonstrates the number of inter-epoch feature correspondences
obtained with patch matching and guided matching before and after the
3D-RANSAC and cross correlation filters. Observe that 3D-RANSAC filter and
cross correlation removed considerable number of feature correspondences, at
the same time enough feature correspondences survived, which guaranteed
robustness of our method. Moreover, the Patch recovered considerably more
feature correspondences than the Guided, which is understandable as SuperGlue
is more invariant over time than SIFT.

| Patch | Guided
---|---|---
| tentative | enhanced | final | tentative | enhanced | final
Pezenas | 1,754,434 | 1,686,486 | 284,067 | 1,671,818 | 390,438 | 129,674
Fréjus | 1,429,510 | 1,006,768 | 93,829 | 1,240,264 | 26,060 | 6,450
Kobe | 616,630 | 147,200 | 12,671 | 54,098 | 4,484 | 1,785
Table 5: Number of inter-epoch feature correspondences in each step. Tentative
stands for the feature correspondences resulting from the patch matching or
the guided matching before applying the robust filters, enhanced stands for
points after filtering with 3D-RANSAC, final means points after filtering with
cross correlation.

##  5 Conclusion

This work proposed to adopt DSMs and tiling schemes to get robust feature
correspondences between historical images taken at different epochs. It opens
up the possibility to unlock the potential of millions of historical images,
which until now have been relatively unexplored due to lack of reliable
processing techniques. Our matching strategy is able to extract reliable
inter-epoch feature correspondences under significant scene changes where SIFT
[5] and SuperGlue [25] fail. Three sets of datasets are tested, including one
that witnessed an earthquake. We showed that with the recovered feature
correspondences we are able to mitigate the camera calibration systematic
errors, leading to more accurate results, whether it is the geo-referenced DSM
or ground displacements.

Our future work will leverage depth maps to train a neural network
architecture in extracting robust features over time.

## Acknowledgment

This work was supported by ANR project DISRUPT (ANR-18-CE31-0012-0). We thank
KIGAM (Korea Institute of Geoscience and Mineral Resources) for providing Kobe
images.

## References

  * Cowley and Stichelbaut [2012] D. C. Cowley, B. B. Stichelbaut,  Historic aerial photographic archives for european archaeology,  European Journal of Archaeology 15 (2012) 217–236.
  * Giordano and Mallet [2019] S. Giordano, C. Mallet,  Archiving and geoprocessing of historical aerial images: current status in europe, official publication no 70,  in: European Spatial Data Research, 2019\.
  * USGS [2019] USGS, earthexplorer, https://earthexplorer.usgs.gov/, 2019\.
  * IGN [2019] IGN, remonterletemps, https://remonterletemps.ign.fr/, 2019\.
  * Lowe [2004] D. G. Lowe,  Distinctive image features from scale-invariant keypoints,  International journal of computer vision 60 (2004) 91–110.
  * MicMac [2021] MicMac, Github, https://github.com/micmacIGN/micmac/tree/master/src/uti_phgrm/TiePHistorical, 2021\.
  * Arandjelović and Zisserman [2012] R. Arandjelović, A. Zisserman,  Three things everyone should know to improve object retrieval,  in: 2012 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 2012, pp. 2911–2918.
  * Bursuc et al. [2015] A. Bursuc, G. Tolias, H. Jégou,  Kernel local descriptors with implicit rotation matching,  in: Proceedings of the 5th ACM on International Conference on Multimedia Retrieval, 2015, pp. 595–598.
  * Dong and Soatto [2015] J. Dong, S. Soatto,  Domain-size pooling in local descriptors: Dsp-sift,  in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 5097–5106.
  * Yu and Morel [2011] G. Yu, J.-M. Morel,  Asift: An algorithm for fully affine invariant comparison,  Image Processing On Line 1 (2011) 11–38.
  * Bay et al. [2006] H. Bay, T. Tuytelaars, L. Van Gool,  Surf: Speeded up robust features,  in: European conference on computer vision, 2006, pp. 404–417.
  * Mikolajczyk and Schmid [2004] K. Mikolajczyk, C. Schmid,  Scale & affine invariant interest point detectors,  International journal of computer vision 60 (2004) 63–86.
  * Alcantarilla et al. [2012] P. F. Alcantarilla, A. Bartoli, A. J. Davison,  Kaze features,  in: European Conference on Computer Vision, 2012, pp. 214–227.
  * Yi et al. [2016] K. M. Yi, E. Trulls, V. Lepetit, P. Fua,  Lift: Learned invariant feature transform,  in: European Conference on Computer Vision, 2016, pp. 467–483.
  * Tian et al. [2017] Y. Tian, B. Fan, F. Wu,  L2-net: Deep learning of discriminative patch descriptor in euclidean space,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 661–669.
  * Mishchuk et al. [2017] A. Mishchuk, D. Mishkin, F. Radenovic, J. Matas,  Working hard to know your neighbor’s margins: Local descriptor learning loss,  in: Advances in Neural Information Processing Systems, 2017, pp. 4826–4837.
  * Luo et al. [2019] Z. Luo, T. Shen, L. Zhou, J. Zhang, Y. Yao, S. Li, T. Fang, L. Quan,  Contextdesc: Local descriptor augmentation with cross-modality context,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2019, pp. 2527–2536.
  * Ono et al. [2018] Y. Ono, E. Trulls, P. Fua, K. M. Yi,  Lf-net: learning local features from images,  in: Advances in Neural Information Processing Systems, 2018, pp. 6234–6244.
  * Noh et al. [2017] H. Noh, A. Araujo, J. Sim, T. Weyand, B. Han,  Large-scale image retrieval with attentive deep local features,  in: Proceedings of the IEEE international conference on computer vision, 2017, pp. 3456–3465.
  * DeTone et al. [2018] D. DeTone, T. Malisiewicz, A. Rabinovich,  Superpoint: Self-supervised interest point detection and description,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, 2018, pp. 224–236.
  * Dusmanu et al. [2019] M. Dusmanu, I. Rocco, T. Pajdla, M. Pollefeys, J. Sivic, A. Torii, T. Sattler,  D2-net: A trainable cnn for joint detection and description of local features,  in: 2019 IEEE Conference on Computer Vision and Pattern Recognition, 2019, pp. 8092–8101.
  * Luo et al. [2020] Z. Luo, L. Zhou, X. Bai, H. Chen, J. Zhang, Y. Yao, S. Li, T. Fang, L. Quan,  Aslfeat: Learning local features of accurate shape and localization,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 6589–6598.
  * Revaud et al. [2019] J. Revaud, C. De Souza, M. Humenberger, P. Weinzaepfel,  R2d2: Reliable and repeatable detector and descriptor,  in: Advances in Neural Information Processing Systems, 2019, pp. 12405–12415.
  * Wiles et al. [2020] O. Wiles, S. Ehrhardt, A. Zisserman,  D2d: Learning to find good correspondences for image matching and manipulation,  arXiv preprint arXiv:2007.08480 (2020).
  * Sarlin et al. [2020] P.-E. Sarlin, D. DeTone, T. Malisiewicz, A. Rabinovich,  Superglue: Learning feature matching with graph neural networks,  in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 4938–4947.
  * Schonberger et al. [2017] J. L. Schonberger, H. Hardmeier, T. Sattler, M. Pollefeys,  Comparative evaluation of hand-crafted and learned local features,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 1482–1491.
  * Jin et al. [2020] Y. Jin, D. Mishkin, A. Mishchuk, J. Matas, P. Fua, K. M. Yi, E. Trulls,  Image matching across wide baselines: From paper to practice,  2020 IEEE Conference on Computer Vision and Pattern Recognition (2020).
  * Trulls et al. [2020] E. Trulls, Y. Jin, K. M. Yi, D. Mishkin, J. Matas, P. Fua, Image Matching Challenge 2020, https://vision.uvic.ca/image-matching-challenge/, 2020\.
  * Pinto et al. [2019] A. T. Pinto, J. A. Gonçalves, P. Beja, J. Pradinho Honrado,  From archived historical aerial imagery to informative orthophotos: A framework for retrieving the past in long-term socioecological research,  Remote Sensing 11 (2019) 1388.
  * Bożek et al. [2019] P. Bożek, J. Janus, B. Mitka,  Analysis of changes in forest structure using point clouds from historical aerial photographs,  Remote Sensing 11 (2019) 2259.
  * Persia et al. [2020] M. Persia, E. Barca, R. Greco, M. Marzulli, P. Tartarino,  Archival aerial images georeferencing: A geostatistically-based approach for improving orthophoto accuracy with minimal number of ground control points,  Remote Sensing 12 (2020) 2232.
  * Micheletti et al. [2015] N. Micheletti, S. N. Lane, J. H. Chandler,  Application of archival aerial photogrammetry to quantify climate forcing of alpine landscapes,  The Photogrammetric Record 30 (2015) 143–165.
  * Mölg and Bolch [2017] N. Mölg, T. Bolch,  Structure-from-motion using historical aerial images to analyse changes in glacier surface elevation,  Remote Sensing 9 (2017) 1021.
  * Giordano et al. [2018] S. Giordano, A. Le Bris, C. Mallet,  Toward automatic georeferencing of archival aerial photogrammetric surveys,  ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences IV-2 (2018) 105–112.
  * Dalal and Triggs [2005] N. Dalal, B. Triggs,  Histograms of oriented gradients for human detection,  in: 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 1, Ieee, 2005, pp. 886–893.
  * Feurer and Vinatier [2018] D. Feurer, F. Vinatier,  Joining multi-epoch archival aerial images in a single sfm block allows 3-d change detection with almost exclusively image information,  ISPRS journal of photogrammetry and remote sensing 146 (2018) 495–506.
  * Filhol et al. [2019] S. Filhol, A. Perret, L. Girod, G. Sutter, T. Schuler, J. Burkhart,  Time-lapse photogrammetry of distributed snow depth during snowmelt,  Water Resources Research 55 (2019) 7916–7926.
  * Cook and Dietze [2019] K. L. Cook, M. Dietze,  A simple workflow for robust low-cost uav-derived change detection without ground control points,  Earth Surface Dynamics 7 (2019) 1009–1017.
  * Parente et al. [2021] L. Parente, J. H. Chandler, N. Dixon,  Automated registration of sfm-mvs multitemporal datasets using terrestrial and oblique aerial images,  The Photogrammetric Record 36 (2021) 12–35.
  * Blanch et al. [2021] X. Blanch, A. Eltner, M. Guinau, A. Abellan,  Multi-epoch and multi-imagery (memi) photogrammetric workflow for enhanced change detection using time-lapse cameras,  Remote Sensing 13 (2021) 1460.
  * Zhang et al. [2020] L. Zhang, E. Rupnik, M. Pierrot-Deseilligny,  Guided feature matching for multi-epoch historical image blocks pose estimation,  in: ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci., 2020.
  * Maiwald and Maas [2021] F. Maiwald, H.-G. Maas,  An automatic workflow for orientation of historical images with large radiometric and geometric differences,  The Photogrammetric Record (2021).
  * Beltrami et al. [2019] C. Beltrami, D. Cavezzali, F. Chiabrando, A. Iaccarino Idelson, G. Patrucco, F. Rinaudo,  3d digital and physical reconstruction of a collapsed dome using sfm techniques from historical images.,  International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences (2019).
  * Bevilacqua et al. [2019] M. Bevilacqua, G. Caroti, A. Piemonte, D. Ulivieri,  Reconstruction of lost architectural volumes by integration of photogrammetry from archive imagery with 3d models of the status quo.,  International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences (2019).
  * Maiwald [2019] F. Maiwald,  Generation of a benchmark dataset using historical photographs for an automated evaluation of different feature matching methods.,  International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences (2019).
  * Fischler and Bolles [1981] M. A. Fischler, R. C. Bolles,  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,  Communications of the ACM 24 (1981) 381–395.
  * Sonka et al. [2014] M. Sonka, V. Hlavac, R. Boyle, Image processing, analysis, and machine vision, Cengage Learning, 2014.
  * Leroy and Rousseeuw [1987] A. M. Leroy, P. J. Rousseeuw,  Robust regression and outlier detection,  Wiley (1987).
  * Torr and Zisserman [2000] P. H. Torr, A. Zisserman,  Mlesac: A new robust estimator with application to estimating image geometry,  Computer vision and image understanding 78 (2000) 138–156.
  * Chum and Matas [2005] O. Chum, J. Matas,  Matching with prosac-progressive sample consensus,  in: 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 1, IEEE, 2005, pp. 220–226.
  * Chum et al. [2005] O. Chum, T. Werner, J. Matas,  Two-view geometry estimation unaffected by a dominant plane,  in: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), volume 1, IEEE, 2005, pp. 772–779.
  * Barath and Matas [2018] D. Barath, J. Matas,  Graph-cut ransac,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 6733–6741.
  * Barath et al. [2019] D. Barath, J. Matas, J. Noskova,  Magsac: marginalizing sample consensus,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2019, pp. 10197–10205.
  * Brachmann et al. [2017] E. Brachmann, A. Krull, S. Nowozin, J. Shotton, F. Michel, S. Gumhold, C. Rother,  Dsac-differentiable ransac for camera localization,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 6684–6692.
  * Moo Yi et al. [2018] K. Moo Yi, E. Trulls, Y. Ono, V. Lepetit, M. Salzmann, P. Fua,  Learning to find good correspondences,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 2666–2674.
  * Pierrot-Deseilligny and Cléry [2012] M. Pierrot-Deseilligny, I. Cléry,  Apero, an open source bundle adjusment software for automatic calibration and orientation of set of images,  ISPRS International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XXXVIII-5/W16 (2012) 269–276.
  * Schonberger and Frahm [2016] J. L. Schonberger, J.-M. Frahm,  Structure-from-motion revisited,  in: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 4104–4113.
  * Moulon et al. [2016] P. Moulon, P. Monasse, R. Marlet, Others, Openmvg, https://github.com/openMVG/openMVG, 2016\.
  * Sweeney [2015] C. Sweeney, Theia multiview geometry library: Tutorial & reference, http://theia-sfm.org, 2015\.
  * Pierrot-Deseilligny and Paparoditis [2006] M. Pierrot-Deseilligny, N. Paparoditis,  A multiresolution and optimization-based image matching approach: An application to surface reconstruction from spot5-hrs stereo imagery,  Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 36 (2006) 1–5.
  * Fraser [1997] C. S. Fraser,  Digital camera self-calibration,  ISPRS Journal of Photogrammetry and Remote sensing 52 (1997) 149–159.
  * Pierrot-Deseilligny et al. [2015] M. Pierrot-Deseilligny, E. Rupnik, L. Girod, J. Belvaux, G. Maillet, M. Deveau, G. Choqueux,  Micmac, apero, pastis and other beverages in a nutshell,  MicMac documentation 4 (2015).
  * Souchon et al. [2010] J.-P. Souchon, C. Thom, C. Meynard, O. Martin, M. Pierrot-Deseilligny,  The ign camv2 system,  The Photogrammetric Record 25 (2010) 402–421.
  * Rosu et al. [2015] A. Rosu, M. Pierrot-Deseilligny, A. Delorme, R. Binet, Y. Klinger,  Measurement of ground displacement from optical satellite image correlation using the free open-source software micmac,  ISPRS Journal of Photogrammetry and Remote Sensing 100 (2015) 48–59.
  * Ian and Uda [1996] A. Ian, S. Uda,  Morphological characteristics of the earthquake surface ruptures on awaji island, associated with the 1995 southern hyogo prefecture earthquake,  Island Arc 5 (1996) 1–15.