$\xi_{\omega}(x_{t,j})=\int_{x}^{x_{t,j}}p^{\ast}\omega=\int_{x}^{\mathsf{g}_{t_{\ast}}(x)}p^{\ast}\omega+\int_{\mathsf{g}_{t_{\ast}}(x)}^{\mathsf{g}_{t_{\ast}-t}(x_{t,j})}p^{\ast}\omega+\int_{\mathsf{g}_{t_{\ast}-t}(x_{t,j})}^{x_{t,j}}p^{\ast}\omega.$

The first term above is the Frobenius function (the geodesic winding cycle) as
in (12); the points $\mathsf{g}_{t_{\ast}}(x)$ and
$\mathsf{g}_{t_{\ast}-t}(x_{t,j})$ in the second integral both lie on a
horocycle orbit segment of length 1, and the last integral is over a segment
of length $\mathcal{O}(\log t)$. Thus, there exists a vector
$\mathbf{e}\in{\mathbb{R}}^{d}$ of norm $\|\mathbf{e}\|\leq{C_{\sharp}}\log t$
so that

$\left\|\xi_{\omega}(x_{t,j})-\int_{x}^{\mathsf{g}_{t_{\ast}}(x)}p^{\ast}\omega\right\|\leq{C_{\sharp}}|\omega\cdot\mathbf{e}|.$

Here, we have used the identification of the vector space $\mathcal{H}$ with
${\mathbb{R}}^{d}$ with the aid of the basis of 1-forms $\omega_{k}$ for
$k=1,\dots,d$ we fixed in Section 4. Let us define
$F_{\ast}=F_{\ast}(x,t_{\ast})$ to be the vector in ${\mathbb{R}}^{d}$ whose
components are
$F_{t_{\ast},\omega_{k}}(x)=\int_{x}^{\mathsf{g}_{t_{\ast}}(x)}p^{\ast}\omega_{k}$.
Then, since for all $\omega\in B(0,\delta t^{-1/4})$ we can bound

$|e^{-2\pi\imath\xi_{\omega}(x_{t,j})}-e^{-2\pi\imath\omega\cdot
F_{\ast}}|\leq{C_{\sharp}}|\omega\cdot\mathbf{e}|\leq{C_{\sharp}}\|\omega\|\log
t,$

we deduce

$\left\lvert\int_{B(0,\delta
t^{-1/4})}e^{-2\pi\imath\xi_{\omega}(x_{t,j})-2\pi^{2}t\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega-\int_{B(0,\delta
t^{-1/4})}e^{-2\pi\imath\omega\cdot
F_{\ast}-2\pi^{2}t\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega\right\rvert\\\
\leq{C_{\sharp}}\log t\int_{B(0,\delta
t^{-1/4})}e^{-t\|\omega\|^{2}}\|\omega\|\mathop{}\\!\mathrm{d}\omega\leq{C_{\sharp}}(\log
t)\,t^{-\frac{d+1}{2}}.$

Thus, from (52) and recalling (48), we obtained

(54)
$\left\lvert\int_{0}^{T}f\circ\mathsf{h}_{s}(x)\mathop{}\\!\mathrm{d}s-e^{h_{\operatorname{top}}t}\left(\sum_{j\in\mathcal{J}}A_{j}(0)\right)\int_{B(0,\delta
t^{-1/4})}e^{-2\pi\imath\omega\cdot
F_{\ast}-2\pi^{2}t\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega\right\rvert\\\
\leq{C_{\sharp}}C(f,x)\frac{T\cdot\log\log T}{(\log T)^{\frac{d+1}{2}}}.$

Now, by Lemma 7.5, we have

$\left\lvert\sum_{j\in\mathcal{J}}A_{j}(0)-\mu(f)\tau(T,t,x)\right\rvert=|\mu(f)|\left\lvert\int_{0}^{\tau(T,t,x)}\psi(s)\mathop{}\\!\mathrm{d}s-\tau(T,t,x)\right\rvert\leq{C_{\sharp}}C(f).$

Furthermore, since $|t_{\ast}-t|\leq{C_{\sharp}}\log t\leq{C_{\sharp}}\log\log
T$, we can also replace $t$ with $t_{\ast}$ in (54). Using also Proposition
7.3 to bound
$|e^{h_{\operatorname{top}}t}\tau(T,t,x)-T|\leq{C_{\sharp}}Te^{-\delta_{*}t}\leq{C_{\sharp}}T^{1-\delta_{*}}$,
we deduce

(55)
$\left\lvert\int_{0}^{T}f\circ\mathsf{h}_{s}(x)\mathop{}\\!\mathrm{d}s-T\mu(f)\int_{B(0,\delta
t_{\ast}^{-1/4})}e^{-2\pi\imath\omega\cdot
F_{\ast}-2\pi^{2}t_{\ast}\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega\right\rvert\leq{C_{\sharp}}C(f,x)\frac{T\cdot\log\log
T}{(\log T)^{\frac{d+1}{2}}}.$

We can now finish the proof by computing the second integral in the left hand
side above. Let $H$ be the symmetric positive definite square root of the
matrix $\Sigma$, then

$\int_{B(0,\delta t_{\ast}^{-1/4})}e^{-2\pi\imath\omega\cdot
F_{\ast}-2\pi^{2}t_{\ast}\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega=t_{\ast}^{-\frac{d}{2}}\int_{B(0,\delta
t_{\ast}^{1/4})}e^{-2\pi\imath\omega\cdot\frac{F_{\ast}}{\sqrt{t_{\ast}}}-2\pi^{2}(H\omega)\cdot(H\omega)}\mathop{}\\!\mathrm{d}\omega\\\
=t_{\ast}^{-\frac{d}{2}}\int_{{\mathbb{R}}^{d}}e^{-2\pi\imath\omega\cdot\frac{F_{\ast}}{\sqrt{t_{\ast}}}-2\pi^{2}(H\omega)\cdot(H\omega)}\mathop{}\\!\mathrm{d}\omega+O(t_{\ast}^{-1})=\frac{1}{t_{\ast}^{\frac{d}{2}}\sqrt{\det\Sigma}}\int_{{\mathbb{R}}^{d}}e^{-2\pi\imath\mathbf{y}\cdot\frac{H^{-1}F_{\ast}}{\sqrt{t_{\ast}}}-2\pi^{2}\mathbf{y}\cdot\mathbf{y}}\mathop{}\\!\mathrm{d}\mathbf{y}+O(t_{\ast}^{-1})\\\
=\frac{1}{(2\pi
t_{\ast})^{\frac{d}{2}}\sqrt{\det\Sigma}}e^{-\frac{1}{2}\left\|\frac{H^{-1}F_{\ast}}{\sqrt{t_{\ast}}}\right\|^{2}}+O(t_{\ast}^{-1})=\frac{1}{(2\pi
t_{\ast})^{\frac{d}{2}}\sqrt{\det\Sigma}}e^{-\frac{1}{2}\left\|\frac{F_{\ast}}{\sqrt{t_{\ast}}}\right\|_{\Sigma}^{2}}+O(t_{\ast}^{-1}).$

where, we recall, the norm $\|\cdot\|_{\Sigma}$ on ${\mathbb{R}}^{d}$ is
defined as $\|\mathbf{x}\|_{\Sigma}^{2}=\mathbf{x}\cdot\Sigma^{-1}\mathbf{x}$.

Using again Proposition 7.3 to bound
$|e^{h_{\operatorname{top}}t_{\ast}}-T|\leq{C_{\sharp}}e^{(h_{\operatorname{top}}-\delta_{*})t_{\ast}}\leq{C_{\sharp}}T^{1-\delta_{*}^{\prime}}$
for some $\delta_{*}^{\prime}>0$, from (55) we conclude

$\left\lvert\int_{0}^{T}f\circ\mathsf{h}_{s}(x)\mathop{}\\!\mathrm{d}s-\frac{h_{\operatorname{top}}^{\frac{d}{2}}T}{(2\pi\log
T)^{\frac{d}{2}}\sqrt{\det\Sigma}}\left(\int_{{\widetilde{M}}}f\mathop{}\\!\mathrm{d}\mu\right)e^{-\frac{1}{2}\left\|\frac{F_{\ast}(p(x),t_{\ast})}{\sqrt{t_{\ast}}}\right\|_{\Sigma}^{2}}\right\rvert\leq{C_{\sharp}}C(f,x)\frac{T\cdot\log\log
T}{(\log T)^{\frac{d+1}{2}}}.$

The vector $F_{\ast}=F_{\ast}(p(x),t_{\ast})$ has components
$F_{t_{\ast},\omega_{k}}(x)$, where $\omega_{k}$, for $k=1,\dots,d$ form the
basis we fixed in Section 4. Thus, the Central Limit Theorem for the geodesic
flow implies

$\frac{F_{\ast}(\cdot,t_{\ast})}{\sqrt{t_{\ast}}}\to\mathcal{N}(0,\Sigma)\qquad\text{in
distribution.}$

The proof of A is therefore complete.

###  7.3. Proof of B

Let us fix $x\in{\widetilde{M}}$ and $\sigma>0$, and recall that we denoted by
$\gamma_{x,\sigma}$ the horocycle segment $s\mapsto\mathsf{h}_{s}(x)$ for
$s\in[0,\sigma]$. Let $\eta$ be a 1-form of class $\mathscr{C}^{2}$ on
${\widetilde{M}}$ with compact support, and let
$f:=f_{\eta}\in\mathscr{C}^{2}_{c}({\widetilde{M}})$ be given by
$f(x)=\langle\eta,U\rangle_{x}$. We also assume that we fixed a fundamental
domain ${\mathcal{F}}\subset{\widetilde{M}}$ containing $x$ and we take
$x=x_{0}$ to be the base point in ${\mathcal{F}}$ as in Section 4.1.

The proof is analogous to the one presented in Section 7.2. By (6), we have

$\int_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}}\eta=\int_{0}^{\sigma}\eta_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}(s)}\left(D\mathsf{g}_{-t}U_{\gamma_{x,\sigma}(s)}\right)\mathop{}\\!\mathrm{d}s=\int_{0}^{\sigma}f\circ\mathsf{g}_{-t}\circ\mathsf{h}_{s}(x)\cdot
J_{-t}(\mathsf{h}_{s}(x))\mathop{}\\!\mathrm{d}s.$

Let $\delta_{*}>0$ be given by Proposition 7.3, and define
$r=d(\min\\{h_{\operatorname{top}},\delta_{*}\\})^{-1}\log t$; in particular,
note that $e^{-h_{\operatorname{top}}r}\leq t^{-d}$. Calling
$x_{-r}=\mathsf{g}_{-r}(x)$, we rewrite the equation above as

$\begin{split}\int_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}}\eta&=\int_{0}^{\sigma}{\mathcal{L}}_{t}f\circ\mathsf{h}_{s}(x)\mathop{}\\!\mathrm{d}s=\int_{0}^{\sigma}{\mathcal{L}}_{t-r}(f)\circ\mathsf{g}_{-r}\circ\mathsf{h}_{s}(x)\cdot
J_{-r}(\mathsf{h}_{s}(x_{-r}))\mathop{}\\!\mathrm{d}s\\\
&=\int_{0}^{\tau(\sigma,r,x)}{\mathcal{L}}_{t-r}(f)\circ\mathsf{h}_{s}(x)\mathop{}\\!\mathrm{d}s,\end{split}$

We note that, by Lemma 3.6, $\tau=\tau(\sigma,-r,x)$ satisfies

$C_{\tau}^{-1}\sigma\,e^{h_{\operatorname{top}}r}\leq\tau\leq
C_{\tau}\sigma\,e^{h_{\operatorname{top}}r}.$

As in the proof of Lemma 7.2, we define a function
$\psi\colon[0,\tau]\to[0,1]$ of class $\mathscr{C}^{1+\operatorname{Lip}}$
with $\|\psi^{\prime}\|_{\operatorname{Lip}}\leq 4$ and equal to 1 in
$[1,\tau-1]$, so that

$\left\lvert
e^{-h_{\operatorname{top}}t}\int_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}}\eta-e^{-h_{\operatorname{top}}r}\int_{0}^{\tau}\hat{{\mathcal{L}}}_{t-r}(f)\circ\mathsf{h}_{s}(x_{-r})\cdot\psi(s)\mathop{}\\!\mathrm{d}s\right\rvert\leq{C_{\sharp}}\|f\|_{\infty}e^{-h_{\operatorname{top}}r}\leq{C_{\sharp}}\|f\|_{\infty}t^{-d},$

where we used the fact that
$\int_{0}^{1}{\mathcal{L}}_{t-r}(|1|)\circ\mathsf{h}_{s}(y)\mathop{}\\!\mathrm{d}s\leq
C_{\tau}e^{h_{\operatorname{top}}(t-r)}$ for any $y\in{\widetilde{M}}$.

As in Section 7.2, let $\\{\rho_{j}\\}_{j\in\mathcal{J}}$ be a smooth
partition of unity on the horocycle segment starting at $x_{-r}$ of length
$\tau$ so that each $\rho_{j}$ is supported on a subarc $I_{j}$ of length 1
starting at $x_{j}$, and
$|\mathcal{J}|\leq{C_{\sharp}}\tau\leq{C_{\sharp}}\sigma
e^{h_{\operatorname{top}}r}$. With the same decomposition of $f$ as in (47),
we obtain

$\begin{split}\int_{0}^{\tau}\hat{{\mathcal{L}}}_{t-r}(f)\circ\mathsf{h}_{s}(x_{-r})\cdot\psi(s)\mathop{}\\!\mathrm{d}s&=\sum_{j\in\mathcal{J}}\int_{{\mathbb{T}}^{d}}\int_{0}^{1}e^{-2\pi\imath\xi_{\omega}(\mathsf{h}_{s}(x_{j}))}\cdot\hat{{\mathcal{L}}}^{(\omega)}_{t-r}u_{\omega}\circ\mathsf{h}_{s}(x_{j})\cdot\psi_{j}(s)\cdot\rho_{j}(s)\mathop{}\\!\mathrm{d}s\mathop{}\\!\mathrm{d}\omega\\\
&=\sum_{j\in\mathcal{J}}\int_{{\mathbb{T}}^{d}}e^{-2\pi\imath\xi_{\omega}(x_{j})}\cdot\mathcal{I}[p(x_{j}),\Psi_{j}(\cdot,\omega)]\Big{(}\hat{{\mathcal{L}}}^{(\omega)}_{t-r}u_{\omega}\Big{)}\mathop{}\\!\mathrm{d}\omega,\end{split}$

where $\psi_{j}$ is the restriction of $\psi$ to the arc $I_{j}$ and, as in
(50), we defined

$\Psi_{j}(s,\omega):=\exp\left(-2\pi\imath\int_{x_{j}}^{\mathsf{h}_{s}(x_{j})}p^{\ast}\omega\right)\cdot\psi_{j}(s)\cdot\rho_{j}(s).$

Clearly, $t_{0}:=t-r>t/2$; hence, following the same steps as in the previous
subsection, we use the spectral decomposition from Proposition 6.5 to reduce
the integral above to the following: define

$A_{j}(\omega):=\mathcal{I}[p(x_{j}),\Psi_{j}(\cdot,\omega)]\left(\Pi_{\omega}u_{\omega}\right),$

then,

$\left\lvert
e^{-h_{\operatorname{top}}t}\int_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}}\eta-e^{-h_{\operatorname{top}}r}\sum_{j\in\mathcal{J}}\int_{B(0,\delta)}e^{-2\pi\imath\xi_{\omega}(x_{j})+z(\omega)t_{0}}A_{j}(\omega)\mathop{}\\!\mathrm{d}\omega\right\rvert\leq{C_{\sharp}}C(f,x)t^{-d}.$

We focus on one fixed $j\in\mathcal{J}$. By Lemma 7.5, we have

$\left\lvert\int_{B(0,\delta)}e^{-2\pi\imath\xi_{\omega}(x_{j})+z(\omega)t_{0}}A_{j}(\omega)\mathop{}\\!\mathrm{d}\omega-
A_{j}(0)\int_{B(0,\delta)}e^{-2\pi\imath\xi_{\omega}(x_{j})-2\pi^{2}t_{0}\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega\right\rvert\\\
\leq{C_{\sharp}}C(f,x)t_{0}^{-\frac{d+1}{2}}\leq{C_{\sharp}}C(f,x)t^{-\frac{d+1}{2}},$

where $\Sigma$ is once again the positive definite matrix associated to the
quadratic form $-\frac{1}{4\pi^{2}}\sigma$ in (45).

We notice that we can restrict the integral above to $B(0,\delta
t_{0}^{-1/4})$. For all $\omega\in B(0,\delta t_{0}^{-1/4})$, we have

$\left\lvert\xi_{\omega}(x_{j})\right\rvert\leq{C_{\sharp}}\|\omega\|(r+\sigma)\leq{C_{\sharp}}\frac{\log
t}{t^{1/4}},$

where we used the fact that $x_{j}$ lies on the segment
$\mathsf{g}_{-r}\circ\gamma_{x,\sigma}$. Therefore,
$|e^{-2\pi\imath\xi_{\omega}(x_{j})}-1|\leq{C_{\sharp}}\|\omega\|\cdot\log t$
and hence

$\left\lvert\int_{B(0,\delta)}e^{-2\pi\imath\xi_{\omega}(x_{j})+z(\omega)t_{0}}A_{j}(\omega)\mathop{}\\!\mathrm{d}\omega-
A_{j}(0)\int_{B(0,\delta
t_{0}^{-1/4})}e^{-2\pi^{2}t_{0}\omega\cdot\Sigma\omega}\mathop{}\\!\mathrm{d}\omega\right\rvert\\\
\leq{C_{\sharp}}C(f,x)t_{0}^{-\frac{d+1}{2}}\leq{C_{\sharp}}C(f,x)\frac{\log
t}{t^{\frac{d+1}{2}}}.$

The integral on the right hand side above is computed as at the end of Section
7.2.

Combining everything together, and using the fact that
$\sum_{j}A_{j}(0)=\mu(f)\cdot\int_{0}^{\tau}\psi(s)\mathop{}\\!\mathrm{d}s$,
we conclude

$\left\lvert
e^{-h_{\operatorname{top}}t}\int_{\mathsf{g}_{-t}\circ\gamma_{x,\sigma}}\eta-\mu(f)\frac{e^{-h_{\operatorname{top}}r}\cdot\tau}{(2\pi
t_{0})^{\frac{d}{2}}\sqrt{\det\Sigma}}\right\rvert\leq{C_{\sharp}}C(f,x)\frac{\log
t}{t^{\frac{d+1}{2}}}.$

Finally, we apply Proposition 7.3 to the point $\mathsf{g}_{-r}(x)$ with $t=r$
so that we obtain

$|\sigma-e^{-h_{\operatorname{top}}r}\cdot\tau|\leq{C_{\sharp}}e^{-\delta_{*}r}\leq{C_{\sharp}}t^{-d}.$

This finishes the proof of B.

##  Appendix A Norm estimates for the winding cycle

This Appendix is devoted to the proof of some bounds regarding the winding
function defined in (12). Firstly, it is convenient to introduce some
notations. For each $\omega\in{\mathbb{T}}^{d}$, we define

(56)
$\begin{split}&A(X,\omega):=\max\\{\|\Phi^{-}\|_{{\mathscr{C}}^{1}}\|\langle
X,\omega\rangle\|_{{\mathscr{C}}^{1}},\|\langle
X,\omega\rangle\|_{{\mathscr{C}}^{1}}^{2},\|\langle
X,\omega\rangle\|_{{\mathscr{C}}^{2}}\\}\\\
&B(\bar{k}):=\max\\{(-\overline{k})^{-\frac{1}{2}},(-\overline{k})^{-1}\\},\end{split}$

where ${\overline{k}}<0$ is given in Lemma 3.2 and $\Phi^{-}$ in (7).

######  Proposition A.1.

Setting

$G_{t,\omega}(x)=\exp\left(2\pi\imath\int_{0}^{t}\langle\omega,X\rangle\circ\mathsf{g}_{-a}(x)\mathop{}\\!\mathrm{d}a\right),$

then for each $t\in{\mathbb{R}}^{+}$,
$I\subset{\mathcal{I}}_{\rho},\omega\in{\mathbb{T}}^{d},\psi\in{\mathscr{C}}^{r}_{c}(I)$
with $\|\psi\|_{{\mathscr{C}}^{r}(I)}\leq C_{*}$ and
$\varphi\in{\mathscr{C}}^{1+\alpha}_{c}(I)$, we have

$\|(\varphi
G_{t,\omega})\circ\mathsf{g}_{t}\circ\mathsf{h}_{(\cdot)}\circ\psi\|_{{\mathscr{C}}^{\varsigma}(I)}\leq
C_{\varsigma}(t)\|\varphi\|_{{\mathscr{C}}^{\varsigma}(I)},\quad\varsigma\in\\{0,\alpha,1+\alpha\\},$

where

(57) $C_{\varsigma}(t)=\begin{cases}1\qquad&\text{if}\qquad\varsigma=0,\\\
2\pi
C_{*}\big{(}\frac{1}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}+e^{-\sqrt{-\overline{k}}\alpha
t}\big{)}\qquad&\text{if}\qquad\varsigma=\alpha,\\\
C_{*}+8\pi^{2}\max\\{C_{*}^{3},C_{*}^{4}\\}B(\overline{k})A(X,\omega)\big{(}1+4e^{-\sqrt{-\overline{k}}(1+\alpha)t}\big{)}\qquad&\text{if}\qquad\varsigma=1+\alpha.\end{cases}$

###### Proof.

The case $\varsigma=0$ is trivial since
$\|G_{t,\omega}\|_{{\mathscr{C}}^{0}}=1$. Let us use the symbol ′ to denote
the derivative with respect to $\eta\in I$. We have

$\begin{split}((\varphi
G_{t,\omega})\circ&\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}\\\
=&(G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}\cdot(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)+G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\cdot(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}\\\
=&2\pi\imath\psi^{\prime}\left(\int_{0}^{t}(D\mathsf{g}_{t-a}(U)\langle\omega,X\rangle)\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi\mathop{}\\!\mathrm{d}a\right)(G_{t,\omega}\cdot\varphi)\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\\\
&+G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\cdot\psi^{\prime}(D\mathsf{g}_{t}(U)\varphi)\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi.\end{split}$

Recalling that $D\mathsf{g}_{t}U=J_{t}U$ and $J_{-t}=D\mathsf{g}_{-t}(U)$, we
note that

(58)
$(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}=\psi^{\prime}(D\mathsf{g}_{t}(U)\varphi)\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi=\psi^{\prime}J_{t}(U\varphi)\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi.$

Therefore, since $\|\psi\|_{{\mathscr{C}}^{r}}\leq C_{*}$ and
$\|G_{t,\omega}\|_{{\mathscr{C}}^{0}}=1$, we have

$\begin{split}|((\varphi
G_{t,\omega})\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}|&\leq
2\pi
C_{*}\|\varphi\|_{{\mathscr{C}}^{0}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}\int_{0}^{t}|J_{t-a}\circ\mathsf{h}_{\eta}\circ\psi|\mathop{}\\!\mathrm{d}a+C_{*}|J_{t}\circ\mathsf{h}_{\eta}|\|\varphi\|_{{\mathscr{C}}^{1}}.\end{split}$

By Lemma 3.2, $|J_{t}(x)|\leq e^{-\sqrt{-\overline{k}}t}$ for each $x\in M$,
and

(59) $\left|\int_{0}^{t}J_{t-a}\circ
h_{\eta}\mathop{}\\!\mathrm{d}a\right|\leq\int_{0}^{t}e^{-\sqrt{-\overline{k}}(t-a)}\mathop{}\\!\mathrm{d}a\leq\dfrac{1}{\sqrt{-\overline{k}}},$

therefore,

(60) $\begin{split}\|(\varphi
G_{t,\omega})\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1}}&\leq\frac{2\pi
C_{*}}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}\|\varphi\|_{{\mathscr{C}}^{0}}+C_{*}e^{-\sqrt{-\overline{k}}t}\|\varphi\|_{{\mathscr{C}}^{1}}\\\
&\leq
C_{*}\left(\frac{2\pi}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}+e^{-\sqrt{-\overline{k}}t}\right)\|\varphi\|_{{\mathscr{C}}^{1}}.\end{split}$

For $\varsigma=\alpha$ it is enough to observe that, for each
$\eta,\tilde{\eta}\in(0,\rho)$,

(61)
$\begin{split}|(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{(\cdot)}\circ\psi)(\tilde{\eta})-(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{(\cdot)}\circ\psi_{\tilde{\eta}})(\eta)|&\leq
d_{I}(\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi,\mathsf{g}_{t}\circ\mathsf{h}_{\tilde{\eta}}\circ\psi)^{\alpha}\|\varphi\|_{{\mathscr{C}}^{\alpha}}\\\
&\leq C_{*}e^{-\sqrt{-\overline{k}}\alpha
t}|\eta-\tilde{\eta}|^{\alpha}\|\varphi\|_{{\mathscr{C}}^{\alpha}},\end{split}$

where $d_{I}(x,y)$ is the distance of two points $x,y$ along $I$. Hence, since
$G_{t,\omega}\in{\mathscr{C}}^{r}$, using the inequality

(62) $\|\varphi
G\|_{{\mathscr{C}}^{\alpha}}\leq\|\varphi\|_{{\mathscr{C}}^{0}}\|G\|_{{\mathscr{C}}^{1}}+\|\varphi\|_{{\mathscr{C}}^{\alpha}}\|G\|_{{\mathscr{C}}^{0}},$

we have

$\begin{split}\|(\varphi
G_{t,\omega})\circ&\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{\alpha}(I)}\leq\\\
&\leq\|G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1}}\|\varphi\|_{{\mathscr{C}}^{0}}+\|G_{t,\omega}\|_{{\mathscr{C}}^{0}}\|\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{\alpha}}\\\
&\leq
C_{*}\left(\frac{2\pi}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}+e^{-\sqrt{-\overline{k}}t}\right)\|\varphi\|_{{\mathscr{C}}^{0}}+C_{*}e^{-\sqrt{-\overline{k}}\alpha
t}\|\varphi\|_{{\mathscr{C}}^{\alpha}}\\\ &\leq
C_{*}\left(\frac{2\pi}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}+e^{-\sqrt{-\overline{k}}\alpha
t}\right)\|\varphi\|_{{\mathscr{C}}^{\alpha}}.\end{split}$

It remains the case $\varsigma=1+\alpha$. Using (62), we can check that if
$G\in{\mathscr{C}}^{2}$ and $\varphi\in{\mathscr{C}}^{1+\alpha}$, then

(63) $\|\varphi
G\|_{{\mathscr{C}}^{1+\alpha}}\leq\|\varphi\|_{{\mathscr{C}}^{1+\alpha}}(1+4\|G\|_{{\mathscr{C}}^{2}}).$

Let us thus compute the second derivative of
$(G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)$: using
again that $D\mathsf{g}_{t}U=J_{t}U$,

(64)
$\begin{split}|(G_{t,\omega}\circ&\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime\prime}|\\\
=&\left|\left(2\pi\imath(G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)\psi^{\prime}\int_{0}^{t}J_{t-a}\circ\mathsf{h}_{\eta}\cdot
U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi)\mathop{}\\!\mathrm{d}a\right)^{\prime}\right|\\\
\leq&2\pi
C_{*}\bigg{(}\left|\int_{0}^{t}(J_{t-a}\circ\mathsf{h}_{\eta})^{\prime}\cdot
U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi)\mathop{}\\!\mathrm{d}a\right|\\\
&+\left|\int_{0}^{t}J_{t-a}\circ\mathsf{h}_{\eta}\cdot(U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi))^{\prime}\mathop{}\\!\mathrm{d}a\right|\bigg{)}\\\
&+2\pi C_{*}\left|\int_{0}^{t}J_{t-a}\circ\mathsf{h}_{\eta}\cdot
U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi)\mathop{}\\!\mathrm{d}a\right|(1+\|G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1}}).\end{split}$

We need to estimate the above terms. By (5), (6), (7) and (59),

$\begin{split}\left|(J_{t-a}\circ{\mathsf{h}_{\eta}})^{\prime}\right|&\leq\left|\int_{0}^{t}(\Phi^{-}\circ\mathsf{g}_{a}\circ\mathsf{h}_{\eta})^{\prime}\right|\left|J_{t-a}\circ\mathsf{h}_{\eta}\right|\leq\dfrac{\|\Phi^{-}\|_{{\mathscr{C}}^{1}}}{\sqrt{-\overline{k}}}e^{-\sqrt{-\overline{k}}(t-a)}.\end{split}$

On the other hand

$\left|(U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi))^{\prime}\right|\leq
C_{*}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{2}}\|J_{t-a}\|_{{\mathscr{C}}^{0}}\leq
C_{*}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{2}}e^{-\sqrt{-\overline{k}}(t-a)}$

and

$\int_{0}^{t}|U(\langle\omega,X\rangle\circ\mathsf{g}_{t-a}\circ\mathsf{h}_{\eta}\circ\psi)|\mathop{}\\!\mathrm{d}a\leq
C_{*}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}\int_{0}^{t}e^{-\sqrt{-\overline{k}}(t-a)}\mathop{}\\!\mathrm{d}a.$

Using the above three inequalities into (64), and recalling (60), we obtain

$\begin{split}\|G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{2}}&\leq
2\pi
C_{*}^{2}\dfrac{\|\Phi^{-}\|_{{\mathscr{C}}^{1}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}}{2{\sqrt{-\overline{k}}}}+\frac{2\pi
C_{*}^{2}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{2}}}{2\sqrt{-\overline{k}}}\\\
&+\frac{2\pi
C_{*}^{2}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}}{2\sqrt{-\overline{k}}}\left(1+C_{*}\left(\frac{2\pi}{\sqrt{-\overline{k}}}\|\langle\omega,X\rangle\|_{{\mathscr{C}}^{1}}+e^{-\sqrt{-\overline{k}}t}\right)\right).\end{split}$

We have thus obtained

$\|G_{t,\omega}\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{2}}\leq
2\pi^{2}\max\left\\{\frac{C_{*}^{2}}{\sqrt{-\overline{k}}},\frac{C_{*}^{3}}{-\overline{k}}\right\\}A(X,\omega)\left(4+e^{-\sqrt{-\overline{k}}t}\right).$

Hence, by (63), we have

$\|(\varphi
G_{t,\omega})\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1+\alpha}}\\\
\leq(1+8\pi^{2}\max\\{C_{*}^{2},C_{*}^{3}\\}B(\overline{k})A(X,\omega)(4+e^{-\sqrt{-\overline{k}}t}))\|\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1+\alpha}}.$

Finally, by (58), and arguing as in (61),

$\|(\varphi\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi)^{\prime}\|_{{\mathscr{C}}^{\alpha}}=\|\psi^{\prime}J_{t}(U\varphi)\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{\alpha}}\leq
C_{*}e^{-\sqrt{\overline{k}}t}e^{-\sqrt{-\overline{k}}\alpha
t}\|\varphi^{\prime}\|_{{\mathscr{C}}^{\alpha}}.$

Inserting the above in the previous equation we conclude:

$\|(\varphi
G_{t,\omega})\circ\mathsf{g}_{t}\circ\mathsf{h}_{\eta}\circ\psi\|_{{\mathscr{C}}^{1+\alpha}}\\\
\leq\left(C_{*}+8\pi^{2}\max\\{C_{*}^{3},C_{*}^{4}\\}B(\overline{k})A(X,\omega)\big{(}1+4e^{-\sqrt{-\overline{k}}(1+\alpha)t}\big{)}\right)\|\varphi\|_{{\mathscr{C}}^{1+\alpha}}.$

The proof is therefore complete. ∎

### Acknowledgements

This work was started when RC was working at University of Pisa and it is
partially supported by the research project PRIN 2022NTKXCX “Stochastic
properties of dynamical systems” funded by the Ministry of University and
Scientific Research of Italy.
This research is part of RC’s activity within the INdAM (Istituto Nazionale di
Alta Matematica) group GNFM, and RC and DR’s activity within the UMI Group
“DinAmicI”.
It is a pleasure to thank Claudio Bonanno, Paolo Giulietti, Carlangelo
Liverani, Omri Sarig, and Milo Viviani for many useful discussions.

## References

  * [1] J. Aaronson. An Introduction to Infinite Ergodic Theory. Math. Surveys and Monographs 50, Amer. Math. Soc. (1997).
  * [2] A. Adam, V. Baladi. Horocycle averages on closed manifolds and transfer operators. Tunisian J. Math 4 387–441 (2022).
  * [3] A. Avila, D. Dolgopyat, E. Duryev, and O. Sarig. The visits to zero of a random walk driven by an irrational rotation. Isr. J. Math. 207, 653–717 (2015).
  * [4] M. Blank, G. Keller, and C. Liverani, Ruelle-Perron-Frobenius spectrum for Anosov maps, Nonlinearity 15 (2002), 1905-1973. MR 1938476. Zbl 1021. 37015.
  * [5] R. Bowen. Periodic orbits for hyperbolic flows, Amer. J. Math. 94 (1972), 1–30.
  * [6] H. Bruin, C. Fougeron, D. Ravotti, and D. Terhesiu. On asymptotic expansions of ergodic integrals for ${\mathbb{Z}}^{d}$-extensions of translation flows. Preprint arXiv:2402.02266 (2024).
  * [7] O. Butterley. A note on operator semigroups associated to chaotic flows–corrigendum. Ergodic Theory Dynam. Systems 36.5, pp. 1409–1410, (2016).
  * [8] O. Butterley. Expanding semiflows on branched surfaces and one-parameter semigroups of operators. Nonlinearity 25(12) (2012), 3487–3503.
  * [9] O. Butterley and C. Liverani. Smooth Anosov flows: correlation spectra and stability. J. Mod. Dyn. 1(2) (2007), 301–322.
  * [10] O. Butterley and L.D. Simonelli. Parabolic flows renormalized by partially hyperbolic maps. Bollettino dell’Unione Matematica Italiana 13, 341–360 (2020).
  * [11] R. Castorrini, C. Liverani. Quantitative statistical properties of two-dimensional partially hyperbolic systems. Adv. Math., 409, Part A, 1–122 (2022).
  * [12] B. Davies. Linear operators and their spectra. Vol. 106. Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, pp. xii+451 (2007).
  * [13] M. F. Demers, N. Kiamari, C. Liverani, Transfer operators in hyperbolic dynamics. An introduction. $33^{o}$ Colóq. Bras. Mat. Instituto Nacional de Matemática Pura e Aplicada (IMPA), Rio de Janeiro, 2021.
  * [14] M.P. do Carmo. Riemannian Geometry. Mathematics: Theory and Applications. Birkhäuser Boston, MA (1992).
  * [15] F. Faure, S. Gouezel, and E. Lanneau. Ruelle spectrum of linear pseudo-Anosov maps. J. Éc. Polytech. Math. 6, 811-877 (2019).
  * [16] L. Flaminio, G. Forni. Invariant distributions and time averages for horocycle flows. Duke Math. J. 119(3) (2003), 465–526.
  * [17] L. Flaminio, D. Ravotti. Abelian covers of hyperbolic surfaces: equidistribution of spectra and infinite mixing asymptotics for horocycle flows. Preprint arXiv:2207.08263 (2022).
  * [18] P. Giulietti and C. Liverani, Parabolic dynamics and anisotropic Banach spaces, JEMS 21 2793–2858 (2019).
  * [19] P. Giulietti, C. Liverani, and M. Pollicott, Anosov flows and dynamical zeta functions, Ann. of Math. (2) 178:2 (2013), 687–773.
  * [20] S. Gouëzel , Limit theorems in dynamical systems using the spectral method (Proceedings of Symposia in Pure Mathematics 89:161–193, 2015)
  * [21] S. Gouëzel and C. Liverani, Banach spaces adapted to Anosov systems, Er- godic Theory Dynam. Systems 26 (2006), 189–217.
  * [22] L.W. Green, Remarks on uniformly expanding horocycle reparametrization, J. Differential Geometry, 13 (1978) 263–271.
  * [23] B.M. Gurevic. The entropy of a horocycle flow, Soviet Math. Doklady 2 (1961), 124–130.
  * [24] Y. Guivarc’h, Y. Le Jan. Asymptotic winding of the geodesic flow on modular surfaces and continuous fractions. Annales scientifiques de l’ École Normale Supérieure, Serie 4, Volume 26 (1993) no. 1, pp. 23–50.
  * [25] F. Ledrappier, O. Sarig. Unique ergodicity for non-uniquely ergodic horocycle flows. Discrete and Continuous Dynamical Systems, 16(2) (2006), 411–433.
  * [26] F. Ledrappier, O. Sarig. Fluctuations of ergodic sums for horocycle flows on ${\mathbb{Z}}^{d}$-covers of finite volume surfaces. Discrete Contin. Dyn. Syst., 22(1-2) (2008), 247–325.
  * [27] C. Liverani, On contact Anosov flows, Ann. of Math. (2) 159:3 (2004), 1275–1312.
  * [28] B. Hasselblatt, A. Katok. Introduction to the Modern Theory of Dynamical Systems. Encyclopedia of Mathematics and its Applications, Issue 54, Cambridge University Press (1995).
  * [29] H. Hennion. Sur un théorème spectral et son application aux noyaux lipchitziens. Proc. Amer. Math. Soc. 118:2 (1993), 627–634.
  * [30] S. Hurder, A. Katok. Differentiability, rigidity and Godbillon-Vey classes for Anosov flows. Publications mathématiques de l’IHES., tome 72 (1990), p. 5–61.
  * [31] T. Kato. Perturbation theory for linears operators. Grundlehren der Mathematicschen Wissenschaften 132 (1996).
  * [32] A. Katsuda, T. Sunada. Closed orbits in homology classes. Publications Mathematiques de l’IHES, Volume 71 (1990), pp. 5–32.
  * [33] G.A. Margulis. Certain measures associated with U-flows on compact manifolds. Functional Anal. Appl. 4 (1) (1970), 55–67.
  * [34] B. Marcus. Unique ergodicity of the horocycle flow: variable negative curvature case. Israel J. Math., 21(2-3) (1975), 133–144.
  * [35] B. Marcus. Ergodic Properties of Horocycle Flows for Surfaces of Negative Curvature. Annals of Mathematics, Second Series, Vol. 105, No. 1 (1977), pp. 81–105.
  * [36] J. Plante. Anosov flows, Amer. J. Math. 94 (1972), 729–754.
  * [37] L. Schwartz, Théorie Des Distributions, Publications de l’Institut de Mathématique de l’Université de Strasbourg, Hermann: Paris, 1966.
  * [38] M. Tsujii, Quasi-compactness of transfer operators for contact Anosov flows, Nonlinearity 23:7 (2010), 1495–1545.


MIDLMedical Imaging with Deep Learning 2022 \- Dave Van Veen, Rogier van der
Sluijs, Batu Ozturkler, Arjun Desai, Christian Bluethgen, Robert D. Boutin,
Marc H. Willis, Gordon Wetzstein, David Lindell, Shreyas Vasanawala, John
Pauly, Akshay S. Chaudhari
{vanveen, sluijs, ozt, arjundd, bluch, boutin, mwillis7,
gordonwz, lindell, vasanawala, pauly<EMAIL_ADDRESS>
Stanford University

# Scale-Agnostic Super-Resolution in MRI using Feature-Based Coordinate
Networks

###### Abstract

We propose using a coordinate network decoder for the task of super-resolution
in MRI. The continuous signal representation of coordinate networks enables
this approach to be scale-agnostic, i.e. one can train over a continuous range
of scales and subsequently query at arbitrary resolutions. Due to the
difficulty of performing super-resolution on inherently noisy data, we analyze
network behavior under multiple denoising strategies. Lastly we compare this
method to a standard convolutional decoder using both quantitative metrics and
a radiologist study implemented in Voxel111https://voxel.im, our newly
developed tool for web-based evaluation of medical images.

###### keywords:

Coordinate networks, super-resolution, MRI.

##  1 Introduction

High-resolution scans in magnetic resonance imaging (MRI) to depict high-
frequency details such as local textures and edges are crucial for many
diagnostic imaging tasks. However, tradeoffs with scan time and signal-to-
noise ratios motivate improved MRI resolution for higher downstream clinical
utility. While deep learning can leverage data-driven priors for encoding
high-frequency information in super-resolution tasks, state-of-the-art methods
are limited to upsampling at fixed, discrete scales due to their convolutional
structure [Lim et al. (2017), Chen et al. (2021a)]. Discrete scales are
undesirable for clinical interpretation [Chaudhari et al. (2021)]. Further,
training such fixed networks places strict limits on acquiring homogeneous
training data. In clinical settings, resolution for MRI acquisitions often
varies based on factors such as patient size [Wargo et al. (2013)]. Hence to
construct a training dataset with uniform resolution, individual images must
be resized, which can lead to poorer model performance [Karras et al. (2019),
Chai et al. (2022)].

We propose a scale-agnostic framework for MRI super-resolution using a
coordinate network as a decoder. Given a coordinate, this decoder queries the
neighboring latent features to predict the pixel value at that location. The
continuous nature of this decoder allows querying at arbitrary resolutions; it
also decouples the training and querying scales, e.g. one can train on a
continuous range of 1-2$\times$ and query at 3$\times$. Additionally, since
super-resolution tasks are prone to augment noise [Singh et al. (2014)], we
analyze the behavior of coordinate networks with various denoising strategies.
This demonstrates the importance of early stopping, a technique commonly used
in convolutional networks to avoid overfitting [Heckel and Soltanolkotabi
(2019)]. Finally we compare the proposed framework’s coordinate decoder
against a standard convolutional decoder, using image quality metrics and a
clinical reader study.

### Contributions

  * •

We propose a framework for MRI super-resolution that is scale-agnostic. This
enables both decoupling between training and querying scales and also querying
at arbitrary resolutions.

  * •

We demonstrate various regularization strategies for coordinate networks in
the presence of noisy data, including early stopping and the use of denoisers
during network training.

  * •

In addition to quantitative metrics, we evaluate results with a radiologist
reader study built on Voxel, our newly developed tool for web-based evaluation
of medical images.

##  2 Related Work

Coordinate networks, sometimes called implicit neural representations, are
recent powerful tools for representing signals such as images with fully-
connected multi-layer perceptrons (MLPs) by mapping image coordinates to their
corresponding pixel values. In contrast to standard pixel-based
representations, this representation is continuous with respect to network
weights, allowing it to model fine detail which is limited by network capacity
instead of grid resolution. Coordinate networks are commonly employed for
unsupervised signal representation by training weights such that the network
output has high fidelity with a single signal of interest [Sitzmann et al.
(2020), Tancik et al. (2020), Lindell et al. (2022)]. These techniques have
also been demonstrated for representing medical images [Wu et al. (2021), Reed
et al. (2021)]. However, unsupervised approaches cannot incorporate novel
high-frequency information for the task of super-resolution. In contrast,
supervised methods learn to represent many signals over a shared function
space, often in concordance with meta-learning [Tancik et al. (2021)], learned
initializations [Shen et al. (2022)], or convolutional structure to encode
features [Chen et al. (2021b), Mehta et al. (2021)]. These supervised methods
are developed for smooth images and do not account for noise present in real-
world medical imaging applications. While it’s been shown that convolutional
networks are prone to overfitting and hence benefit from early stopping [Van
Veen et al. (2018), Heckel and Soltanolkotabi (2019)], the behavior of
coordinate networks in noisy settings is a developing area of study.
Concurrently with our work, Kim et al. (2022a) emphasize the importance of
early stopping for coordinate networks in an unsupervised setting.
Tangentially, we note another concurrent work which demonstrates a scale-
agnostic method for super-resolution of audio signals Kim et al. (2022b).

##  3 Methods

fig:diagram

Figure 1:  Left: System overview, described in Section 3.1. Right: Image
comparison given input cropped from the red box at left. Coord and conv are
each trained at 2$\times$. Coord can also be queried at 4$\times$ without re-
training because it is scale-agnostic (F). Coord benefits from the denoising
regularization term i.e. when $\lambda\neq 0$ (D vs. E), and conv does not
(see Table LABEL:tab:example).

###  3.1 System Overview

We now provide a description of our method following the diagram in Figure
LABEL:fig:diagram. High-resolution image
$\textbf{x}_{hr}\in\mathbb{R}^{n\times n}$ is obtained from fully-sampled
k-space measurements y characterized by the forward model
$\textbf{y}_{i}=\textbf{F}\textbf{S}_{i}\textbf{x}_{hr}$, where F denotes the
Fourier transform and $\textbf{S}_{i}$ the $i^{th}$ coil sensitivity map.
Given $\textbf{x}_{hr}$, a low-resolution network input
$\textbf{x}_{lr}\in\mathbb{R}^{\frac{n}{s}\times\frac{n}{s}}$ is created via
bicubic downsampling by scale $s$. Subsequently an encoder
$f_{\psi}:\textbf{x}_{lr}\rightarrow\textbf{C}\in\mathbb{R}^{d\times l\times
w}$ parameterized by weights $\psi$ maps this to an equidistant feature grid.
Hence there are $l\times w$ grid points, each containing a feature vector
$\textbf{c}\in\mathbb{R}^{d}$.

These feature vectors are then input to the coordinate decoder
$g_{\theta}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ defined by:

$g_{\theta}\coloneqq\textbf{W}_{m}(\alpha_{m-1}\circ\alpha_{m-2}\circ\cdots\circ\alpha_{0})(\textbf{c})+\textbf{b}_{m},\;\;\;\;\textbf{c}_{j}\rightarrow\alpha_{j}(\textbf{c}_{j})\coloneqq\phi(\textbf{W}_{j}\textbf{c}_{j}+\textbf{b}_{j}),$

where $\alpha_{j}$ corresponds to the $j^{th}$ network layer composed of
weight matrix $\textbf{W}_{j}$ and biases $\textbf{b}_{j}$ operating on the
input $\textbf{c}_{j}$, followed by application of a nonlinear function
$\phi$.

Given a particular feature vector $\textbf{c}^{*}\in\mathbb{R}^{d}$ at a
coordinate of interest, the decoder
$g_{\theta}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ queries the four neighboring
latent codes
$[\textbf{c}^{(1)},\textbf{c}^{(2)},\textbf{c}^{(3)},\textbf{c}^{(4)}]$ such
that the decoder’s output predicts the grayscale pixel value at each of those
four locations. The pixel value at $\textbf{c}^{*}$ is estimated as a linear
combination of surrounding pixel values based on relative distance
$\textbf{w}^{(k)}$, i.e.
$g_{\theta}(\textbf{c}^{*})=\sum_{k=1}^{4}\textbf{w}^{(k)}\cdot
g_{\theta}(\textbf{c}^{(k)})$, consequently preventing discontinuities in the
output image. Querying over many coordinates produces the predicted image
$\hat{\textbf{x}}$, so in summary we have $\hat{\textbf{x}}=(g_{\theta}\circ
f_{\psi})(\textbf{x}_{lr})$.

We train this system end-to-end with $p$ pairs
$\\{\textbf{x}_{lr},\textbf{x}_{hr}\\}_{i=1}^{p}$ to find encoder and decoder
weights $\psi$, $\theta$ which minimize the loss $L_{c}+\lambda
L_{d}=\|\hat{\textbf{x}}-\textbf{x}_{hr}\|_{1}+\lambda\|\hat{\textbf{x}}-D_{\sigma}(\textbf{x}_{hr})\|_{2}^{2}$.
Here $L_{c}$ is a consistency loss and $L_{d}$ a denoising loss, such that
$D_{\sigma}$ is a denoiser with strength $\sigma$ and $\lambda$ is the
relative weight of the regularization term. We use the Adam optimizer over a
batch of randomly drawn image tiles for each training step and perform early
stopping after $T$ iterations, as discussed below.

###  3.2 Implementation

This framework allows for many choices of encoders $f_{\psi}$ or denoisers
$D_{\sigma}$. For simplicity we choose an EDSR convolutional encoder Lim et
al. (2017) and BM3D denoiser Dabov et al. (2007), respectively. Both the
denoiser strength $\sigma$ and number of training steps $T$ were tuned
according to radiologist preference, resulting in $\sigma=0.03$,
$T_{\textit{coord}}=10^{3}$, and $T_{\textit{conv}}=10^{5}$. Alternatively one
could determine $T$ by maintaining a hold-out validation set, e.g. similar to
Yaman et al. (2020). However, this relies on quantitative metrics, which can
be ambiguous, as discussed in Section 4.

We compare our decoder’s continuous representation (“coord”), to the same
framework with a convolutional decoder, i.e. the original EDSR (“conv”), which
is not scale-agnostic. Coord is similar to conv but modifies the decoder to be
a five-layer MLP containing 256 hidden units and ReLU activations. Each model
contains roughly 1.6 million parameters and runs on a single Quadro RTX 8000
GPU. We also compare against bicubic interpolation, which can be queried at
arbitrary scales but has no prior to incorporate higher frequency information.
To create the dataset, we extract 2D sagittal slices from SKM-TEA Desai et al.
(2021), randomly sampling five of the central 80 slices from each 3D echo-1
scan and partitioning into a 80%/10%/10% train/test/validation split by
volume.

fig:voxel

Figure 2: User interface for Voxel, our web-based tool used in the reader
study.

###  3.3 Reader study details

We perform a reader study with radiologists comparing coord and conv, both
trained on 2$\times$ super-resolution. In clinical applications one would want
to scale beyond ground-truth; hence at inference we bypass downsampling and
perform 2$\times$ super-resolution on the ground-truth itself. Readers used a
five-point Likert scale to evaluate randomized side-by-side image pairs in
terms of sharpness and noise.

The reader study was built using Voxel, a medical imaging tensor viewer for
radiologist evaluation (see Fig. LABEL:fig:voxel). This allows the reader to
review images from her/his browser, forgoing the cumbersome file sharing and
annotation commonly required with reader studies. Voxel provides similar
features to standard DICOM viewers such as zooming, panning through image
volumes, and adjusting window levels. Furthermore, responses are time-stamped,
allowing for review-time analysis; while we do not investigate it here, this
could enable answering questions such as, “Would radiologists be faster
diagnosing images from method A or method B?” Broadly we hope this tool will
enhance the impact of deep learning in medical imaging by streamlining the
feedback process between researchers and radiologists. To try a demo of Voxel,
please visit https://voxel.im or contact the authors.

##  4 Results and Discussion

###  4.1 Quantitative Results

Per Table LABEL:tab:example, coord performs comparably when trained on a range
of scales (1-2$\times$, row 2) vs. a fixed scale (2$\times$, row 5). Because
coord is scale-agnostic, it can be queried at a resolution which is both
arbitrary and independent of its training scales. Conversely, conv—without
additional interpolation—queries only at fixed integer upsampling according to
its training scales. Conv was slightly superior to coord in PSNR but slightly
inferior in VIF. While PSNR is biased toward smooth images, VIF is perhaps
better suited to evaluate high-frequency detail augmented during super-
resolution; further, it is more indicative of clinical diagnostic quality
Mason et al. (2019).

tab:example

Table 1: Left: Quantitative scores (VIF/PSNR). Coord obtains similar
performance to conv (slightly better VIF, slightly worse PSNR) and has the
benefit of agnosticism with respect to training and querying scales. Below the
midline we provide ablations at various training scales and values for
$\lambda$, i.e. the denoising regularization term (see discussion, Section 4).
Right: Reader study scoring criteria (top) and results (bottom) demonstrating
a slight preference for coord.

fig:early-stop-coord

Figure 3: Early stopping is an important regularization technique for
coordinate networks. Quantitative metrics for coord are evaluated at each
training step. Similar to the known phenomena in convolutional networks,
coordinate networks overfit to high frequency components given too many
training steps (see image insets), corresponding to a significant decrease in
PSNR. When analyzing conv, early stopping shows a similar benefit, albeit less
severe.

###  4.2 Reader Study

We note the challenge of evaluating clinical potential using image metrics
alone, as they frequently do not correspond to the images preferred by end-
user radiologists Kastryulin et al. (2022). Consider instances where higher
metrics scores seemingly do not pertain to better quality: compared to coord,
$\lambda=10$ (Fig. LABEL:fig:diagram, E), bicubic interpolation (B) achieves
higher PSNR while coord, $\lambda=0$ (D) achieves higher VIF; however, these
are perceptually undesirable in terms of sharpness and noise, respectively.
Furthermore, quantitative metrics require a ground-truth reference; yet in a
clinical setting, the goal is to scale larger than ground-truth resolution.
Hence to gain insight beyond these limitations, we present a reader study in
Table LABEL:tab:example. This demonstrates that coord is equivalent or
slightly preferable to conv in terms of perceived sharpness and noise. It’s
interesting to consider that despite coord receiving slighly lower scores on
PSNR, radiologists actually prefer these images with respect to noise.
Fundamentally we advocate for the use of reader studies when comparing
reconstruction methods to facilitate more clinically relevant research. One
tool to streamline this process is Voxel, as discussed in Section 3.3.

###  4.3 Denoising

Given the noise inherent to MR images, super-resolution methods such as coord
and conv are prone to augment noise in the output Singh et al. (2014). To
mitigate this effect, we leverage two techniques: denoising regularization and
early stopping during network training. Unlike coord which benefited from
denoising regularization (Fig. LABEL:fig:diagram, D vs. E), conv did not
(Table LABEL:tab:example, row 6). Meanwhile both benefit from early stopping,
although the effect on coord (Fig. LABEL:fig:early-stop-coord) is larger than
that on conv (Appendix). We suspect this difference is due to the fact that
convolutional structure encodes global image features, inherently acting as a
denoiser. Coordinate networks do not exhibit this structural property;
furthermore, their receptive field–––size of the input region that produces
the output—is typically much smaller compared to convolutional networks.
Overall, the comparison of coordinate and convolutional networks as image
priors has not been well studied, although concurrent work has examined their
behavior in unsupervised regimes Kim et al. (2022a). Similar to our findings,
this demonstrates that coordinate networks have an inductive bias toward lower
frequencies and consequently benefit from early stopping.

##  5 Conclusion

We propose a scale-agnostic framework for MRI super-resolution to circumvent
the requirement of homogeneous training data in convolutional networks. We
also demonstrate various regularization strategies for coordinate networks in
the presence of noise. Lastly we evaluate results using both quantitative
metrics and a reader study built on Voxel.

In the future, we plan to extensively evaluate across different encoding,
decoding, and denoising methods and also assess downstream impact via pixel-
level quantitative MRI metrics and reader studies on various pathological
regions. Crucial to this work is the behavior of coordinate networks in the
presence of noise, which warrants further study both for supervised and
unsupervised regimes.

## References

  * Chai et al. (2022) Lucy Chai, Michael Gharbi, Eli Shechtman, Phillip Isola, and Richard Zhang.  Any-resolution training for high-resolution image synthesis.  _arXiv preprint arXiv:2204.07156_ , 2022.
  * Chaudhari et al. (2021) Akshay S Chaudhari, Christopher M Sandino, Elizabeth K Cole, David B Larson, Garry E Gold, Shreyas S Vasanawala, Matthew P Lungren, Brian A Hargreaves, and Curtis P Langlotz.  Prospective deployment of deep learning in mri: a framework for important considerations, challenges, and recommendations for best practices.  _Journal of Magnetic Resonance Imaging_ , 54(2):357–371, 2021.
  * Chen et al. (2021a) Hanting Chen, Yunhe Wang, Tianyu Guo, Chang Xu, Yiping Deng, Zhenhua Liu, Siwei Ma, Chunjing Xu, Chao Xu, and Wen Gao.  Pre-trained image processing transformer.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 12299–12310, 2021a.
  * Chen et al. (2021b) Yinbo Chen, Sifei Liu, and Xiaolong Wang.  Learning continuous image representation with local implicit image function.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8628–8638, 2021b.
  * Dabov et al. (2007) Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian.  Image denoising by sparse 3-d transform-domain collaborative filtering.  _IEEE Transactions on image processing_ , 16(8):2080–2095, 2007.
  * Desai et al. (2021) Arjun D Desai, Andrew M Schmidt, Elka B Rubin, Christopher Michael Sandino, Marianne Susan Black, Valentina Mazzoli, Kathryn J Stevens, Robert Boutin, Christopher Re, Garry E Gold, Brian Hargreaves, and Akshay Chaudhari.  SKM-TEA: A dataset for accelerated MRI reconstruction with dense image labels for quantitative clinical evaluation.  In _Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2)_ , 2021.
  * Heckel and Soltanolkotabi (2019) Reinhard Heckel and Mahdi Soltanolkotabi.  Denoising and regularization via exploiting the structural bias of convolutional generators.  _arXiv preprint arXiv:1910.14634_ , 2019.
  * Karras et al. (2019) Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 4401–4410, 2019.
  * Kastryulin et al. (2022) Segrey Kastryulin, Jamil Zakirov, Nicola Pezzotti, and Dmitry V Dylov.  Image quality assessment for magnetic resonance imaging.  _Medical Image Analysis_ , 2022.
  * Kim et al. (2022a) Chaewon Kim, Jaeho Lee, and Jinwoo Shin.  Zero-shot blind image denoising via implicit neural representations.  _arXiv preprint arXiv:2204.02405_ , 2022a.
  * Kim et al. (2022b) Jaechang Kim, Yunjoo Lee, Seunghoon Hong, and Jungseul Ok.  Learning continuous representation of audio for arbitrary scale super resolution.  In _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 3703–3707. IEEE, 2022b.
  * Lim et al. (2017) Bee Lim, Sanghyun Son, Heewon Kim, Seungjun Nah, and Kyoung Mu Lee.  Enhanced deep residual networks for single image super-resolution.  In _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , pages 136–144, 2017.
  * Lindell et al. (2022) David B. Lindell, Dave Van Veen, Jeong Joon Park, and Gordon Wetzstein.  Bacon: Band-limited coordinate networks for multiscale scene representation.  In _CVPR_ , 2022.
  * Mason et al. (2019) Allister Mason, James Rioux, Sharon E Clarke, Andreu Costa, Matthias Schmidt, Valerie Keough, Thien Huynh, and Steven Beyea.  Comparison of objective image quality metrics to expert radiologists’ scoring of diagnostic quality of mr images.  _IEEE transactions on medical imaging_ , 39(4):1064–1072, 2019.
  * Mehta et al. (2021) Ishit Mehta, Michaël Gharbi, Connelly Barnes, Eli Shechtman, Ravi Ramamoorthi, and Manmohan Chandraker.  Modulated periodic activations for generalizable local functional representations.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 14214–14223, 2021.
  * Reed et al. (2021) Albert W Reed, Hyojin Kim, Rushil Anirudh, K Aditya Mohan, Kyle Champley, Jingu Kang, and Suren Jayasuriya.  Dynamic ct reconstruction from limited views with implicit neural representations and parametric motion fields.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 2258–2268, 2021.
  * Shen et al. (2022) Liyue Shen, John Pauly, and Lei Xing.  Nerp: implicit neural representation learning with prior embedding for sparsely sampled image reconstruction.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2022\.
  * Singh et al. (2014) Abhishek Singh, Fatih Porikli, and Narendra Ahuja.  Super-resolving noisy images.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2846–2853, 2014.
  * Sitzmann et al. (2020) Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  _Advances in Neural Information Processing Systems_ , 33:7462–7473, 2020.
  * Tancik et al. (2020) Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  _Advances in Neural Information Processing Systems_ , 33:7537–7547, 2020.
  * Tancik et al. (2021) Matthew Tancik, Ben Mildenhall, Terrance Wang, Divi Schmidt, Pratul P Srinivasan, Jonathan T Barron, and Ren Ng.  Learned initializations for optimizing coordinate-based neural representations.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2846–2855, 2021.
  * Van Veen et al. (2018) Dave Van Veen, Ajil Jalal, Mahdi Soltanolkotabi, Eric Price, Sriram Vishwanath, and Alexandros G Dimakis.  Compressed sensing with deep image prior and learned regularization.  _arXiv preprint arXiv:1806.06438_ , 2018.
  * Wargo et al. (2013) Christopher J Wargo, Jay Moore, and John C Gore.  A comparison and evaluation of reduced-fov methods for multi-slice 7 t human imaging.  _Magnetic resonance imaging_ , 31(8):1349–1359, 2013.
  * Wu et al. (2021) Qing Wu, Yuwei Li, Lan Xu, Ruiming Feng, Hongjiang Wei, Qing Yang, Boliang Yu, Xiaozhao Liu, Jingyi Yu, and Yuyao Zhang.  Irem: High-resolution magnetic resonance image reconstruction via implicit neural representation.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 65–74. Springer, 2021.
  * Yaman et al. (2020) Burhaneddin Yaman, Seyed Amir Hossein Hosseini, Steen Moeller, Jutta Ellermann, Kâmil Uğurbil, and Mehmet Akçakaya.  Self-supervised learning of physics-guided reconstruction neural networks without fully sampled reference data.  _Magnetic resonance in medicine_ , 84(6):3172–3191, 2020.

## Appendix

Figure A1: Early stopping for conv: Quantitative metrics are evaluated at each
training step. The maximum PSNR occurs at roughly $10^{5}$ steps before
decreasing slightly, consistent with the tendency of convolutional networks to
overfit to noise. While this decrease in PSNR is non-trivial, it is less
severe than the behavior of coord, discussed in Section 4.3.