\section{Introduction}
Mass transportation problem has its origin in Monge-Kantorovich \cite{Monge1781, Kantorovich} and is stated as follows.
Let $X$ be a topological space and let $\mu_1$ and $\mu_2$ be two Borel probability measures on $X$.
Let $c\colon X^2\to [0,+\infty]$ be a function and call it a \emph{cost function}.
The Monge mass transportation problem, the Monge problem for short, with the cost function $c$ asks whether the infimum
\begin{equation}\label{Monge}
\inf_{T}\int_{X}c(x,T(x))\, d\mu_1(x) \tag{M}
\end{equation}
is achieved, where $T$ runs over all Borel measurable maps from $X$ to $X$ satisfying $T_{\sharp}\mu_1=\mu_2$.
We call such a map $T$ a \emph{transport map} from $\mu_1$ to $\mu_2$.
The measure $T_{\sharp}\mu_1$ is the push-forward measure of $\mu_1$ under $T$.
We call the objective functional in (\ref{Monge}) the \emph{total cost functional}.
A minimizer in (\ref{Monge}) is called a $c$-\emph{optimal transport map} from $\mu_1$ to $\mu_2$.
The lack of linearity and coercivity of the total cost functional makes the Monge problem difficult.
Kantorovich formulated a relaxation of the problem and overcame these difficulties.
His idea is to use so-called \emph{transport plans}, also known as joint distributions, instead of maps.
The set of transport plans, $\Pi(\mu_1,\mu_2)$, is defined by
\[
\Pi(\mu_1,\mu_2):=\left\{ \pi\in\mathcal{P}({X^2}) \mathrel{}\middle|\mathrel{} (\text{pr}_i)_{\sharp}\pi=\mu_i, i=1,2\right\},
\]
where the map $\text{pr}_i\colon X^2\to X$ is the canonical projection for $i=1,2$ and the set $\mathcal{P}(X^2)$ denotes the set of all Borel probability measures on $X^2$.
Kantorovich proposed to study the attainability of the infimum
\begin{equation}\label{Kantorovich}
\inf_{\pi\in\Pi(\mu_1,\mu_2)}\int_{X^2}c(x,y)\, d\pi(x,y). \tag{K}
\end{equation}
A minimizer in (\ref{Kantorovich}) is called a $c$-\emph{optimal transport plan} from $\mu_1$ to $\mu_2$.
If $X$ is a Polish space and if the cost function $c$ is lower semi-continuous, then an optimal transport plan  exists.
In order to obtain an optimal transport map, it is sufficient to see that a certain optimal transport plan is induced by a map.

Let  $(X,d)$ be a complete separable metric space. We focus on the case that the cost function is the distance function $d$.
Under a certain lower bound condition of curvature, the existence of an optimal transport map is proved by many researchers.
For metric spaces without non-branching geodesics, see \cite{BiCa11, BiCa13, Cava17}.
For normed spaces, see e.g., \cite{AKP04, AP03, CP11}.

In this paper, inspired by the work of Champion-Pascale \cite{CP11},
we consider the Monge problem on a metric space $(\Omega,\rho)$, where $\Omega$ is a bounded convex $G_{\delta}$-set in an Euclidean space with non-empty interior and $\rho$ is a metric on $\Omega$ satisfying some projectivity conditions.
\begin{assump*}
We assume the following three conditions on a metric $\rho$ on $\Omega$.
\begin{enumerate}
\renewcommand{\labelenumi}{{\rm(\roman{enumi})}}
\item The topology induced by $\rho$ coincides with  the Euclidean one.
\item Any line segment is a geodesic in $(\Omega,\rho)$.
\item The $n$-dimensional Lebesgue measure $\mathcal{L}^n$ on $(\Omega,\rho)$ is locally doubling.
\end{enumerate}
\end{assump*}
Note that we do not assume that $(\Omega,\rho)$ is non-branching. Under these settings, we have the following main theorem.
\begin{thm}\label{thm: main theorem}
Let $\Omega$ be a bounded convex $G_{\delta}$-set in the $n$-dimensional Euclidean space $\mathbb{R}^n$
with non-empty interior and $\rho$ be a metric on $\Omega$ satisfying {\rm Assumption A}. Let $\mu_1$ and $\mu_2$ be two Borel probability measures on $\Omega$ with compact support.
If $\mu_1$ is absolutely continuous with respect to the $n$-dimensional Lebesgue measure $\mathcal{L}^n$, then there exists a $\rho$-optimal transport map from $\mu_1$ to $\mu_2$.
\end{thm}
The idea of the proof of Theorem \ref{thm: main theorem} is based on \cite{CP11}.
We modify a variational approximation of transport plans, introduced in \cite{CP11}, and use it.
To see that our approximation scheme works well, we need an appropriate grid for a projective metric.
Unlike the case of normed spaces, it is not clear that there exists such a grid.
If we assume (iii), then this difficulty is resolved.
In this case, the doubling dimension may be different from $n$.
Note that our approximation scheme depends on the doubling dimension.

The Hilbert metric on an bounded convex open set in an Euclidean space satisfies Assumption A.
The Hilbert metric is a generalization of the Cayley-Klein model of the hyperbolic geometry.
For a bounded convex open set $\Omega$ in an Euclidean space endowed with the Hilbert metric is called the Hilbert geometry for $\Omega$.
Under some regularity assumptions on the boundary of a domain, the Hilbert geometry is regarded as a Finsler manifold of constant flag curvature $-1$.
We obtain the existence result of an optimal transport map for a Hilbert geometry as a corollary.
\begin{cor}\label{cor: Hilb geom}
Let $\Omega\subset \mathbb{R}^n$ be a bounded convex open set and $h_{\Omega}$ the Hilbert metric on $\Omega$.
Let $\mu_1$ and $\mu_2$ be two Borel probability measures on $\Omega$ with compact support.
If $\mu_1$ is absolutely continuous with respect to the $n$-dimensional Lebesgue measure $\mathcal{L}^n$, then there exists an $h_{\Omega}$-optimal transport map from $\mu_1$ to $\mu_2$.
\end{cor}
\begin{ack*}
The author would like to thank Professor Takashi Shioya for helpful comments.
\end{ack*}

\section{Preliminaries}

\subsection{General facts from optimal transport theory}
In this section, we recall some basics on optimal transport theory.
We refer to \cite{AG13, Vil09} for more details.
For a topological space $X$, denote by $\mathcal{P}(X)$ the set of all Borel probability measures on $X$.
\begin{defn
For two topological spaces $X$ and $Y$, a Borel measurable map $f\colon X\to Y$ and a Borel probability measure $\mu\in\mathcal{P}(X)$,
we define the \emph{push-forward measure} $f_{\sharp}\mu$ of $\mu$ under $f$ by
\[
f_{\sharp}\mu(B):=\mu\left(f^{-1}(B)\right)
\]
for any Borel set $B$ in $Y$.
\end{defn}
Let $\mu_1, \mu_2\in \mathcal{P}(X)$ be two Borel probability measures.
For a Borel measurable map $T\colon X\to X$, we say that $T$ is a \emph{transport map} from $\mu_1$ to $\mu_2$ if
$T_{\sharp}\mu_1=\mu_2$ holds. The set of all transport maps from $\mu_1$ to $\mu_2$ is denoted by $\mathcal{M}(\mu_1,\mu_2)$.
For a Borel probability measure $\pi\in \mathcal{P}(X^2)$, we say that $\pi$ is a \emph{transport plan} from $\mu_1$ to $\mu_2$ if
$(\text{pr}_i)_{\sharp}\pi=\mu_i$ for $i=1,2$ holds, where the map $\text{pr}_i\colon X^2\to X$ is the canonical projection for $i=1,2$.
In this case, we say that $\mu_1$ is \emph{the first marginal of $\pi$} and that $\mu_2$ is \emph{the second marginal of $\pi$}.
The set of all transport plans from $\mu_1$ to $\mu_2$ is denoted by $\Pi(\mu_1,\mu_2)$.
For a transport map $T\in\mathcal{M}(\mu_1,\mu_2)$, we denote $(\text{id}, T)_{\sharp}\mu_1$ by $\pi_T$, where $(\text{id}, T)\colon X\to X^2$ is the map that assigns each $x\in X$ to $(x, T(x))$.

\begin{defn}
Let $c\colon X^2\to [0,+\infty]$ be a function, which we call a \emph{cost function}. A transport map $T_0\in \mathcal{M}(\mu,\mu_2)$ is a \emph{$c$-optimal transport map} if it minimizes
\[
I_{\mathrm{M}}(T):=\int_{X}c(x,T(x))\, d\mu(x),\, T\in\mathcal{M}(\mu,\mu_2).
\]
A transport plan $\pi_0\in\Pi(\mu_1,\mu_2)$ is a \emph{$c$-optimal transport plan} if it minimizes
\[
I_{\mathrm{K}}(\pi):= \int_{X^2}c(x,y)\, d\pi(x,y),\, \pi\in \Pi(\mu_1,\mu_2).
\]
\end{defn}
The minimization problem for $I_M$ and $I_K$ is called the Monge problem and the Kantorovich problem for the cost function $c$, respectively.
By a direct argument, the existence of a $c$-optimal transport plan is ensured if $X$ is a Polish space and if $c$ is lower semi-continuous.
Kantorovich problem admits a dual formulation. We recall the notion of $c$-transform and $c$-concavity.
\begin{defn
Let $\varphi\colon X\to [-\infty,+\infty]$ be a function.
\begin{enumerate}
\item The \emph{$c$-transform} $\varphi^c$ of $\varphi$ is defined by
\[
\varphi^c(y):= \inf_{x\in X}(c(x,y)-\varphi(x)),\, y\in X.
\]
\item $\varphi$ is \emph{$c$-concave} if there exists a function $\psi\colon X\to[-\infty,+\infty]$ on $X$ such that $\varphi=\psi^c$ holds.
\end{enumerate}
\end{defn}

\begin{example}
Let $(X,d)$ be a metric space and $\varphi\colon X\to\mathbb{R}$ a function.
Then, $\varphi$ is $d$-concave if and only if $\varphi$ is $1$-Lipschitz continuous.
\end{example}

\begin{defn
A subset $\Gamma\subset X^2$ is said to be \emph{c-cyclically monotone} if
for any finitely many  points $(x_1,y_1),\ldots,(x_n,y_n)\in \Gamma$, we have
\[
\sum_{i=1}^n c(x_i,y_i)\leq \sum_{i=1}^n c(x_i,y_{i+1}),
\]
where $y_{n+1}:=y_1$. A Borel probability measure $\pi$ on $X^2$ is said to be \emph{$c$-cyclically monotone}
if $\pi$ is concentrated on a $c$-cyclically monotone set.
\end{defn}

In the case that $(X,d)$ is the $n$-dimensional Euclidean space $\mathbb{R}^n$ and $c$ is the square of the Euclidean metric, the $c$-cyclical monotonicity of $\Gamma$ yields the monotonicity of $\Gamma$, i.e., for any $(x,y), (x',y')\in \Gamma$, we have
\[
(y-y')\cdot (x-x') \geq 0,
\]
where the dot is the Euclidean inner product of $\mathbb{R}^n$.
Next we recall the notion of Kantorovich potential.
Denote by $\Phi_c$ the set of all pairs of functions $(\varphi,\psi)\in L^1(X,\mu_1)\times L^1(X,\mu_2)$ such that
\[
\varphi(x)+\psi(y)\leq c(x,y)
\]
holds for $\mu_1$-almost all $x\in X$ and $\mu_2$-almost all $y\in X$.
Denote by $\Psi_c$ the set of all pairs of functions $(\varphi,\psi)\in C_{b}(X)\times C_{b}(X)$ such that
\[
\varphi(x)+\psi(y)\leq c(x,y)
\]
holds for all $x, y\in X$.
For $(\varphi,\psi)\in\Phi_c$, we define $J(\varphi,\psi)$ by
\[
J(\varphi,\psi):=\int_{X}\varphi\, d\mu_1 + \int_{X}\psi\, d\mu_2.
\]
We call every $(\varphi,\psi)\in\Phi_c$ maximizing  $J$ a \emph{Kantorovich potential}.
The following theorem is known as the Kantorovich duality.
\begin{thm}[Kantorovich duality]\label{KRduality}
Let $c\colon X^2 \to [0,+\infty]$ be a lower semi-continuous function. Then we have the equalities
\[
\min_{\pi\in\Pi(\mu_1,\mu_2)}I_{\mathrm{K}}(\pi)
=\sup_{(\varphi,\psi)\in \Phi_c}J(\varphi,\psi)
=\sup_{(\varphi,\psi)\in \Psi_c}J(\varphi,\psi).
\]
\end{thm}

\begin{thm}
Let $c\colon X^2\to [0,+\infty]$ be a lower semi-continuous function.
If $\pi\in\Pi(\mu_1,\mu_2)$ is optimal and if $I_{\mathrm{K}}(\pi)$ is finite, then $\pi$ is concentrated on a $c$-cyclically monotone set.
Moreover, if $c$ is real-valued, then there exists a $c$-cyclically monotone set $\Gamma\subset X^2$ such that for any $\pi\in\Pi(\mu_1,\mu_2)$, the following are equivalent to each other.
\begin{enumerate}
\item $\pi$ is a $c$-optimal transport plan from $\mu_1$ to $\mu_2$.
\item $\pi$ is $c$-cyclically monotone.
\item There exists a $c$-concave Borel measurable function $\varphi\colon X\to \mathbb{R}$ such that
for $\pi$-almost every point \ $(x,y)\in X\times Y$, we have
\[
\varphi(x)+\varphi^c(y)=c(x,y).
\]
\item $\pi$ is concentrated on $\Gamma$.
\end{enumerate}
\end{thm}

\begin{thm}\label{thm: potential}
Let $c\colon X^2\to [0,+\infty)$ be a lower semi-continuous function.
Assume that $\mathcal{I}_c(\mu_1,\mu_2)<+\infty$,
\[
\mu_1\left(\left\{ x\in X \mathrel{}\middle|\mathrel{} \int_Y c(x,y)\, d\mu_2(y) < +\infty\right\}\right) > 0,
\]
and
\[
\mu_2\left(\left\{ y\in X \mathrel{}\middle|\mathrel{} \int_X c(x,y)\, d\mu_1(x) < +\infty\right\}\right) > 0.
\]
Then, there exists a $c$-concave function $\varphi\colon X \to \mathbb{R}$ such that the pair $(\varphi, \varphi^c)$ is a Kantorovich potential.
\end{thm}

Let $(X,d)$ be a complete and separable metric space and let $p\geq 1$.
We say that a measure $\mu\in \mathcal{P}(X)$ has \emph{finite $p$-th moment with respect to $d$} if there exists a point $x_0\in X$ such that
\[
\int_X d^p(x,x_0)\, d\mu(x)<+\infty.
\]
Denote by $\mathcal{P}_p(X,d)$ the set of all probability measures with finite $p$-th moment with respect to $d$.
Note that for $\mu_1,\mu_2\in \mathcal{P}_p(X)$, the assumptions in Theorem \ref{thm: potential} are fulfilled for $c=d^p$.
\emph{The $L^p$-Wasserstein metric} $W_p$ on $\mathcal{P}(X)$ is defined by
\[
W_p(\mu_1,\mu_2):=\inf\left\{\left(\int_{X^2}d^p(x,y)\, d\pi(x,y)\right)^{1/p} \mathrel{}\middle|\mathrel{} \pi\in\Pi(\mu_1,\mu_2)\right\},
\]
i.e., $W_p(\mu_1,\mu_2)$ is the $p$-th root of the $d^p$-optimal transport cost from $\mu_1$ to $\mu_2$. The function $W_p$ is indeed a metric on $\mathcal{P}_p(X,d)$.

\subsection{Some facts on metric spaces}
In this section, we enumerate some facts and prove some claims on metric spaces.

Let $(X,d)$ be a metric space.
We call a discrete subset of $X$ a \emph{net}.
For a net $\mathcal{N}\subset X$ and $\varepsilon>0$,
we say that $\mathcal{N}$ is an $\varepsilon$-\emph{net} if the $\varepsilon$-neighborhood $B_{\varepsilon}(\mathcal{N})$ of $\mathcal{N}$
coincides with the whole $X$.
We say that $\mathcal{N}$ is $\varepsilon$-\emph{discrete} if $d(x,y)>\varepsilon$ for any two distinct points $x,y\in \mathcal{N}$.
The set of all $\varepsilon$-discrete nets equipped with the inclusion relation forms a poset and we see that  any maximal element of the set is an $\varepsilon$-net.
\begin{defn}
Let $\mathfrak{m}$ be a non-negative Borel measure on $X$.
We say that $\mathfrak{m}$ is \emph{locally doubling} if for any $R>0$, there exists a constant $C=C_R\geq 1$ such that
$\mathfrak{m}(B(x,2r))\leq C\mathfrak{m}(B(x,r))$ for any $x\in X$ and $0<r\leq R$.
\end{defn}
The locally doubling condition is equivalent to
\begin{equation}\label{ineq: doubling}
\mathfrak{m}(B(x,r'))\leq C^{\lceil \log_2(r'/r) \rceil}\mathfrak{m}(B(x,r))
\end{equation}
for any $x\in X$ and $0<r<r'\leq R$,
where $\log_2$ is the base $2$ logarithm function and $\lceil\cdot\rceil$ is the ceil function.
\begin{lem}\label{lem: numberofnets}
Let $X$ be a metric space, let $B\subset X$ be a bounded Borel set with diameter $D$ and let $\mathcal{N}\subset B$ be an $\varepsilon$-discrete net of $B$. Assume that $X$ admits a locally doubling Borel measure $\mathfrak{m}$. Then the number $\#\mathcal{N}$ of $\mathcal{N}$ is bounded above by $C^{\lceil \log_2(2D/\varepsilon)\rceil}$.
\end{lem}

\begin{proof}
The $\varepsilon$-discreteness of $\mathcal{N}$ yields that a family of $\varepsilon/2$-balls \\
$\{B(x,\varepsilon/2) \mid x\in \mathcal{N}\}$ is pairwise disjoint. Combining this fact with the doubling condition (\ref{ineq: doubling}) for $R=D$ and $r=\varepsilon/2$, we obtain
\[
\mathfrak{m}(B)\geq \mathfrak{m}\left(\bigcup_{x\in \mathcal{N}}B(x,\varepsilon/2)\right)
 = \sum_{x\in \mathcal{N}} \mathfrak{m}(B(x,\varepsilon/2))
\geq \frac{\#\mathcal{N}\mathfrak{m}(B)}{C^{\lceil \log_2 (2D/\varepsilon)\rceil}},
\]
which yields the lemma.
\end{proof}

\begin{defn}
Let $A\subset X$ be any subset of $X$. A map $p\colon X\to A$ is called a \emph{nearest point projection to} $A$
if $d(x,p(x)) = d(x,A)$ for any $x\in X$, where $d(x,A)=\inf_{a\in A} d(x,a)$.
\end{defn}

\begin{lem}
If $A=\{a_1,\ldots,a_n\}\subset X$ is a finite subset, then there exists a Borel measurable nearest point projection to $A$.
\end{lem}

\begin{proof}
We put
\[
i(x):=\min\{i\mid d(x,a_i)\leq d(x,a_j)\ \text{for any}\ j\}
\]
and $p(x):=a_{i(x)}$. The assertion follows from
\[
p^{-1}(a_i)=\left(\bigcap_{j=1}^n\{f_i\leq f_j\}\right)\cap \left(\bigcap_{j<i}\bigcup_{k=1}^n\{f_j > f_k\}\right), i=1,\ldots,n,
\]
where $f_i$ is the distance function from $a_i$.
\end{proof}

\section{Proof of Main Theorem}
In order to prove the existence of an optimal transport map, it suffices to show that a certain optimal transport plan is induced by a transport map.
The following lemma gives a criterion according to which a measure in a product space is induced by a map.
\begin{lem}\label{lem: map}
Let $X$ and $Y$ be two Polish spaces and let $\pi\in \mathcal{P}(X\times Y)$.
Let $\mu_1$ be the first marginal of $\pi$.
Then, $\pi$ is induced by a map if and only if
there exists a Borel measurable set $\Gamma\subset X\times Y$ of full $\pi$-measure such that
for $\mu_1$-almost every $x\in X$ there exists a unique $y=T(x)\in Y$ with $(x,y)\in \Gamma$.
In this case, the map $T$ is $\mu_1$-measurable and $\pi=({\rm id}, T)_{\sharp}\mu_1$.
\end{lem}

Throughout this section, let $\Omega\subset\mathbb{R}^n$ be a bounded and convex $G_{\delta}$-set with non-empty interior and let $\rho$ be a metric on $\Omega$ satisfying Assumption A (see Section 1).
The idea that we treat $G_{\delta}$-sets is based on Mazurkiewictz's theorem (for the proof, see e.g.\cite{Wil70}).
\begin{thm}[Mazurkiewictz]
Let $X$ be a completely metrizable topological space and $A$ a subspace of $X$.
Then $A$ is completely metrizable if and only if $A$ is a $G_{\delta}$-set.
\end{thm}
Let $\lambda$ be a Borel probability measure on $\Omega$ that is absolutely continuous with respect to the $n$-dimensional Lebesgue measure $\mathcal{L}^n$.
We call the function $f$ defined by
\[
f(x):= \limsup_{r\to 0+}\frac{\lambda(B(x,r))}{\mathfrak{m}(B(x,r))}
\]
the \emph{density} of $\lambda$
and denote it by $d\lambda/d\mathcal{L}^n$.

For two probability measures $\mu_1,\mu_2\in \mathcal{P}(\Omega)$, let $\mathcal{O}_{0}(\mu_1,\mu_2)$ be the set of all $\rho$-optimal transport plans from $\mu_1$ to $\mu_2$.
We construct two families $\mathcal{O}_1(\mu_1,\mu_2)$ and $\mathcal{O}_2(\mu_1,\mu_2)$ of optimal transport plans and show that each of elements in $\mathcal{O}_2(\mu_1,\mu_2)$ is induced by a map. As a result, the set $\mathcal{O}_2(\mu_1,\mu_2)$ consists of a single element.
The construction of $\mathcal{O}_2(\mu_1,\mu_2)$ is similar to that of \cite{San09} and \cite{CP11}.
We first define the set $\mathcal{O}_1(\mu_1,\mu_2)$.
We consider the secondary variational problem:
\[ \label{svp1}
\inf \left\{ \int_{\Omega^2}\lvert x-y \rvert^2\, d\pi \mathrel{}\middle|\mathrel{} \pi \in \mathcal{O}_0(\mu_1,\mu_2)\right\}. \tag{SVP1}
\]
We denote by $\mathcal{O}_1(\mu_1,\mu_2)$ the set of all solutions of (\ref{svp1}).
Since the topology of $(\Omega,\rho)$ is Euclidean, it is easy to see that $\mathcal{O}_0(\mu_1,\mu_2)$ is non-empty.
On the other hand, we cannot deduce the monotonicity of a measure in $\mathcal{O}_1(\mu_1,\mu_2)$ due to the lack of elements in the feasible region in (\ref{svp1}). Nevertheless, we will see that a weak version of monotonicity holds for any measure in $\mathcal{O}_1(\mu_1,\mu_2)$.
We say that a subset $\Gamma$ of $\Omega^2$ is \emph{restrictedly monotone} if for any $(x,y), (x',y')\in \Gamma$ with $x\in [x',y']$, we have
\[
(y-y')\cdot (x-x')\geq 0,
\]
where the set $[x',y']\subset \Omega$ is the line segment between $x'$ and $y'$.
A Borel measure $\pi$ on $\Omega^2$ is said to be \emph{restrictedly monotone} if $\pi$ is concentrated on a restrictedly monotone subset of $\Omega^2$.
We will see that any measure in $\mathcal{O}_1(\mu_1,\mu_2)$ is restrictedly monotone.
We replace the problem (\ref{svp1}) with an equivalent one whose feasible region is the whole $\Pi(\mu_1,\mu_2)$.
Accordingly, we must change the objective functional. We fix a Kantorovich potential $\varphi$ for $W_1(\mu_1,\mu_2)$ and
put
\[
\beta(x,y) =
\begin{cases}
\lvert x-y\rvert^2 & \varphi(x)-\varphi(y)=\rho(x,y), \\
+\infty & otherwise.
\end{cases}
\]
The ordinary optimal transport problem with cost function $\beta$
\[ \label{svp2}
\inf \left\{ \int_{\Omega^2}\beta(x,y)\, d\pi \mathrel{}\middle|\mathrel{} \pi \in \Pi(\mu_1,\mu_2)\right\} \tag{SVP2}
\]
is equivalent to the secondary variational problem (\ref{svp1}). More precisely, we have the following lemma.
\begin{lem}\label{equivalence}
Let $\mu_1, \mu_2, \varphi$ and $\beta$ be as above.
Then, for any $\pi\in\mathcal{O}_0(\mu_1,\mu_2)$, $\pi$ is a solution of {\rm (SVP1)} if and only if it is a solution of {\rm (SVP2)}.
\end{lem}
By virtue of Lemma \ref{equivalence}, we see that for any $\pi\in \mathcal{O}_1(\mu_1,\mu_2)$,
there exists a $\beta$-cyclically monotone set $\Gamma$ on which $\pi$ concentrates.
We also observe that for any $\pi\in \mathcal{O}_{0}(\mu_1,\mu_2)$, $\beta$ takes finite values on $\Gamma$.
Furthermore, the $\beta$-cyclical monotonicity is a weak version of the monotonicity.
More generally, we have the following. In the proof, we use the condition on  geodesics for $\rho$.
\begin{prop}\label{lem: monotone}
Let $\varphi\colon \Omega\to\mathbb{R}$ be a $1$-Lipschitz continuous function.
Define a extended-real valued function $\beta\colon \Omega^2\to [0,+\infty]$ by
\[
\beta(x,y):=
\begin{cases}
\lvert x-y \rvert^2 & {\rm if}\ \varphi(x)-\varphi(y)=\rho(x,y), \\
+\infty & otherwise.
\end{cases}
\]
Let $\Gamma\subset \Omega^2$ be a $\beta$-cyclically monotone set on which $\beta$ is finite.
Then $\Gamma$ is restrictedly monotone.
\end{prop}
\begin{proof}
By the $\beta$-cyclical monotonicity of $\Gamma$ and the finiteness of $\beta$, we have
\[
\lvert x-y \rvert^2 + \lvert x'- y' \rvert^2 \leq \beta(x,y') + \beta(x',y).
\]
Due to the $1$-Lipschitz continuity of $\varphi$, it suffices to prove that
\begin{equation}\label{lem: monotonicity ineq 1}
\varphi(x)-\varphi(y') \geq \rho(x,y')
\end{equation}
and
\begin{equation}\label{lem: monotonicity ineq 2}
\varphi(x')-\varphi(y) \geq \rho(x',y).
\end{equation}
Note that the argument is not symmetric.
We first prove (\ref{lem: monotonicity ineq 2}).
Since $x', x$ and $y$ are aligned, we have
\begin{equation}\label{eq: align}
\rho(x', y') = \rho(x', x) + \rho(x, y').
\end{equation}
Combining (\ref{eq: align}) with the $1$-Lipschitz continuity of $\varphi$ yields
\begin{align*}
\varphi(x')-\varphi(x) &\leq \rho(x,x') \\
                             &= \rho(x',y') - \rho(x,y') \\
                             &\leq \varphi(x')-\varphi(y')-(\varphi(x)-\varphi(y')) \\
                             &= \varphi(x')-\varphi(x).
\end{align*}
In particular, we have
\begin{equation}\label{eq: aaa}
\varphi(x') - \varphi(x) = \rho(x',x).
\end{equation}
Combining (\ref{eq: aaa}) with the triangle inequality, we arrive at
\begin{align*}
\rho(x',y) &\leq \rho(x',x) + \rho(x,y) \\
                      &= \varphi(x')-\varphi(x)+\varphi(x)-\varphi(y) \\
                      &= \varphi(x')-\varphi(y),
\end{align*}
which proves (\ref{lem: monotonicity ineq 2}).
We next prove (\ref{lem: monotonicity ineq 1}). By using the triangle inequality for $x', x$ and $y$, we obtain
\begin{align*}
\rho(x,y') &\leq \rho(x,y')+\rho(x',x)+\rho(x,y)-\rho(x',y)  \\
                      &= \rho(x,y)+\rho(x',y')-\rho(x',y) \\
                      &= \varphi(x)-\varphi(y)+\varphi(x')-\varphi(y')-(\varphi(x')-\varphi(y)) \\
                      &= \varphi(x)-\varphi(y').
\end{align*}
The proof is completed.
\end{proof}

Next, we construct a subset $\mathcal{O}_2(\mu_1,\mu_2)$ of $\mathcal{O}_1(\mu_1,\mu_2)$, which plays a crucial role in the proof of the main theorem.
We assume that the support of $\mu_1$ and $\mu_2$ is compact, respectively.
We denote by $\mathcal{P}_c(\Omega)$ the set of all compactly supported Borel probability measures on $\Omega$.
Let $\varepsilon>0$ be a positive real number. We define a functional $C_{\varepsilon}\colon \mathcal{P}(\Omega^2)\to [0,+\infty]$ by
\[
C_{\varepsilon}(\pi):=\frac{1}{\varepsilon}W_1((\text{pr}_2)_{\sharp}\pi,\mu_2) + \int_{\Omega^2}c_{\varepsilon}\, d\pi+\varepsilon^{3d+2}\#\mathop{\rm supp}\nolimits((\text{pr}_2)_{\sharp}\pi),
\]
where $d$ is the the base 2 logarithm of the doubling constant of $\mathop{\rm supp}\nolimits\mu_2$ and $c_{\varepsilon}(x,y):=\rho(x,y)+\varepsilon\lvert x-y\rvert^2$, $x,y\in \Omega$.
The functional $C_{\varepsilon}$ is lower semi-continuous with respect to the weak convergence of measures.
We denote by $D_{\varepsilon}$ the set of all minimizers of $C_{\varepsilon}$ whose first marginal is $\mu_1$:
\[
D_{\varepsilon}:=\arg\min\left\{ C_{\varepsilon}(\pi) \mathrel{}\middle|\mathrel{} \pi\in \mathcal{P}(\Omega^2), (\text{pr}_1)_{\sharp}\pi=\mu_1 \right\}.
\]
We observe that the set $D_{\varepsilon}$ is not empty.
$\mathcal{O}_2(\mu_1,\mu_2)\subset\Pi(\mu_1,\mu_2)$ is defined to be the set of all cluster points of any sequence $\{\pi_{\varepsilon}\in D_{\varepsilon}\}_{\varepsilon>0}$.
By the compactness of the support of $\mu_1$ and $\mu_2$, we see that $\mathcal{O}_2(\mu_1,\mu_2)$ is not empty.
Furthermore, we see that any measure in $\mathcal{O}_2(\mu_1,\mu_2)$ solves (\ref{svp1}).
\begin{lem}\label{well-definedness of O2}
For any $\varepsilon>0$, let $\pi_{\varepsilon}\in D_{\varepsilon}$.
Then the second marginal of  $\pi_{\varepsilon}$ converges to $\mu_2$ weakly as $\varepsilon\to 0$.
Moreover, every limit point $\pi_*$ of $\{\pi_{\varepsilon}\}_{\varepsilon>0}$ is an element of $\mathcal{O}_1(\mu_1, \mu_2)$.
\end{lem}

\begin{proof}
For every natural number $j$, we take a maximal $(1/j)$-net $\mathcal{N}_j$ of  $\mathop{\rm supp}\nolimits\mu_2$.
Put $d:=\log_2 C$. By Lemma \ref{lem: numberofnets}, we have
$\#\mathcal{N}_j\leq O(j^d)$ for any $j$. Since $\mathcal{N}_j$ is finite, there exists a Borel measurable nearest point projection
$p_j\colon \mathop{\rm supp}\nolimits\mu_2\to\mathcal{N}_j$.
Pick and fix a transport plan $\pi\in \mathcal{O}_1(\mu_1,\mu_2)$. Set $\pi'_j:=(\text{id}\times p_j)_{\sharp}\pi$, where the map $\text{id}\times p_j$ is defined by
$(\text{id}\times p_j)(x,x'):=(x,p_j(x'))$. The measure $\pi'_j$ is a transport plan from $\mu_1$ to $\mu_2$.
Since $\pi_{\varepsilon}\in D_{\varepsilon}$, we have
\begin{align}
C_{\varepsilon}(\pi_{\varepsilon})
&= \frac{1}{\varepsilon}W_1((\text{pr}_2)_{\sharp}\pi_{\varepsilon},\mu_2)+\int_{\Omega^2}c_{\varepsilon}\, d\pi_{\varepsilon}+\varepsilon^{3d+2}\#\mathop{\rm supp}\nolimits(\text{pr}_2)_{\sharp}\pi_{\varepsilon} \label{ineq: masterineq} \\
&\leq \frac{1}{\varepsilon}W_1((\text{pr}_2)_{\sharp}\pi'_j,\mu_2)+\int_{\Omega^2}c_{\varepsilon}\, d\pi'_j+\varepsilon^{3d+2}\#\mathop{\rm supp}\nolimits(\text{pr}_2)_{\sharp}\pi'_j. \nonumber
\end{align}
Since the map $p_j$ is a nearest point projection to $\mathcal{N}_j$, we have
\[
W_1((\text{pr}_2)_{\sharp}\pi'_j, \mu_2)=W_1((p_j)_{\sharp}\mu_2,\mu_2) \leq \int_{\Omega}\rho(x,p_j(x))\, d\mu_2(x)\leq \frac1j
\]
and
\[
\#\mathop{\rm supp}\nolimits(\text{pr}_2)_{\sharp}\pi'_j=\mathop{\rm supp}\nolimits(p_j)_{\sharp}\mu_2\leq \#\mathcal{N}_j,
\]
so we obtain
\[
C_{\varepsilon}(\pi_{\varepsilon})\leq \frac{1}{j\varepsilon}+\int_{\Omega^2}c_{\varepsilon}\, d\pi'_j + \varepsilon^{3d+2}\#\mathcal{N}_j.
\]
In particular, this yields that
\[
\frac{1}{\varepsilon}W_1((\text{pr}_2)_{\sharp}\pi_{\varepsilon},\mu_2)\leq C_{\varepsilon}(\pi_{\varepsilon}) \leq \frac1{j\varepsilon} + \int_{\Omega^2}c_{\varepsilon}\, d\pi'_j +
\varepsilon^{3d+2}\#\mathcal{N}_j.
\]
By multiplying $\varepsilon>0$ and letting $\varepsilon\to 0$, we have
\[
\limsup_{\varepsilon\to 0}W_1((\text{pr}_2)_{\sharp}\pi_{\varepsilon},\mu_2) \leq \frac1j.
\]
Since $j$ is arbitrary, $(\text{pr}_2)_{\sharp}\pi_{\varepsilon}$ converges to $\mu_2$ weakly as
$\varepsilon\to 0$.
By letting $j_{\varepsilon}\approx \varepsilon^{-2}$, (\ref{ineq: masterineq}) yields
\[
\int_{\Omega^2}\rho\, d\pi_{\varepsilon}
\leq \int_{\Omega^2}c_{\varepsilon}\, d\pi_{\varepsilon}
\leq \int_{\Omega^2}\rho\, d\pi'_{j_{\varepsilon}}+O(\varepsilon).
\]
Letting $\varepsilon\to 0$ yields
\[
\int_{\Omega^2}\rho\, d\pi_* \leq \liminf_{\varepsilon\to 0}\int_{\Omega^2}\rho\, d\pi_{\varepsilon}\leq W_1(\mu_1,\mu_2).
\]
Since $\pi_*\in\Pi(\mu_1,\mu_2)$, we have $\pi_*\in \mathcal{O}_0(\mu_1,\mu_2)$. Moreover, by using
\[
\int_{\Omega^2}\rho\, d\pi_{\varepsilon} \geq W_1(\mu_1,(\text{pr}_2)_{\sharp}\pi_{\varepsilon})\geq W_1(\mu_1,\mu_2)-W_1(\mu_2,(\text{pr}_2)_{\sharp}\pi_{\varepsilon}),
\]
\begin{align*}
\int_{\Omega^2}\rho(x,y)\, d\pi'_j
&= \int_{\Omega^2}\rho(x,p_j(y))\, d\pi \\
&\leq \int_{\Omega^2}\rho(x,y)\, d\pi +\int_{\Omega^2}\rho(y,p_j(y)) \leq W_1(\mu_1,\mu_2)+\frac1j
\end{align*}
and (\ref{ineq: masterineq}), we get
\[
\frac{1-\varepsilon}{\varepsilon}W_1(\mu_2,(\text{pr}_2)_{\sharp}\pi_{\varepsilon})+\varepsilon\int_{\Omega^2}\lvert x-y\rvert^2\, d\pi_{\varepsilon}
\leq \frac{1+\varepsilon}{j\varepsilon}+\varepsilon\int_{\Omega^2}\lvert x-y\rvert^2\, d\pi'_j +O(\varepsilon^{3d+2}).
\]
In particular, by letting $\varepsilon<1$ and $j_{\varepsilon}\approx\varepsilon^{-3}$, we obtain
\[
\int_{\Omega^2}\lvert x-y\rvert^2\, d\pi_{\varepsilon} \leq \int_{\Omega^2}\lvert x-y\rvert^2\, d\pi'_{j_{\varepsilon}} + O(\varepsilon^2),
\]
which leads us to
\[
\int_{\Omega^2}\lvert x-y\rvert^2\, d\pi_* \leq \int_{\Omega^2}\lvert x-y\rvert^2\, d\pi.
\]
Since $\pi\in\mathcal{O}_1(\mu_1,\mu_2)$, we have $\pi_*\in\mathcal{O}_1(\mu_1,\mu_2)$.
\end{proof}

We consider restricted measures and interpolated measures.
For a Borel measure $\pi$ on $\Omega^2$ and a Borel set $B\subset \Omega^2$ of positive measure,
we denote by $\pi\lfloor_B$ the restriction of $\pi$ to $B$.
We interpolate measures along line segments in $\Omega$. For $0\leq t\leq 1$, we put
\[
P_t(x,y):=(1-t)x+ty,\, x,y\in\Omega.
\]
\begin{lem}\label{lem: densitybound}
Let $B\subset \Omega^2$ be a Borel set, let $\varepsilon>0$ and $\pi_{\varepsilon} \in D_{\varepsilon}$.
Denote by $\mu_1^{\varepsilon ,B}$ and $\mu_2^{\varepsilon,B}$ the first and the second marginal of $\pi_{\varepsilon}\lfloor_B$ respectively.
Then, $\mu_1^{\varepsilon, B}$ is absolutely continuous with respect to $\mathcal{L}^n$ and we have
the following {\rm (1), (2)} and {\rm (3)}.
\begin{enumerate}
\item The measure $\pi_{\varepsilon}\lfloor_B$ is a $c_{\varepsilon}$-optimal transport plan from $\mu_1^{\varepsilon, B}$ to
$\mu_2^{\varepsilon, B}$.
\item For any $t\in [0,1)$, the interpolated measure $(P_t)_{\sharp}(\pi_{\varepsilon}\lfloor_B)$ is absolutely continuous with respect to $\mathcal{L}^n$.
\item If the density $\rho_{\varepsilon, B}$ of $\mu_1^{\varepsilon,B}$ is essentially bounded, then
so is the density of $(P_t)_{\sharp}(\pi_{\varepsilon}\lfloor_{B})$. In this case, we have
\[
\left\|\frac{d(P_t)_{\sharp}(\pi_{\varepsilon}\lfloor_B)}{d\mathcal{L}^n}\right\|_{L^\infty} \leq (1-t)^{-n}\|\rho_{\varepsilon,B}\|_{L^\infty}.
\]
\end{enumerate}
\end{lem}

\begin{proof}
The absolute continuity of $\mu_1^{\varepsilon,B}$ is clear.
The assertion (i) follows from the stability of optimality under the restriction of measures.
Let us prove (ii). If $t=0$, then the assertion is trivial. We may assume $t>0$.
Let $\{y_i\}_{i\in I}$ be the support of $\mu_2^{\varepsilon,B}$ and put $\Omega_i:=\mathop{\rm supp}\nolimits(\pi_{\varepsilon}\lfloor_{\Omega\times\{y_i\}})$ and $\Omega_i(t):= P_t(\Omega\times \{y_i\})$. By a simple calculation, for any Borel set $A\subset \Omega$,
\begin{align}
(P_t)_{\sharp}(\pi_{\varepsilon}\lfloor_B)(A) \nonumber
&= \pi_{\varepsilon}\lfloor_B(P_t^{-1}(A)) \nonumber \\
&\leq \sum_{i\in I} \pi_{\varepsilon}\lfloor_B\left(P_t^{-1}(A\cap \Omega_i(t))\right) \nonumber \\
&= \sum_{i\in I} \pi_{\varepsilon}\lfloor_B\left(\frac{1}{1-t}(A\cap \Omega_i(t)-ty_i)\times \Omega\right) \nonumber \\
&= \sum_{i\in I}  \mu_1^{\varepsilon,B}\left(\frac{1}{1-t}(A\cap \Omega_i(t)-ty_i)\right) \nonumber \\
&\leq (1-t)^{-n}\sum_{i\in I} \|\rho_{\varepsilon,B}\|_{L^{\infty}}\mathcal{L}^n(A\cap\Omega_i(t)) \label{ineq: densitybound},
\end{align}
where we have used the translation invariance of $\mathcal{L}^n$ in the last inequality.
We now claim that the family $\{\Omega_i(t)\}_{i\in I}$ is pairwise disjoint for any $0<t<1$.
Assume the contrary.
Then, there exist two points $x_i,x_j\in \Omega$ such that
\[
(1-t)x_i+ty_i=(1-t)x_j+y_j.
\]
From this, we have
\begin{align*}
x_i-y_j
&= (1-t)x_j-(1-t)y_j+t(x_i-y_i) \\
&= (1-t)(x_j-y_j)+t(x_i-y_i)
\end{align*}
and
\[
c_{\varepsilon}(x_i,y_j) < (1-t)c_{\varepsilon}(x_j,y_j)+tc_{\varepsilon}(x_i,y_i)
\]
by the strict convexity of $c_{\varepsilon}$. Exchanging $i$ and $j$, we also have
\[
c_{\varepsilon}(x_j,y_i) < (1-t)c_{\varepsilon}(x_i,y_i)+tc_{\varepsilon}(x_j,y_j).
\]
On the other hand, the $c_{\varepsilon}$-cyclical monotonicity of $\mathop{\rm supp}\nolimits\pi_{\varepsilon}$ yields
\[
c_{\varepsilon}(x_i,y_i)+c_{\varepsilon}(x_j,y_j)
\leq
c_{\varepsilon}(x_i,y_j)+c_{\varepsilon}(x_j,y_i).
\]
This is a contradiction and the assertion (iii) follows.
\end{proof}

\begin{lem}\label{lem: restriction}
Let $\pi_k,\pi\in \mathcal{P}(\Omega^2)$, $k=1,2,\ldots$, be Borel probability measures with a common first marginal.
If $\{\pi_k\}_{k=1}^{\infty}$ converges to $\pi$ weakly, then for any Borel set $G\subset \Omega$, the sequence $\{\pi_k\lfloor_{G\times \Omega}\}_{k=1}^{\infty}$ converges to $\pi\lfloor_{G\times \Omega}$ weakly.
\end{lem}

\begin{defn}
Let $\Gamma\subset \Omega^2$ be a subset.
We put
\[
\mathbb{T}(\Gamma):=\left\{P_t(x,y)\mathrel{}\middle|\mathrel{} t\in [0,1], (x,y)\in \Gamma\right\}
\]
and call it \emph{the transport set associated to $\Gamma$}.
\end{defn}
Note that if $\Gamma$ is $\sigma$-compact, then so is $\mathbb{T}(\Gamma)$.
We also observe that if a measure $\lambda\in\mathcal{P}(\Omega^2)$ is concentrated on $\Gamma$, then
the measure $(P_t)_{\sharp}\lambda$ is concentrated on $\mathbb{T}(\Gamma)$ for any $t\in [0,1]$.
Now we recall the notion of density-regular points and regular points, introduced in \cite{CP11}.

\begin{defn}\label{def: densityregularity
Let $\pi\in\mathcal{P}(\Omega^2)$ be a probability measure concentrated on a $\sigma$-compact set $\Gamma\subset \Omega^2$.
Assume that the first marginal $\mu_1$ of $\pi$ is absolutely continuous with respect to $\mathcal{L}^n$ and let $f$ be the density.
A point $(x,y)\in \Gamma$ is said to be \emph{density-regular} if  for any $r>0$, there exists a point $\tilde{y}\in\Omega$ and a positive number $\tilde{r}>0$ such that
\begin{itemize}
\item $y\in B(\tilde{y},\tilde{r})\subset B(y,r)$,
\item The point $x$ is a Lebesgue point $f$ and $\tilde{f}$.
\item $f(x)<+\infty$, $\tilde{f}(x)>0$,
\end{itemize}
where $\tilde{f}$ is the density of $(\text{pr}_1)_{\sharp}(\pi\lfloor_{\Omega\times B(\tilde{y},\tilde{r})})$. We denote by $\mathop{\rm DR}\nolimits(\Gamma)$ the set of all density-regular points of $\Gamma$.
\end{defn}

\begin{defn}\label{def: densityregularity
Let $\Gamma\subset \Omega^2$ be a $\sigma$-compact set.
A point $(x,y)\in \Gamma$ is said to be \emph{regular} if for any $r>0$, $x$ is a Lebesgue point of $\Gamma^{-1}(B(y,r))$,
where the set $\Gamma^{-1}(B(y,r))$ is defined by
\[
\Gamma^{-1}(B(y,r)):=\left\{x \in \Omega \mathrel{}\middle|\mathrel{} (x,z)\in\Gamma\ \text{for some}\ z\in B(y,r)\right\}.
\]
\end{defn}

\begin{remark}
Any density-regular point is regular.
\end{remark}

\begin{lem}\label{lem: densityregularity}
Under the same settings as in Definition \ref{def: densityregularity},
$\pi$ is concentrated on $\mathop{\rm DR}\nolimits(\Gamma)$.
\end{lem}

\begin{prop}\label{prop: density regularity}
Let $\mu_1, \mu_2\in\mathcal{P}_c(\Omega)$ and let $\pi\in\mathcal{O}_2(\mu_1,\mu_2)$.
Assume that the measure $\mu_1$ is absolutely continuous with respect to $\mathcal{L}^n$ and let $f$ be the density of
$\mu_1$.
Then, for any point $(x,y)\in \mathop{\rm DR}\nolimits(\Gamma)$ with $x\neq y$ and positive real number $r>0$,
\[
\liminf_{\delta\to0}\frac{\mathcal{L}^n\left(\mathbb{T}\left(\Gamma\cap(B(x,{\delta/4})\times B(y,r))\right)\cap B(x,{\delta})\right)}{\mathcal{L}^n\left(B(x,{\delta})\right)}>0.
\]
\end{prop}

\begin{proof}
Define a Borel set $F\subset \Omega$ by
\[
F:=\left\{z\in \Omega \mathrel{}\middle|\mathrel{} {\tilde{f}(x)}/{2}\leq \tilde{f}(z)\right\} \cap
\left\{z\in \Omega \mathrel{}\middle|\mathrel{} f(z)\leq f(x)+1\right\}.
\]
Since $x$ is a Lebesgue point of $f$ and $\tilde{f}$, we have
\[
\lim_{s\to 0}\frac{\mathcal{L}^n\left(F\cap B(x,s)\right)}{\mathcal{L}^n\left(B(x,s)\right)}=1.
\]
Fix a positive real number $\delta>0$ so that $\delta/2(\lvert x-y\rvert + r)<1$ and
\begin{equation}\label{ineq: lebesguediffthm}
\frac{\mathcal{L}^n\left(F\cap B(x,s)\right)}{\mathcal{L}^n\left(B(x,s)\right)}\geq \frac12
\end{equation}
hold for any $0<s<\delta$.
By putting $F_{\delta}:=B(x,{\delta/4})\cap F$, we have
\begin{equation}\label{ineq: diff}
\frac{\tilde{f}(x)}{4}\mathcal{L}^n\left(B(x,{\delta/4})\right)
\leq \frac{\tilde{f}(x)}{2}\mathcal{L}^n(F_{\delta}).
\end{equation}
Put $A_{\delta}:=F_{\delta}\times B(\tilde{y},{\tilde{r}})$ and let $f_{\delta}$ be the density of $(\text{pr}_1)_{\sharp}\left(\pi\lfloor_{A_{\delta}}\right)$. By the definition of $\tilde{f}$ and $F_{\delta}$, we have
\begin{align*}
f_{\delta}(z) &= \limsup_{s\to 0}\frac{\pi\left((B(z,s)\cap F_{\delta})\times B(\tilde{y},{\tilde{r}})\right)}{\mathcal{L}^n\left(B(z,s)\right)} \\
&= \limsup_{s\to 0}\frac{1}{\mathcal{L}^n\left(B(z,s)\right)}\int_{B(z,s)\cap F_{\delta}} \tilde{f}\, d\mathcal{L}^n \\
&\geq \frac{\tilde{f}(x)}{2}\limsup_{s\to 0}\frac{\mathcal{L}^n\left(B(z,s)\cap F_{\delta}\right)}{\mathcal{L}^n\left(B(z,s)\right)}.
\end{align*}
Lebesgue's differentiation theorem tells
\[
f_{\delta}(z)\geq \frac{\tilde{f}(x)}{2}
\]
for $\mathcal{L}^n$-almost all $z\in F_{\delta}$.
Integrating this on $F_{\delta}$, we obtain
\begin{equation}\label{ineq: integrate}
\frac{\tilde{f}(x)}{2}\mathcal{L}^n(F_{\delta})\leq \int_{F_{\delta}}f_{\delta}(z)\, dz =(\text{pr}_1)_{\sharp}(\pi\lfloor_{A_{\delta}})(F_{\delta}).
\end{equation}
We also observe that $P_t(z,w)\in B(x,\delta)$ for any $0<t<\delta/(2(\lvert x-y\rvert+r))$, $z\in B(x,\delta/4)$ and $w\in B(y,r)$, which yields
\begin{equation}\label{ineq: ball}
(\text{pr}_1)_{\sharp}(\pi\lfloor_{A_{\delta}})(F_{\delta}) \leq (P_t)_{\sharp}(\pi\lfloor_{A_{\delta}})(B(x,{\delta})).
\end{equation}
Combining (\ref{ineq: diff}), (\ref{ineq: integrate}) and (\ref{ineq: ball}), we reach
\begin{equation}\label{ineq: interpolate}
\frac{\tilde{f}(x)}{4}\mathcal{L}^n\left(B(x,{\delta/4})\right)
\leq (P_t)_{\sharp}(\pi\lfloor_{A_{\delta}})(B(x,{\delta})).
\end{equation}

Since $\pi\in\mathcal{O}_2(\mu_1,\mu_2)$, there exist a sequence $\{\varepsilon_k\}_{k=1}^{\infty}$ of positive real numbers converging to $0$ and a sequence $\{\pi_{k}\in D_{\varepsilon_k}\}_{k=1}^{\infty}$ of  Borel probability measures converging to $\pi$ weakly.
From Lemma \ref{lem: densitybound} (iii), for any natural number $k$, we have
\begin{align*}
\|(P_t)_{\sharp}\left(\pi_{k}\lfloor_{F_{\delta}\times \Omega}\right)\|_{L^{\infty}} &\leq (1-t)^{-n}\|(\text{pr}_1)_{\sharp}\left(\pi_k\lfloor_{F_{\delta}\times \Omega}\right)\|_{L^{\infty}} \\
&\leq (1-t)^{-n}\|f\lfloor_{F_{\delta}}\|_{L^{\infty}} \leq 2^{n}(f(x)+1).
\end{align*}
By Lemma \ref{lem: restriction}, we know that the sequence $\{\pi_{k}\lfloor_{F_{\delta}\times \Omega}\}_{k=1}^{\infty}$
 converges to $\pi\lfloor_{F_{\delta}\times \Omega}$ weakly. Since the $L^{\infty}$-bound is stable under the weak convergence, we obtain
\[
\|(P_t)_{\sharp}\left(\pi\lfloor_{A_{\delta}}\right)\|_{L^{\infty}}
\leq \|(P_t)_{\sharp}(\pi\lfloor_{F_{\delta}\times \Omega})\|_{L^{\infty}}
\leq 2^n(f(x)+1).
\]
The fact that
$(P_t)_{\sharp}(\pi\lfloor_{A_{\delta}})$ is concentrated on
$\mathbb{T}\left(\Gamma \cap (B(x,{\delta/4})\times B(y,r))\right)$ implies that
\begin{align*}
&(P_t)_{\sharp}(\pi\lfloor_{A_{\delta}})(B(x,{\delta})) \\
&= (P_t)_{\sharp}(\pi\lfloor_{A_{\delta}})\left(\mathbb{T}\left(\Gamma \cap (B(x,{\delta/4})\times B(y,r))\right) \cap B(x,{\delta})\right) \\
 &\leq 2^n(f(x)+1)\mathcal{L}^n\left(\mathbb{T}\left(\Gamma \cap (B(x,{\delta/4})\times B(y,r))\right) \cap B(x,{\delta})\right).
\end{align*}
Combining this with (\ref{ineq: interpolate}), we arrive at
\[
\frac{\mathcal{L}^n\left(\mathbb{T}\left(\Gamma\cap(B(x,{\delta/4})\times B(y,r))\right)\cap B(x,{\delta})\right)}
{\mathcal{L}^n\left(B(x,{\delta/4})\right)}\geq \frac{\tilde{f}(x)}{2^{n+2}(f(x)+1)},
\]
which is equivalent to
\[
\frac{\mathcal{L}^n\left(\mathbb{T}\left(\Gamma\cap(B(x,{\delta/4})\times B(y,r))\right)\cap B(x,{\delta})\right)}
{\mathcal{L}^n\left(B(x,{\delta})\right)}\geq \frac{\tilde{f}(x)}{2^{3n+2}(f(x)+1)}.
\]
This completes the proof.
\end{proof}
The following theorem is the main result of this section and Theorem \ref{thm: main theorem} follows from this.
\begin{thm}\label{main theorem
Let $\mu_1,\mu_2\in \mathcal{P}_c(\Omega)$.
If the measure $\mu_1$ is absolutely continuous with respect to $\mathcal{L}^n$, then $\mathcal{O}_2(\mu_1,\mu_2)$ consists of a single element $\pi_0$ and the transport plan $\pi_0$ is induced by a map.
\end{thm}

\begin{proof}
Let $\pi\in\mathcal{O}_2(\mu_1,\mu_2)$. By virtue of Lemma \ref{lem: map} and Lemma \ref{lem: densityregularity}, it suffices to prove that the set $\mathop{\rm DR}\nolimits(\Gamma)$ is concentrated on the graph of a map.
Assume that there exist two points $(x_0,y_0), (x_0,y_1)\in \mathop{\rm DR}\nolimits(\Gamma)$ with $y_0\neq y_1$.
Then we have either $(y_1-y_0)\cdot (y_0-x_0)<0$ or $(y_0-y_1)\cdot (y_1-x_0)<0$. By symmetry, we may assume
$(y_1-y_0)\cdot (y_0-x_0)<0$. By continuity, for sufficiently small $r>0$, any points $x\in B(x_0,r)$, $y'\in B(y_0,r)$ and $y\in B(y_1,r)$ satisfy
\[
(y-y')\cdot (y'-x)<0.
\]
On the other hand, by the $\Gamma$-regularity of $(x_0,y_1)\in \Gamma$, we have
\[
\lim_{\delta\to 0}\frac{\mathcal{L}^n\left(\text{pr}_1(\Gamma \cap (\Omega \times B(y_1,r)))\cap B(x_0,{\delta})\right)}{\mathcal{L}^n(B(x_0,{\delta}))}=1
\]
and by the density-regularity of $(x_0,y_0)\in \mathop{\rm DR}\nolimits(\Gamma)$, we have
\[
\liminf_{\delta\to 0}\frac{\mathcal{L}^n\left(\mathbb{T}(\Gamma\cap (B(x_0,{\delta/4}))\times B(y_0,r))\cap B(x_0,{\delta})\right)}{\mathcal{L}^n(B(x_0,{\delta}))}>0.
\]
Thus for sufficiently small $0<\delta<r$, there exists a point $x\in B(x_0,{\delta})$ such that
\[
x\in \text{pr}_1(\Gamma \cap (\Omega \times B(y_1,r)))\cap \mathbb{T}(\Gamma\cap (B(x_0,{\delta/4}))\times B(y_0,r))).
\]
This implies the existence of points $y\in B(y_1,r)$, $x'\in B(x_0,{\delta/4})$ and $y'\in B(y_0,r)$ such that $(x,y), (x',y')\in \Gamma$ and $x\in [x',y']$. From the restricted monotonicity, we have
\[
(y-y')\cdot (x-x')\geq 0.
\]
Since $x$ lies in the line segment $[x',y']$, we have
\[
x-x'=a(y'-x')
\]
for some $a>0$. Thus we obtain
\[
(y-y')\cdot (x-x') = a(y-y')\cdot (y'-x)<0,
\]
which is a contradiction.
\end{proof}

\section{Hilbert geometries}
In this section, we recall the notion of Hilbert geometries, introduced by Hilbert \cite{Hilbert1895} himself.
We refer to \cite{Papa14} for details.

Let $\Omega\subset \mathbb{R}^n$ be a bounded convex open domain.
For any distinct two points $x,y\in \Omega$,
the \emph{Hilbert metric} $h_{\Omega}(x,y)$ from $x$ to $y$ is given by
\[
h_{\Omega}(x,y):=\frac{1}{2}\log\frac{\lvert y-x'\rvert\lvert x-y'\rvert}{\lvert x-x'\rvert\lvert y-y'\rvert},
\]
where the points $x'=x+s(y-x), s<0$ and $y'=x+t(y-x), t>0$ are the intersection of the affine straight line passing through $x$ and $y$ and the boundary $\partial \Omega$ (see Figure 1).
The function $h_{\Omega}$ defines a complete metric on $\Omega$ and the metric space $(\Omega, h_{\Omega})$ is called the \emph{Hilbert geometry} for $\Omega$.
The Hilbert geometry for a unit open disk coincides with the Cayley-Klein model of the hyperbolic geometry.
If the boundary $\partial \Omega$ is $C^2$-regular and strictly convex, then the Hilbert metric on $\Omega$ is induced by the smooth Finsler structure $F$, given by
\[
F(x,v)=\frac{\lvert v\rvert}{2}\left(\frac{1}{\lvert x-a\rvert}+\frac{1}{\lvert x-b\rvert}\right),
\]
where the points $a=x+sv, s<0,$ and $b=x+tv, t>0$, are the intersection of the affine straight line passing through $x$ with the direction $v$ (see Figure 1). In this case, the Finsler manifold $(\Omega, F)$ has the constant flag curvature $-1$ (cf. \cite{Okada83}).
\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{hilbertgeometry.png}
\caption{Hilbert metric}
\end{figure}
For any bounded convex open set $\Omega$, the Hilbert metric $h_{\Omega}$ satisfies (i) and (ii) of Assumption A in  Section 1.
To verify (iii), we use the following result, proved by Ohta \cite{Ohta13}.
\begin{thm}[\cite{Ohta13}]
Let $\Omega\subset \mathbb{R}^n$ be a bounded convex open domain and $h_{\Omega}$ the Hilbert metric.
Then the metric-measure space $(\Omega,h_{\Omega},\mathcal{L}^n)$ satisfies the curvature-dimension condition $\mathop{\rm CD}\nolimits(K,N)$ for
\[
K=-(n-1)-\frac{(n+1)^2}{N-n}\ \text{and}\ N>n.
\]
\end{thm}
The curvature-dimension condition $\mathop{\rm CD}\nolimits(K,N)$, which we do not explain in this paper, is a generalization of the condition that the weighted Ricci curvature is bounded below by $K$ and the dimension is bounded above by $N$.
If a metric-measure space $(X,d,\mathfrak{m})$ satisfies the curvature-dimension condition $\mathop{\rm CD}\nolimits(K,N)$ for $K\in \mathbb{R}$ and finite $N>1$, then the base measure $\mathfrak{m}$ satisfies the local doubling property.
In particular, the Hilbert metrics satisfy (iii) of Assumption A, so we obtain Corollary \ref{cor: Hilb geom}.

\bibliographystyle{plain}
\begin{bibdiv}
\begin{biblist}

\bib{AKP04}{article}{
      author={Ambrosio, L.},
      author={Kirchheim, B.},
      author={Pratelli, A.},
       title={Existence of optimal transport maps for crystalline norms},
        date={2004},
        ISSN={0012-7094},
     journal={Duke Math. J.},
      volume={125},
      number={2},
       pages={207\ndash 241},
         url={https://doi.org/10.1215/S0012-7094-04-12521-7},

}

\bib{AG13}{incollection}{
      author={Ambrosio, Luigi},
      author={Gigli, Nicola},
       title={A user's guide to optimal transport},
        date={2013},
   booktitle={Modelling and optimisation of flows on networks},
      series={Lecture Notes in Math.},
      volume={2062},
   publisher={Springer, Heidelberg},
       pages={1\ndash 155},
         url={https://doi.org/10.1007/978-3-642-32160-3_1},

}

\bib{AP03}{incollection}{
      author={Ambrosio, Luigi},
      author={Pratelli, Aldo},
       title={Existence and stability results in the {$L^1$} theory of optimal
  transportation},
        date={2003},
   booktitle={Optimal transportation and applications ({M}artina {F}ranca,
  2001)},
      series={Lecture Notes in Math.},
      volume={1813},
   publisher={Springer, Berlin},
       pages={123\ndash 160},
         url={https://doi.org/10.1007/978-3-540-44857-0_5},

}

\bib{BiCa11}{incollection}{
      author={Bianchini, Stefano},
      author={Cavalletti, Fabio},
       title={The {M}onge problem in geodesic spaces},
        date={2011},
   booktitle={Nonlinear conservation laws and applications},
      series={IMA Vol. Math. Appl.},
      volume={153},
   publisher={Springer, New York},
       pages={217\ndash 233},
         url={https://doi.org/10.1007/978-1-4419-9554-4_10},

}

\bib{BiCa13}{article}{
      author={Bianchini, Stefano},
      author={Cavalletti, Fabio},
       title={The {M}onge problem for distance cost in geodesic spaces},
        date={2013},
        ISSN={0010-3616},
     journal={Comm. Math. Phys.},
      volume={318},
      number={3},
       pages={615\ndash 673},
         url={https://doi.org/10.1007/s00220-013-1663-8},

}

\bib{Cava17}{incollection}{
      author={Cavalletti, Fabio},
       title={An overview of {$L^1$} optimal transportation on metric measure
  spaces},
        date={2017},
   booktitle={Measure theory in non-smooth spaces},
      series={Partial Differ. Equ. Meas. Theory},
   publisher={De Gruyter Open, Warsaw},
       pages={98\ndash 144},

}

\bib{CP11}{article}{
      author={Champion, Thierry},
      author={De~Pascale, Luigi},
       title={The {M}onge problem in {$\Bbb R^d$}},
        date={2011},
        ISSN={0012-7094},
     journal={Duke Math. J.},
      volume={157},
      number={3},
       pages={551\ndash 572},
         url={https://doi.org/10.1215/00127094-1272939},

}

\bib{Hilbert1895}{article}{
      author={Hilbert, David},
       title={Ueber die gerade linie als k{\"u}rzeste verbindung zweier
  punkte},
        date={1895Mar},
        ISSN={1432-1807},
     journal={Mathematische Annalen},
      volume={46},
      number={1},
       pages={91\ndash 96},
         url={https://doi.org/10.1007/BF02096204},
}

\bib{Kantorovich}{article}{
      author={Kantorovitch, L.},
       title={On the translocation of masses},
        date={1958},
        ISSN={0025-1909},
     journal={Management Sci.},
      volume={5},
       pages={1\ndash 4},
         url={https://doi.org/10.1287/mnsc.5.1.1},

}

\bib{Monge1781}{book}{
      author={Monge, Gaspard},
       title={M{\'e}moire sur la th{\'e}orie des d{\'e}blais et des remblais},
   publisher={De l'Imprimerie Royale},
        date={1781},
}

\bib{Ohta13}{article}{
      author={Ohta, Shin-ichi},
       title={Weighted {R}icci curvature estimates for {H}ilbert and {F}unk
  geometries},
        date={2013},
        ISSN={0030-8730},
     journal={Pacific J. Math.},
      volume={265},
      number={1},
       pages={185\ndash 197},
         url={https://doi.org/10.2140/pjm.2013.265.185},

}

\bib{Okada83}{article}{
      author={Okada, Tsutomu},
       title={On models of projectively flat {F}insler spaces of constant
  negative curvature},
        date={1983},
        ISSN={0040-3504},
     journal={Tensor (N.S.)},
      volume={40},
      number={2},
       pages={117\ndash 124},

}

\bib{Papa14}{book}{
      author={Papadopoulos, Athanase},
       title={Metric spaces, convexity and non-positive curvature},
     edition={Second},
      series={IRMA Lectures in Mathematics and Theoretical Physics},
   publisher={European Mathematical Society (EMS), Z\"{u}rich},
        date={2014},
      volume={6},
        ISBN={978-3-03719-132-3},
         url={https://doi.org/10.4171/132},

}

\bib{San09}{article}{
      author={Santambrogio, Filippo},
       title={Absolute continuity and summability of transport densities:
  simpler proofs and new estimates},
        date={2009},
        ISSN={0944-2669},
     journal={Calc. Var. Partial Differential Equations},
      volume={36},
      number={3},
       pages={343\ndash 354},
         url={https://doi.org/10.1007/s00526-009-0231-8},

}

\bib{Vil09}{book}{
      author={Villani, C\'{e}dric},
       title={Optimal transport},
      series={Grundlehren der Mathematischen Wissenschaften [Fundamental
  Principles of Mathematical Sciences]},
   publisher={Springer-Verlag, Berlin},
        date={2009},
      volume={338},
        ISBN={978-3-540-71049-3},
         url={https://doi.org/10.1007/978-3-540-71050-9},
        note={Old and new},

}

\bib{Wil70}{book}{
      author={Willard, Stephen},
       title={General topology},
   publisher={Addison-Wesley Publishing Co., Reading, Mass.-London-Don Mills,
  Ont.},
        date={1970},

}

\end{biblist}
\end{bibdiv}

\end{document}


\section{Introduction}
\label{sec:introduction}
Magnetic Resonance Imaging (MRI) is an important tool for examination and diagnosis of a wide array of pathologies under diverse circumstances. Its analysis has several applications, such as the study of Alzheimer's disease evolution~\cite{Wang16-JAD}, multiple sclerosis lesions~\cite{Garcia13-MIA}, and bipolar disorder~\cite{Hibar18-MP}, to name a few examples.

However, manual analysis is often subjective~\cite{Morel16-NR} and time-intensive. Thus, fast and automatic methods of MRI analysis are required. As several studies are based on the segmentation of different tissues or structures in the body, and the amount of patient data available grows significantly, automatic or semi-automatic segmentation of regions of interest becomes essential.

While several methods for the analysis of adult MRI have been developed over the years, relatively few have been developed for pediatric MRI~\cite{Devi15-CBM}. Pediatric images present particular challenges~\cite{Devi15-CBM,Neubauer13-PR}, such as the inversion of the contrasts in white and gray matter, patient motion artifacts~\cite{Devi15-CBM} and significant corporal variation among individuals with small age differences~\cite{Huppi98-AN}.

Additionally, most segmentation methods focus on healthy patients, or on a single anomaly (such as a specific type of tumor)~\cite{harish2017-Extensive}. There is a demand for methods that are capable of dealing with a wide range of abnormalities -- extraneous elements such as tumors or missing elements such as organs removed by surgery -- in the corporal structure without significant retraining or adjustment.

When dealing with structured scenes, that is, scenes where the objects of interest obey some sort of predictable spatial structure (such as the interior of the body), the usage of structural information to aid segmentation and analysis has been shown to allow for significant increases in accuracy~\cite{Colliot06-PR,Fouquier12-CVIU,Morel16-EMBC,Nempont13-IS}.

The goal of this project is to develop a structural (graph-based) method for the analysis of pediatric images. This method will be applied to real MRI data in order to solve active problems in the area.

\section{Methods and Approach}
\label{sec:methods}
The project will involve the design of a suitable graph-based structural model for image segmentation and body parts recognition. In order to accurately model and represent the knowledge found in MRI data, proper understanding of the structural components of the images is required. As shown by \cite{Colliot06-PR,Fouquier12-CVIU,Morel16-EMBC,Nempont13-IS}, structural knowledge is of primary importance for the interpretation of images where intensity and shape may not be sufficient data for an adequate analysis, such as abdominal and brain MRIs. To mitigate this deficiency, we propose the usage of graphs (and, potentially, hypergraphs) for the representation of anatomical structures, where each vertex represents a region or structure of interest and the edges represent relations between these (such as relative spatial positioning). Our challenge will be the application of this knowledge in a way that guides and enhances the image interpretation process.  Figure \ref{fig:pipe} presents the proposed pipeline to be followed in order to segment and recognize structures of interest in the MRI sequence.

\begin{figure}[h!]
\centering
\includegraphics[width = 1\linewidth]{new_pipeline}
\caption{Pipeline showing the steps of the segmentation task.}
\label{fig:pipe}
\end{figure}

\subsection{The Statistical-Relational Graph}
\label{sec:srg}
The Statistical-Relational Graph (SRG)~\cite{Graciano12-THESIS} is an hybrid structured scene representation technique based on Attributed Relational Graphs, where vertices of a graph represent primitive elements of the scene and edges represent their structural relationship. Both vertices and edges contain attributes which describe the element or relationship they represent.

For recognition of structured scenes, we first build a graph which represents the template for a given scene, with a set of primitive and relational attributes, called the \emph{model graph}. Object recognition is thus a matter of building an equivalent or comparable graph from the observed input, the \emph{observation graph}, and attempting to match both graphs.

Structural-relational graphs have been successfully applied to object tracking in digital image sequences~\cite{Graciano07-SIBGRAPI,Morimitsu15-THESIS}, automatic image segmentation~\cite{Graciano12-THESIS} and interactive image segmentation~\cite{Noma08-ARXIV, Noma12-PR}.

In their thesis, Graciano~\cite{Graciano12-THESIS} demonstrates this approach's robustness in the segmentation of abdominal organs in normal, adult patient MRIs, utilizing classic probabilistic pattern recognition techniques to describe organ tissues while simultaneously considering the spatial relationship between distinct organs in order to improve segmentation quality. We believe this method may be further developed for the automatic segmentation and recognition of regions of interest in pediatric MRIs.

Primitive elements, in the context of MRI analysis, may be organs, veins, tumors, or similar regions of interest such as hepatic divisions. Discriminative attributes for these elements may include mean intensity, total volume, and position in the body. Structural relationships describes the connections between distinct primitive elements, such as relative position and distance, relative intensity and relative size.

\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{srg_example}
\caption{Schema of a simple, example SRG. Vertex attributes are between square brackets. Edge attributes are between curly braces.}
\label{fig:srg_example}
\end{figure}

Figure \ref{fig:srg_example} presents an example of a SRG. The vertices represent some abdominal structures of interest: Segment I (SI), Segment V (SV) and the Vena Cava (VC). The vertex attributes are the `centroid' and the `mean intensity', calculated from the MRI for each of the structures (vertices). The edge attribute is `distance', which represents the vector between the centroid of the two structures in the edge's vertices.

\subsection{Model Description}
\label{sec:model_description}

In order to perform our tasks of segmentation and recognition, we must first build a model capable of extracting and properly representing the structures and objects of interest. This model must be able to discriminate precisely between structures, using both statistic attributes from the structures and features derives from the relationships between them.

\subsubsection{Conceptual modeling}
\label{ssec:conceptual}

In order to properly segment and recognize the structures of interest, the following attributes of each structure must be considered: \textbf{i)} Image intensity values; \textbf{ii)} volume; and \textbf{iii)} absolute position in the body.

Additionally, the relationships between a structure and all the others must also be taken into consideration, with the following attributes being most relevant: \textbf{i)} relative distance, position and orientation; \textbf{ii)} differences in volume; and \textbf{iii)} contrast in intensity values.

We believe the combination of these statistical and relational attributes represent an accurate approximation of the discriminative factors required for proper segmentation and recognition of the relevant structures. Previous works have shown successful results produced by such attributes~\cite{Graciano12-THESIS, Colliot06-PR,Graciano14-ICPR,Morimitsu16-CVIU}.

\subsubsection{Computational modeling}
\label{ssec:computational_model}

Considering the conceptual model proposed in Section~\ref{ssec:conceptual}, we then define a single data structure, known as Statistical-Relational Graph (SRG), which will be used to represent the conceptual model.

Formally, the SRG is defined as a quadruple $G = (V, E, A_V, A_E)$ where:
\begin{itemize}
\item $V$ is the set of vertices in $G$. Each vertex represents a single part of the object of interest;
\item $E$ is the set of edges in $G$. Each edge represents a relationship between a pair of parts of the object;
\item $A_V$ is the set of vertex attributes and their respective probability distributions. Each element in $A_V$ maps to a single vertex in $V$ and contains that vertex's attributes;
\item $A_E$ is the set of edge attributes and their respective probability distributions. Each element in $A_E$ maps to a single edge in $E$ and contains that edge's attributes.
\end{itemize}

Three vertex attributes are computed for each vertex. Considering the structure they represent, these attributes are \textbf{i)} the average of the intensity values in the structure; \textbf{ii)} the total volume of the structure; and \textbf{iii)} the position of the centroid of the structure.

The volume attribute, in particular, may hinder initial matching efforts due to the difference between the volumes of the model vertices and the observation vertices, as the latter are much less voluminous than the former. However, this project aims to explore a way to utilize volume and, potentially, orientation information to aid the matching of a SRG.

Three edge attributes are computed for each edge. Considering the structures represented by the connected edges, these attributes are \textbf{i)} the vector between the centroids of each structure; \textbf{ii)} the proportional difference in volume between the structures; and \textbf{iii)} the contrast (that is, difference) between the average intensities of each structure.

\subsubsection{Graph creation}
\label{sssec:graph_creation}
From an annotated volume, an SRG is created as follows: for each distinct label, a vertex is created. The corresponding attributes are then computed for each label.
Once all vertices are created, the graph is then fully connected, and edge attributes are computed.

A model graph is created from manually annotated images. One of this project's goals is to find the best possible model graph from a set of such images. The sets of vertices and edges ($(V, E)$) will be defined \textit{a priori}. For example, $V$ will be provided from an anatomical atlas and $E$ will have a specified topology (e.g. fully connected). Then, learning algorithms will be applied in order to learn  the attributes' distributions in the statistical model. Initially, we'll consider a few multivariate distributions, such as Gaussian, for both parameters.

A super-observation graph is created from an automatic super-segmentation algorithm such as watershed. From the super-observation graph, regions may be joined to form an observation graph of the same size as the model graph.

\subsubsection{Graph matching}
\label{sssec:graph_matching}
In order to compare two SRGs, a matching solution must be generated and evaluated. A solution $S$, which matches a super observation graph $G_{super}$ with $|V_{super}| = n_{super}$ vertices to a model graph $G_M$ with $|V_M| = n$ vertices, is defined as
\begin{equation}
	S = [s_1, s_2, ..., s_{n_{super}}], s_i \in \{1,n\}
\end{equation}

From $S$, all regions with the same prediction are joined and the observation graph $G_{obs}$, with $|V_{obs}| = n$, is created. To evaluate the quality of $S$, a cost function must be established. The cost of a solution $S$ is defined as
\begin{equation}
	C(S) = \alpha \frac{1}{n} \sum_{j=1}^{n}c_V(s_j) + (1-\alpha)\frac{1}{{n}^2} \sum_{j=1}^{n}\sum_{k=1}^{n}c_E(s_j,s_k)
\end{equation}
where $\alpha \in [0,1]$ is a the weight of the verticial cost; the costs between vertices $j$ in $G_{obs}$ and $s_j$ in $G_M$, $c_V(s_j)$, is defined as
\begin{equation}
	c_V(s_j) = \sum_{a\in A_V}{\alpha_a d_a({s_j}_a)}
\end{equation}
and the cost between the edges $\{j,k\}$ in $G_{obs}$ and $\{s_j,s_k\}$ in $G_M$, $c_E(s_j,s_k)$, is defined as
\begin{equation}
	c_E(s_j,s_k) = \sum_{a\in A_E}{\alpha_a d_a(\{s_j\}_a, \{s_k\}_a)}
\end{equation}
with $d_a$ being the distance function corresponding to each attribute $a$ and $\alpha_a$ being a customizable parameter, used for tuning the weight of the distance between each attribute $a$ in the final cost. Note that $\sum_{a\in A_E}{\alpha_a} = 1$.

Given a super-observation SRG $G_{super}$, an initial solution may be built using a simple greedy algorithm,
where each super-observation vertex $j \in V_{super}$ is matched to the model vertex $i \in V_M$ with the lowest cost, that is, $s_j = argmin_{i=0}^{n}(c_V(j,i)), j=0,1,...,n_{super}$.

\section{Experiments}
\label{sec:experiments}
\subsection{Data Acquisition}
The main source of data for this project is the ICr-HC-FMUSP. All acquisitions were performed using ICr's Philips Ingenia (1.5T) machine. Abdominal T2-weighted MRI of nine patients, ranging from four to eighteen years old, were used for experiments in liver segmentation; among these patients, six were affected with hepatic iron buildup and three were healthy at the time of acquisition. From the same source, there is also brain T1, T2 and Diffusion MRI of three healthy patients. Their ages were 7, 8 months and 13 years at the time of acquisition. Further patients will be introduced to the dataset, with ages ranging from 6 months old to 13 years old, which may improve the method's efficiency.

\subsection{Exploratory SRG Experiments}
Exploratory experiments were conducted utilizing only \emph{mean intensity}, \emph{centroid} and \emph{distance} features and data from a single patient from the abdominal MRI dataset. The goal of these experiments was to determine the usefulness of the method for solving the problem in the given data, and spot and solve any potential problems that may arise in a more complete implementation.

The model graph was created from the annotated patient's data. The observation graph was built from an watershed segmentation of the morphological gradient of the full volum
; the regions produced by this segmentation were used as the vertices of the observation graph. We then matched both graphs using the simple greedy algorithm described in Section~\ref{sssec:graph_creation}. We explored several parameters for the $\alpha_a$ tuning parameters.

\subsection{Results and Discussion}
The results of the exploration of different vertex attribute weights, and their effects both on establishing an initial, greedy solution and computing the total cost of the solution, are shown on Table~\ref{tab:preliminary_results}. Selected results are displayed in Figure~\ref{fig:preliminary_results}.

\begin{table}
	\centering
	\caption{Results of the exploratory experiments on SRG.}
	\label{tab:preliminary_results}
    \begin{tabular}{c|c|c}
    	\textbf{Centroid $\alpha$} & \textbf{Intensity $\alpha$} & \textbf{Cost}\\
        \hline
        0 & 1 & $39.659.347$ \\
        0.001 & 0.999 & $29.912.702$ \\
        0.005 & 0.995 & $23.693.941$ \\
        0.01 & 0.99 & $23.064.714$ \\
        0.02 & 0.98 & $22.803.427$ \\
        0.1 & 0.9 & $22.720.578$ \\
        0.2 & 0.8 & $22.718.200$ \\
        0.5 & 0.5 & $22.718.200$ \\
        1 & 0 & $22.718.200$ \\
    \end{tabular}
\end{table}

\begin{figure}
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.45\linewidth]{pred01.png}&\includegraphics[width=0.45\linewidth]{pred001999.png}\\
        (a)&(b)\\
		\includegraphics[width=0.45\linewidth]{pred0199.png}&\includegraphics[width=0.45\linewidth]{pred19.png}\\
        (c)&(d)\\
		\includegraphics[width=0.45\linewidth]{pred28.png}&\includegraphics[width=0.45\linewidth]{truth11.png}\\
        (e)&(f)
	\end{tabular}
    \caption{Selected results from the exploratory experiments on SRG. All images are from the same slice. Weights given are for the \emph{centroid} and the \emph{mean intensity}, respectively. All results for \emph{centroid} weights above 0.2 are similar.}
    \label{fig:preliminary_results}
\end{figure}

We can notice that the \emph{centroid} attribute dominates the \emph{mean intensity}; indeed, for all weights of the \emph{intensity} attribute below $0.8$, including $0.0$, there are no differences in cost. While this may point to a significant advantage in using the \emph{centroid}, this may also be the result of a deficiency of more discriminant attributes (such as volume), of improper segmentation in the watershed step, and of improper normalization of the attributes. These results have paved the way for the next set of experiments, which will explore different sets of segmentation algorithms, region joining in the observation graph, and attribute normalization techniques.

\section{Conclusion}
\label{sec:conclusion}
We proposed a method for automatic segmentation and recognition of structures of interest in pediatric MRI. As stated before, this is an active and interesting problem with very clear real-world applications. The success of this project will both help pave the way for further applications using SRGs in medical imaging, and provide the scientific and medical community with a useful product, capable of aiding the radiologist in diagnostic and analysis of liver MRIs.

The proposed framework is already under implementation, and relevant data has been acquired and prepared for use. Previous work on similar methods show that this approach has potential for fast and accurate segmentation and recognition of structures, as well as high extensibility and capability for improvement and modification, such as transposing the proposed method to a different domain.

However, as the preliminary results of the exploratory experiments show, there are significant avenues for improvement that must be explored. Improving the quality of the observation graph, through the use of more discriminant attributes, better algorithms for building solutions, and improved super-segmentation approaches, is of utmost importance.

\section*{Acknowledgment}
The authors would like to thank the Funda\c{c}\~ao de Apoio \`a Pesquisa do Estado de S\~ao Paulo (FAPESP), grants \#2017/09465-7, \#2018/07386-5, \#2017/50236-1 and \#2015/22308-2; and also the Coordena\c{c}\~ao de Aperfei\c{c}oamento de Pessoal de N\'ivel Superior (CAPES), grant \#1767508.

\bibliographystyle{IEEEtran}