\section{Introduction}
\noindent {\it Introduction.$-$}
The interaction between vector mesons and the nucleon ($N$) is one of the most fundamental quantities to study the meson properties in nuclear matter and the modification of QCD condensates at finite baryon density ~\cite{Hayano-Hatsuda2010}.
Among others, the $\phi$-meson attracts particular interest, since the interactions of the $s\bar{s}$ pair with the nucleon and nucleus are not well understood at low energies.

Experimentally, the production and absorption of the $\phi$-meson  in nuclear matter has been actively studied through photon and proton($p$)-induced reactions with various nuclear targets
(Spring8-LEPS~\cite{Ishikawa2005}, KEK-PS E325~\cite{KEK-PS2007}, JLab-CLAS~\cite{CLAS2010}, and COSY-ANKE~\cite{Polyanskiy2011}) as summarized in Ref.~\cite{Tolos2020}.
Recently, the ALICE Collaboration at LHC presented a new result on the $p$-$\phi$ interaction  through the femtoscopic analysis of the $p$-$\phi$ pair produced in the $pp$ reaction:
They reported a spin-averaged scattering length $a_{p\phi} = -0.85(34)_{\rm stat.}(14)_{\rm syst.}$ fm~\cite{ALICE2021}, an order of magnitude larger than that obtained from the photo-production data by the CLAS Collaboration at JLab~\cite{Dey:2014tfa} combined with the vector meson dominance \cite{Strakovsky2020}.

Theoretically, the properties of the $\phi$-meson in nuclear matter and its relation to the in-medium $\bar{s}s$ condensate have been extensively discussed~\cite{Hatsuda:1991ez,Asakawa:1994tp,Gubler:2016itj,Gubler:2018ctz}.
The microscopic origin of the force between $\phi$ and $N$ is also an important open problem.
In particular, the two-pion exchange (TPE) as the major long-range contribution between a color-dipole and the nucleon~\cite{TarrusCastella:2018php} is of strong relevance to this problem and will shed new light on the important role of pion dynamics in hadron interactions.

 Under these experimental and theoretical circumstances, it is most desirable to carry out realistic lattice QCD simulations of the $N$-$\phi$ interaction.
In this Letter, we report a first simulation of the $N$-$\phi$ system in a large lattice volume $\simeq$ (8.1 fm)$^3$  with light dynamical quarks near the physical point.
We focus on the highest-spin $N$-$\phi$ system, the $^4S_{3/2}$ channel,  with the notation $^{2s+1}L_J$ ($s$ = total spin, $L$ = orbital angular momentum, $J$ = total angular momentum). This is because  its coupling to two-body open channels $\Lambda K(^2D_{3/2})$ and $\Sigma K(^2D_{3/2})$ are kinematically suppressed at low energies  due to their $D$-wave nature, similar to the situation for the highest-spin   $N$-$\Omega$ system (the $^5S_{2}$ channel)  \cite{Iritani2019PLB}.
Moreover, decay processes with three or more particles in the final state such as $\Sigma \pi K$, $\Lambda \pi K$, and $\Lambda \pi\pi K$ are expected to be suppressed due to the small phase space.
Indeed, we will show below that visible effects from such decays cannot be seen from our lattice data in the $^4S_{3/2}$  channel.
By using the HAL QCD method~\cite{Ishii2007, Ishii2012, Aoki2020}, which converts the spacetime correlation of hadrons to the physical observables,
we obtain the $N$-$\phi$ scattering phase shift, the $S$-wave scattering length $a_{0}^{(3/2)}$ and the effective range $r_{\rm eff}^{(3/2)}$ in the $^4S_{3/2}$ channel. (For recent applications of the HAL QCD method to the baryon-baryon interactions, see Refs.~\cite{Gongyo2018, Iritani2019PLB, Sasaki2020, Lyu2021}.)

\noindent {\it HAL QCD Method.$-$}
Thanks to the Haag-Nishijima-Zimmermann reduction formula for composite particles~\cite{Zimmermann1987},
the interaction between hadrons can be extracted as the energy-independent non-local potential $U(\bm r,\bm r')$ ~\cite{Ishii2007}.
Such a potential can be defined from the normalized four-point function,
 \begin{eqnarray}\label{eq-R}
    R(\bm{r},t)&=&
   \frac{\sum_{\bm{x}}\langle 0| N(\bm{r}+\bm{x},t)\phi(\bm{x},t)\overline{\mathcal{J}}(0)|0\rangle}
    {\sqrt{Z_{\phi}Z_N}  e^{-(m_{N}+m_{\phi})t}} \nonumber \\
    &=&\sum_n a_n \psi_n(\bm{r})e^{-(\Delta E_n)t}+O(e^{-(\Delta E^*)t}).
\end{eqnarray}
Here the energy eigenvalue and the energy shift of each elastic scattering state are given by $E_n=\sqrt{m_N^2+\bm{k}_n^2}+\sqrt{m_\phi^2+\bm{k}_n^2}$
and   $\Delta E_n=E_n-(m_N+m_\phi)$, respectively, with $\bm{k}_n$ being the relative momentum  in the center of mass frame.
The equal-time Nambu-Bethe-Salpeter wave function for each eigenstate is denoted by $\psi_n(\bm r)$.
The wave function renormalization constant for $\phi$ ($N$) is  given by $Z_{\phi}$ ($Z_N$).
 The contributions from the inelastic scattering states are exponentially suppressed as $O(e^{-(\Delta E^*)t})$ with $\Delta E^*$ being the energy of the inelastic threshold relative to $m_N+m_{\phi}$.
The overlapping factor between the $N$-$\phi$   source operator  $\overline{\mathcal{J}}(0)$
 and the $n$-th eigenstate $|n\rangle$ is given by  $a_n=\langle{n}|\overline{\mathcal{J}}(0)|0\rangle$.
The wall-type  operator  with the Coulomb gauge fixing is used for the source, while the local operators are adopted for the sink,
\begin{equation}\label{eq-op}
 \begin{split}
  & N_\alpha(x)=\epsilon_{abc}\left[{{u^a}^T}(x)C\gamma_5d^b(x)\right]u^c_\alpha(x), \\
  & \phi_i(x)=\delta_{ab}\bar{s}^a(x)\gamma_is^b(x).
 \end{split}
\end{equation}
 Here $a$, $b$, and $c$ are the color indices, $i$ is the vector index, $\alpha$ is the Dirac index restricted to the upper two components, $C=\gamma_4\gamma_2$ is the charge conjugation, and $T$ denotes the transposition for the Dirac index.
Misner's method of approximate partial wave decomposition on a cubic grid, originally developed in the field of numerical relativity~\cite{Misner:1999ab}, is used to make the $S$-wave projection on the lattice~\cite{Miyamoto2020}.
The normalized four-point functions $R(\bm{r},t)$ with the spin projection onto the states with $J=3/2$ are used to extract the potential.

It has been shown that $R(\bm{r},t)$ satisfies the integro-differential equation ~\cite{Ishii2012},
\begin{equation}\label{eq-TDHAL}
 \begin{split}
  &\left[ \frac{1+3\delta^2}{8\mu}\frac{\partial^2}{\partial t^2}-\frac{\partial}{\partial t}-H_0 +O(\delta^2\partial_t^3) \right] R(\bm{r},t) \\
 &= \int d\bm{r}'U(\bm{r},\bm{r}')R(\bm{r}',t),
 \end{split}
\end{equation}
with $H_0=-\nabla^2/(2\mu)$, $\mu=m_N m_\phi/(m_N+m_\phi)$ being the reduced mass,  $\delta=(m_N-m_\phi)/(m_N+m_\phi)$ being the
 mass asymmetry.  We have neglected the term of $O(\delta^2\partial_t^3)$ in the above, since it is found to be consistent with zero within the statistical error in our simulation.  In practical calculations, a derivative expansion of the non-local potential is employed; $U(\bm{r},\bm{r}')=V(r)\delta(\bm{r}-\bm{r}')+ \sum\limits_{n=1}V_{n}(\bm{r})\nabla^{n}\delta(\bm{r}-\bm{r}')$.
 Then, the leading-order  central potential $V(r)$ reads,
\begin{equation}\label{eq-V}
 V(r)= R^{-1}(\bm{r},t)\left[ \frac{1+3\delta^2}{8\mu}\frac{\partial^2}{\partial t^2}-\frac{\partial}{\partial t}-H_0 \right] R(\bm{r},t).
\end{equation}
The truncation error of the derivative expansion is found to be small at low energies by using the method of finite volume spectra analysis~\cite{Iritani2019Jhep, Lyu2022} (Fig. S1 in the Supplemental Material).

\noindent{\it Lattice setup.$-$}
The ($2+1$)-flavor gauge configurations are generated with the Iwasaki gauge action at $\beta=1.82$ and the nonperturbatively $O(a)$-improved Wilson quark action with stout smearing at nearly physical quark masses~\cite{Ishikawa2016}.
The lattice spacing is $a\simeq0.0846$ fm ($a^{-1}\simeq2333$ MeV), and the lattice volume is $L^4=96^4$, corresponding to $La\simeq8.1$ fm.
The isospin-averaged masses of $\pi$, $K$, $\phi$, and $N$
 from the single-state fitting  in the interval $t/a=15$-25 for the mesons and 11-18 for the nucleon are
 listed in Table ~\ref{tab-mass}  with statistical errors in the parentheses.
Corresponding experimental values~\cite{PDG2020} are shown in the third column.
The lattice results are about 6\%  larger than the experimental values for the pseudo-scalar mesons, while about 3\% (2\%) larger for $\phi$ ($N$).
 With our light quark masses slightly heavier than the physical masses, $\phi$ is below the $K\bar{K}$ threshold,
which prohibits a fall-apart decay of  $\phi$.

\begin{table}[htbp]
\caption{ Isospin-averaged hadron masses with statistical errors obtained from (2+1)-flavor lattice QCD simulations together with the experimental values.
}
\begin{tabular}{ccc}
  \hline\hline
    Hadron &~~Lattice [MeV]  &~~Expt. [MeV] \\
  \hline
    $\pi$ &\ \ 146.4(4)  & \ \ 138.0\\
    $K$  & \ \ 524.7(2) & \ \ 495.6\\
    $\phi$ &\ \ 1048.0(4) &\ \ 1019.5\\
    $N$  &\ \ 954.0(2.9)  & \ \ 938.9 \\
  \hline\hline
\end{tabular} \label{tab-mass}
\end{table}

We use 200 configurations separated by 10 trajectories.
For each configuration, the forward and backward propagations are averaged, the hypercubic symmetry on the lattice (4 rotations) is utilized, and 80 measurements are performed by shifting the source position in a temporal direction in order to reduce the statistical fluctuation.
In total, the number of measurements is 128000.
With the  periodic boundary condition for all directions, the quark propagators are calculated by the domain-decomposed solver~\cite{Ishikawa2021},
and the hadron correlation functions are obtained by the unified contraction algorithm~\cite{Doi2013}.
The contraction diagram with the pair annihilation of $s\bar{s}$  is not considered,
since it is known to be small, e.g. less than 0.4\% correction to the
 $\phi$-meson mass \cite{Chakraborty:2017hry}.
The statistical errors are evaluated by the jackknife method with a bin size of 20 configurations throughout this Letter.
A comparison with a bin size of 40 configurations shows that the bin size dependence is small.
\begin{figure}[htbp]
  \centering
  \includegraphics[width=8cm]{Fig1}
  \caption{(Color online).
  The $N$-$\phi$ potential $V(r)$ in the ${^4S_{3/2}}$ channel as a function of separation $r$ at
 Euclidean time $t/a= 12$ (red squares), $13$ (green circles) and $14$ (blue triangles).
  } \label{Fig1}
\end{figure}

\noindent{\it Numerical results.$-$}
The $N$-$\phi$ potential $V(r)$ in the $^4S_{3/2}$ channel defined in Eq.~(\ref{eq-V}) with the
lattice measurement of $R(\bm r, t)$ is shown in Fig.~\ref{Fig1} for Euclidean times, $t/a=12$, $13$, and $14$. (See Fig. S2 in the Supplemental Material for the $t$-dependence of $V(r)$ in a wider range of $t$.)
These Euclidean times are chosen such that they are large enough to suppress contaminations from excited states in the single-hadron correlator and simultaneously small enough to avoid exponentially
increasing statistical errors.
The potentials from different $t$ agree with each other within statistical errors, which indicates that the $N$-$\phi$ correlation function is dominated by the elastic scattering states in the $^4S_{3/2}$ channel; we do not find the effects of the two-body open channels  ($\Lambda K(^2D_{3/2})$ and $\Sigma K(^2D_{3/2})$) as well as the  three-body open channels including $N\phi \rightarrow \{\Sigma^* K, \Lambda(1405) K\} \rightarrow \{\Lambda \pi K, \Sigma \pi K$\}.
This is in sharp contrast to the  $^2S_{1/2}$ case where we found that the $N$-$\phi$ potential shows a clear $t$-dependence, as expected from the $S$-wave fall-apart decay into $\Lambda K(^2S_{1/2})$ and $\Sigma K(^2S_{1/2})$.

The potential $V(r)$ in the $^4S_{3/2}$ channel shown in Fig.~\ref{Fig1} is attractive for all distances and has a characteristic two-component structure, the
attractive core at short distance and the attractive tail at long distance,  similar to the case of the $N\Omega(^5S_{2})$ potential~\cite{Iritani2019PLB}.  We note that the Pauli exclusion principle between quarks, which partially gives rise to the repulsive core in the $NN$ interaction~\cite{Inoue2012,Oka2000},  does not operate in the present case, since $N$ and $\phi$ have no common valence quarks.

As has been discussed  for the interaction between color dipoles~\cite{Bhanot:1979vb,Fujii:1999xn,Brambilla:2015rqa},
  non-perturbative gluon exchange is  expected to appear in the form of  the TPE at long distance.
  The idea was generalized to the interaction between a color-dipole and the nucleon with the result,
 $V(r \gg (2m_{\pi})^{-1} ) = -\alpha  \frac{\exp (-2m_\pi r)}{r^{2}}$, where $\alpha$ is proportional to $m_{\pi}^4$ \cite{TarrusCastella:2018php}.
To check such a long distance behavior of $V(r)$,  we
 show in Fig.~\ref{Fig2} the spatial effective energy, $E_{\rm eff}(r)=-\frac{\ln[-V(r)r^{2}/\alpha]}{r}$,
as a function of $r$ with $\alpha \simeq 91~{\rm MeV}\cdot{\rm fm}^{2}$ determined by fitting the lattice data of $V(r)$ at long distance.
We find that $E_{\rm eff}(r)$ has a plateau at  $2m_\pi = 292.8$ MeV  for $r > 1.0$ fm, which indicates that the long-range part of the $N$-$\phi$ potential is indeed dominated by the TPE.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=8cm]{Fig2}
  \caption{(Color online).
  The spatial effective energy $E_{\rm eff}(r)$ as a function of separation $r$ at
 Euclidean time $t/a= 12$ (red squares), $13$ (green circles) and $14$ (blue triangles).
The orange dashed line corresponds to $2m_\pi$ with lattice pion mass $m_\pi=146.4$ MeV.
  } \label{Fig2}
\end{figure}

In order to convert the potential to physical observables, we perform an uncorrelated fit of the lattice QCD potential
by using two different functional forms,
\begin{eqnarray}\label{eq-fit}
 \! \! {\rm{A}}: V_{\rm{fit}}(r) &=& \sum_{i=1,2}a_i e^{-(r/b_i)^2} + a_3 m_\pi^4 f(r;b_3)  \frac{e^{-2m_\pi r}}{r^{2}},  \nonumber \\
\! \! {\rm{B}}:  V_{\rm{fit}}(r)&=&\sum_{i=1,2,3}a_i e^{-(r/b_i)^2}.
\end{eqnarray}

The fit A is motivated by the TPE tail at long distance with an overall strength proportional to $m_{\pi}^4$~\cite{TarrusCastella:2018php},
  while the fit B is  a purely  phenomenological Gaussian form  for comparison.
  In fit A, we consider two types of  form factors  commonly used in the $NN$ potentials: (i) The Nijmegen-type form factor \cite{Stoks:1994wp},
  $f_{\rm erfc}(r;b_3)= \left[  {\rm erfc}\left( \frac{m_{\pi}}{\Lambda}-\frac{\Lambda r}{2}  \right) - e^{2m_{\pi}r} {\rm erfc}\left(  \frac{m_{\pi}}{\Lambda}+\frac{\Lambda r}{2}   \right)  \right]^2/4$
with $\Lambda = 2/b_3$ and the complementary error function ${\rm erfc}(x)= \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-z^2} dz$.
This is motivated by the exponential-type regularization of the pion propagator in the momentum space, $1/(k^2+m_{\pi}^2)  \rightarrow e^{-(k/\Lambda)^2 } /(k^2+m_{\pi}^2 )$.
 (ii) The Argonne-type  form factor~\cite{Wiringa1995}, $ f_{\rm exp}(r;b_3) =  \left(1-e^{-(r/b_3)^{2}}\right)^2$.
We refer to fit A with (i) and (ii) as fit A$_{\rm erfc}$ and fit A$_{\rm exp}$, respectively.
The pion mass in fit A$_{\rm erfc, exp}$ is taken to be  $m_\pi=146.4$ MeV, and the fit range is  chosen as  $0 < r < 3.0$ fm.
We found that the fit A$_{\rm erfc, exp}$ and the fit B provide an equally good fit
 ($\chi^2_{\rm d.o.f}=$0.3-0.4) and are stable against the choice of $t$.
 In Table~\ref{tab-prm}, we show the fit results for $t/a=14$, which are expected to have least contamination from the inelastic states.
Changing the fit range of the potential to  $0.1<r<2.5$ fm does not affect the results within statistical errors.
Also we found that the simple fitting functions  such as the
Yukawa  form $\sim  - \frac{\exp (-\mu r)}{r}$ ~\cite{Brodsky1990,Gao2001} and  the van der Waals (Casimir-Polder) form $\sim - \frac{1}{r^k}$ with $k=6$ ($7$)~\cite{Appelquist1978}  cannot reproduce the lattice data.

\begin{table}[hbtp]
\caption{The fit parameters in Eq.~(\ref{eq-fit}) with statistical errors quoted in the parentheses at $t/a=14$.
The fit range is $0<r<3.0$ fm. In $a_3 m_\pi^{4n}$, we take  $n=1$ and $n=0$ for fit A and B, respectively.
A$_{\rm erfc}$ (A$_{\rm exp}$) denotes fit A with the Nijmegen-type (Argonne-type) form factor.}
\begin{tabular}{cccc}
  \hline\hline
    fit~~&A$_{\rm erfc}$~~~&A$_{\rm exp}$~~~&B\\
  \hline
    $a_1$ [MeV] ~~~&-376(20)~~~&-371(27)~~~&-371(19)\\
    $b_1$ [fm] ~~~&0.14(1)~~~&0.13(1)~~~&0.15(3)\\
    $a_2$ [MeV] ~~~&306(122)~~~&-119(39)~~~&-50(35) \\
    $b_2$ [fm] ~~~&0.46(4)~~~&0.30(5)~~~ &0.66(61)\\
    $a_3 m_\pi^{4n}$ [MeV$\cdot{\rm fm}^{2n}$]~~~&-95(13)~~~&-97(14)~~~&-31(53)\\
    $b_3$ [fm] ~~&0.41(7)~~~&0.63(4)~~~&1.09(41) \\
  \hline\hline
\end{tabular} \label{tab-prm}
\end{table}

\begin{figure}[hbtp]
  \centering
  \includegraphics[width=8cm]{Fig3}
  \caption{
    The $N$-$\phi$ scattering phase shifts $\delta^{(3/2)}_0$ in the $^4S_{3/2}$ channel obtained from $V_{\rm fit}(r)$ at $t/a=12$  (red squares),  $13$ (green circles), and $14$ (blue triangles).
  } \label{Fig3}
\end{figure}

Figure~\ref{Fig3} shows the $N$-$\phi$ scattering phase shifts $\delta^{(3/2)}_0$ in the $^4S_{3/2}$ channel as a function of the center of mass kinetic energy $E_{\mathrm{CM}}=\sqrt{m_N^2+{k}^2}+\sqrt{m_\phi^2+{k}^2}-(m_N+m_\phi)$ obtained by using $V_{\rm fit}(r)$ with the fit A$_{\rm erfc}$.
The scattering phase shifts from different $t$ are consistent with each other within statistical errors.
 The scattering length $a^{(3/2)}_0$ and the effective range $r^{(3/2)}_{\rm eff}$ can be extracted from the
 effective range expansion for small $k$ as $k\cot\delta^{(3/2)}_0(k)=-\frac{1}{a^{(3/2)}_0}+\frac12 r^{(3/2)}_{\rm eff}k^2 +O(k^4)$.
 In Table~\ref{tab-scattering}, $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ are shown for the present  pion mass $m_{\pi}=146.4$ MeV;
 the central values and the statistical errors are obtained from the data at $t/a=14$ with the  fit A$_{\rm erfc}$,  while
  the systematic errors in the second parentheses are estimated by comparing results for $t/a=12$-14 with A$_{\rm erfc, exp}$ and B.
Other possible systematic errors are as follows; (i) the finite volume effect of the potential, which is expected to be small as  $\exp (-2m_{\pi} (L/2)a) \lesssim 0.3\%$ due to the large volume, and (ii) the finite cutoff effect, which is also expected to be small as $O((a\Lambda_{\rm QCD})^2)\sim O(1)\%$ due to the non-perturbative $O(a)$-improvement.
As an alternative estimate, we study how much scattering parameters are sensitive to the potential in short range, where the finite cutoff effect should be most prominent.
It is found that $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ change only $\sim 2\%$ even when we completely cut the short-range potential at $r < 0.1$ fm.

\begin{table}[tbhp]
\caption{The scattering length $a^{(3/2)}_0$ and the effective range $r^{(3/2)}_{\rm eff}$ obtained by using $V_{\rm fit}(r)$ at $m_\pi=146.4$ MeV
 with statistical and systematic errors. Estimated central values using  a model-dependent extrapolation  of $V_{\rm fit}(r)$ to $m_\pi=138.0$ MeV are
  also shown for comparison.}
\begin{tabular}{ccc}
  \hline\hline
    $m_\pi$ [MeV]~&$a^{(3/2)}_0$ [fm]~~~&$r^{(3/2)}_{\rm eff}$ [fm]\\
  \hline
    146.4 ~~~&$-1.43(23)_{\rm stat.}\left(^{+36}_{-06}\right)_{\rm syst.}$~~~&$2.36(10)_{\rm stat.}\left(^{+02}_{-48}\right)_{\rm syst.}$\\
      138.0 ~~~&$\simeq -1.25$~~~&$\simeq 2.49$\\

 \hline\hline
\end{tabular} \label{tab-scattering}
\end{table}

To estimate how the scattering parameters change toward the physical quark mass,
we keep $a_{1,2,3}$ and $b_{1,2,3}$ in $V_{\rm fit}(r)$ fixed in fit A$_{\rm erfc, exp}$ and smoothly change the long-range potential by taking the isospin-averaged physical pion mass $m_{\pi} =138.0$ MeV  in the region where the TPE is dominated ($r>1.0$ fm from Fig.\ref{Fig2}).
By calculating the scattering phase shifts with such a potential with the physical masses of  $\phi$ and $N$, we obtain estimated values of  $a^{(3/2)}_0$ and $r^{(3/2)}_{\rm eff}$ for  $m_{\pi} =138.0$ MeV in Table~\ref{tab-scattering}.
 Although the range of the TPE is increased by the smaller pion mass, the characteristic $m_{\pi}^4$ behavior  of the TPE strength makes the overall attraction weaker.
Note that this is only a model-dependent qualitative estimate and needs to be confirmed by future physical-point simulations.

\begin{figure}[tbhp]
  \centering
  \includegraphics[width=8cm]{Fig4}
  \caption{
    The scattering length $a^{(3/2)}_0$ obtained from $V(r;r_{\rm c})=\theta (r_{\rm c}-r)  V_{\rm fit}(r)$ as a function of  cutoff length $r_{\rm c}$ at $t/a=12$  (red squares),  $13$ (green circles), and $14$ (blue triangles).
  } \label{Fig4}
\end{figure}

Since we do not have reliable information on the $N\phi(^2S_{1/2})$ potential from lattice QCD at the moment due to the effect of the open channels,
a comparison of our results with spin-averaged scattering parameters should be made with caution. With this reservation in mind, our $a^{(3/2)}_0$ is found to be one or two orders of magnitude larger than the previous theoretical results in QCD sum rules~\cite{Koike1997,Klingl1997kf}.
Such a discrepancy may be due to the difficulty of obtaining the long-range TPE contribution from the low-order truncation of the operator product expansion in QCD sum rules.  In fact,  the magnitude of $a^{(3/2)}_0$ becomes considerably smaller when the long-range potential is cut off.
Shown in Fig.~\ref{Fig4} is $a^{(3/2)}_0$ as a function of cutoff length $r_{\rm c}$ obtained by the potential,
$  V(r;r_{\rm c}) = \theta (r_{\rm c}-r)  V_{\rm fit}(r)$ with the fit A$_{\rm erfc}$.
Considerable decrease of  $|a^{(3/2)}_0|$ from $1.43$ fm at $r_{\rm c}=\infty$ to about $0.1$ fm at $r_{\rm c}=0.5$ fm can be seen.

\noindent{\it Summary.$-$}
In this Letter, we present a first lattice QCD calculation on the interaction of the  $N$-$\phi$ system in the $^4S_{3/2}$
channel based on the  ($2+1$)-flavor simulations with nearly physical quark masses.
The interaction potential in the  $N\phi(^4S_{3/2})$ channel is extracted from lattice data of the
hadronic spacetime correlation using the HAL QCD method.
The potential is found to be attractive for all distances and appears to be a combination of an attractive core at short distances and a two-pion exchange (TPE) tail at long distances ($r >$ 1 fm).
The latter is well fitted by the characteristic form of the TPE obtained by the interaction of a color-dipole and the nucleon.
The scattering parameters obtained from our potential at $m_{\pi}=146.4$ MeV is summarized in Table \ref{tab-scattering}.
By examining the potential fitted to the lattice data, we find that the scattering length $a^{(3/2)}_0$ is sensitive to the length scale of $r > 0.5$ fm.
Also, we suggest that the $N$-$\phi$ attraction could be weaker at the physical pion mass due to the characteristic $m^4_{\pi}$ dependence of the strength of the TPE.

Our $a^{(3/2)}_0$ is substantially larger in magnitude than the previous calculations of the spin-averaged $a_0$ using QCD sum rules but is comparable to the spin-averaged $a_0$  by ALICE Collaboration within the error bar ~\cite{ALICE2021}.
Also, our $r^{(3/2)}_{\rm eff}$ is about three times smaller than the spin-averaged $r_{\rm eff}$ by ALICE Collaboration.
To make a firm comparison between theory and experiment, we need to perform a $\Lambda K$-$\Sigma K$-$N\phi$ coupled-channel analysis to reveal the complex dynamics of the $^2S_{1/2}$ and $^4S_{3/2}$
channels as well as to perform physical-point simulations.

The present lattice QCD study near the physical point provides a first step to exploring the interaction of $s\bar{s}$ with the nucleon
 from first principle.  The heavier system such as $c\bar{c}$ interacting with the nucleon pioneered in ~\cite{Kawanai:2010ev,Sugiura:2019pye}
  is also an important problem to be studied near the physical point.

\begin{acknowledgments}
{\it Acknowledgments.$-$}
We thank members of the HAL QCD Collaboration for stimulating discussions.
Y. L. thanks Xu Feng for the helpful discussions.  We thank members of the PACS Collaboration for the gauge configuration generation conducted on the K computer at RIKEN.
The lattice QCD measurements have been performed on Fugaku and HOKUSAI supercomputers at RIKEN.
We thank ILDG/JLDG \cite{ldg}, which serves as essential infrastructure in this study.
This work was partially supported by HPCI System Research Project (hp120281, hp130023, hp140209, hp150223, hp150262, hp160211, hp170230, hp170170, hp180117, hp190103, hp200130, hp210165, and hp210212),  the National Key R\&D Program of China (Contract Nos. 2017YFE0116700 and 2018YFA0404400), the National Natural Science Foundation of China (Grant Nos. 11935003, 11975031, 11875075, and 12070131001), the JSPS (Grant Nos. JP18H05236, JP16H03978, JP19K03879, JP18H05407, and JP21K03555), the MOST-RIKEN Joint Project ``Ab initio investigation in nuclear physics'', ``Priority Issue on Post-K computer'' (Elucidation of the Fundamental Laws and Evolution of the Universe), ``Program for Promoting Researches on the Supercomputer Fugaku'' (Simulation for basic science: from fundamental laws of particles to creation of nuclei), and Joint Institute for Computational Fundamental Science (JICFuS).

\end{acknowledgments}

\


\section{Introduction}

\IEEEPARstart{I}{schemic} heart disease is one of the leading causes of death in western countries \cite{kishore2011global}.
It refers to ischemia of the myocardium due to stenosis of the coronary arteries, which supply oxygenated blood to the myocardium.
If a stenosis develops to completely occlude the vessel, the patient undergoes a myocardial \emph{infarction}, i.e., heart attack.
The part of the myocardium which undergoes a prolonged shortage of oxygen is damaged and can be either non-viable or hibernating after the infarction.
Since the hibernating myocardium has the potential to resume contraction after re-vascularisation \cite{hibernatingMyocardium1989}, viability assessment of the myocardium is essential for diagnosis and therapy planning.

Among various cardiac scans, late gadolinium enhanced (LGE) cardiac magnetic resonance (CMR) imaging can directly visualize and thus discriminate non-viable myocardium (known as infarcts) from normal myocardium \cite{LGEmedical1999}.
In an LGE CMR scan, a gadolinium-based contrast agent is injected and a single-frame sequence is acquired 10 to 20 minutes later, by which time infarcts exhibit hyper-enhanced intensities compared to healthy myocardium.
This phenomenon has been hypothesized to be the result of delayed wash-out kinetics of the contrast agent in non-viable myocardium.
Post-mortem histologic staining of the myocardium using animal models has shown that the hyper-enhanced regions in LGE CMR images correlate well with the location and extent of non-viable tissue \cite{LGEmedical1999,fieno2000contrast}.
However, there is one exception called \emph{microvascular obstruction} (MVO).
MVO refers to the no-reflow phenomenon in which the involved sub-endocardial infarcts do not get enhanced and appear as dark as normal myocardium.
It is mostly observed in acute infarctions~\cite{MVO2010}.
An LGE image with a hyper-enhanced infarct and MVO is shown in Fig.~\ref{fig:LGEillustration4Introduction}.
LGE CMR is usually triggered with electrocardiography (ECG) gating at mid-diastole of the cardiac cycle.

\begin{figure}
    \centerline{
    \hfil
    \subfloat[]{
        \includegraphics[width=.3\columnwidth]{LGEillustration4Introduction_a}
        \label{fig_first_case}
    }
    \hfil
    \subfloat[]{
        \includegraphics[width=.3\columnwidth]{LGEillustration4Introduction_b}
        \label{fig_second_case}
    }
    \hfil
    \subfloat[]{
        \includegraphics[width=.3\columnwidth]{LGEillustration4Introduction_c}
        \label{fig_second_case}
    }
    \hfil
    }\caption{An example LGE image: (a) the original image; (b) the hyper-enhanced infarct is outlined. (c) the MVO is outlined.}\label{fig:LGEillustration4Introduction}
\end{figure}

Quantification of LGE CMR images can be divided into two stages: segmentation of the myocardium and classification of infarcts within the segmented myocardium.
Although the quantification can be done manually by experts, it is not only time-consuming but also subject to inter- and intra-observer variation.
Therefore, computer aided \mbox{(semi-)} automatic techniques are highly desired.
However, the automation is not straightforward.
First, automatic segmentation of the myocardium is often difficult due to the intensity heterogeneity of the myocardium and intensity similarity between the infarcts and blood pool (BP).
Second, misclassification of infarcts can happen because of the intensity inconsistencies and misalignment artifacts of a stack of slices, image noise and other artifacts.

Our contribution in this paper is a complete and comprehensive 3D framework for automatic quantification of LGE CMR images, which is rare in the literature.
This framework includes two crucial pre-processing measures -- corrections of misalignment artifacts and intensity inconsistencies -- that are omitted or improperly handled in previous works.
Furthermore, both segmentation and classification stages in the framework are 3D in nature, hence utilizing intensity and spatial continuity beyond individual 2D images.
The rest of this paper is organised as follows.
Section \ref{sec:background} introduces the background knowledge and reviews related works.
Section \ref{sec:method} describes the methods used in every step of our 3D framework.
The experimental results are presented and discussed in Section \ref{sec:results}, followed by the conclusions in Section \ref{sec:conclusion}.

\section{Background and Related Works}
\label{sec:background}

\subsection{Misalignment Correction}

Misalignment artifacts are common for multi-slice CMR data.
Clinically, a few consecutive short-axis (SA) slices and 1-3 standard long-axis (LA) slices are acquired.
Slices are usually acquired in multiple breath-holds, introducing differences in lung volume which in turn cause displacement of the heart \cite{heartMotionStudy}.
Misalignment artifacts can also be caused by patient motion and heartbeat;
the former can be considered similar to the breath-related artifacts, while the latter is often neglected due to the use of ECG gating.
Misalignment artifacts cause slices to move away from their true positions relative to each other (Fig.~\ref{fig:misalignment_n_correction_n_intersection}, first row) and, if not tackled, may lead to distortion of the 3D reconstruction from the multi-slice data.

\begin{figure}
  \centering
  \includegraphics[width=.95\columnwidth]{misalignment_n_correction_n_intersection}\\
  \caption{First row: misalignment artifacts illustrated with intersections between SA and LA slices.
  Second row: the same slices after misalignment correction.
  The right column shows sampled line segments along the intersection line from the slices in the mid column.
  4C: 4-chamber; 2C: 2-chamber.}\label{fig:misalignment_n_correction_n_intersection}
\end{figure}

Researchers have proposed various approaches to realign a set of CMR slices.
\citet{A_Comprehensive_Approach} registered LGE slices to a whole-heart coronary angiography volume slice by slice by maximizing similarity measures in the overlapping image parts.
However, the high resolution whole-heart volume data is not widely available under clinical settings.
\citet{alignManyLA2004} proposed a method more suitable for clinically acquired CMR, where slices were realigned by maximizing normalized mutual information (NMI) of the intersecting parts between slices.
However, no regularization term was included to avoid unrealistic movements.
Recently, \citet{alignEqualSpacing2010} introduced an effective regularization term favoring original spacing among SA slices, but its underlying assumption that the original spacing be correct is sometimes broken by superior-inferior motion of the heart.
A common deficiency of \cite{alignEqualSpacing2010,alignManyLA2004} is that they merely used similarity at the intersections between slices but ignored the anatomical continuity of the heart within the stack of SA slices.
Furthermore, both \cite{alignManyLA2004} and \cite{alignEqualSpacing2010} worked with multi-frame slices.
As LGE data is single-frame, the available information that can be used for its realignment is $N_\mathrm{f}-1$ times fewer, where $N_\mathrm{f}$ is the number of available frames.

\subsection{Myocardium Segmentation}

The delineation of myocardial contours is a prerequisite to subsequent automatic localization and quantification of infarcts -- nearly all such works~\cite{infarctSeg2005,A_Comprehensive_Approach,infarctSeg2010a,infarctSeg2010b,infarctSeg2011} in the literature assume that high quality myocardium segmentation is given, either manually or (semi-) automatically.
However, automatic segmentation of the myocardium in LGE images is a difficult task.
There has been little research aimed at automatic myocardium segmentation for LGE CMR images, and there is no commercially or publicly available automatic segmentation tool for clinical use.
Most existing approaches utilize pre-delineated myocardial contours in cine CMR data of the same patient as \emph{a priori} knowledge~\cite{relatedworks01_Ciofolo2008,relatedworks02_Dikici2004,myownMICCAI2011}.
Such an approach is reasonable because the patient is asked to stay still during the entire acquisition process and there are many methods available for (semi-)automatic segmentation of the cine data, e.g., \cite{cineSegmentation_chen2008}.
The major difficulties of this approach include:
(i) displacement and nonrigid deformation between the cine and LGE data due to respiratory and body motion and imperfectness of ECG gating; and
(ii) different intensity characteristics of the cine and LGE data.

\subsection{Infarct Classification}

Early works on infarct classification were purely intensity-based with a global threshold, e.g., \cite{infarctSeg2005,FWHM2004}.
These methods suffer from the lack of spatial continuity, resulting in separated false positives probably caused by noise or artifacts.
To overcome this problem, \citet{FACT_I2006} proposed the feature analysis and combined thresholding (FACT) method, which enforces spatial constraints by assuming that the infarcts are generally sub-endocardial and reasonably large.
Recent works applied more advanced techniques to this problem.
\citet{heiberg2005semi} proposed a level set algorithm to regularize the thresholding and exclude small regions that constitute noise rather than infarction.
Similarly, \citet{metwally2010improved} proposed to enhance the thresholding using k-means clustering.
In \cite{A_Comprehensive_Approach} and \cite{infarctSeg2010a}, advanced signal intensity analysis is carried out based on the Rician distribution of noisy MRI data \cite{ricianMRI1995}, where the myocardium is modeled with a mixture model comprising a Rayleigh distribution for normal myocardium plus a Gaussian for infarcts.
However, these methods are generally 2D and hence do not take advantage of 3D spatial and intensity continuity.

Recently, \citet{infarctSeg2010b} proposed a method which first globally thresholds the image followed by false-positive and false-negative removal combining intensity and spatial information.
The method has three novelties:
(i) The LV BP is included as the same class as the infarcts during classification.
(ii) All SA slices of a patient are used together in the intensity distribution analysis to derive the global threshold.
This strategy works effectively when there is no obvious intensity inconsistency in a set of LGE images but \emph{degenerates} in the presence of significant intensity inconsistency.
(iii) The Otsu threshold \cite{OtsuThreshold1979} which maximizes the ratio of between-class and within-class variances is adopted, so there is no need to make any assumption regarding the underlying true distribution.
Nevertheless, this work suffers from two drawbacks: there is neither correction for intensity inconsistencies nor correction for misalignment artifacts across slices, and it is still based on a global threshold.

The only work in the literature that handled both misalignment artifacts and intensity inconsistencies was presented in \cite{infarctSeg2011}.
SA slices are realigned by aligning myocardial centers of all slices in a line;
however, this is not always correct since anatomically the slices' myocardial centers do not necessarily form a line.
The intensity of each slice is adjusted to the same range as in the first slice.
Although doing so can make the overall intensity characteristics of the entire images more consistent, it does not necessarily normalize the intensities of the LV, which despite being the real target only occupies a relatively small area in each image.
Another drawback of this work is that again global thresholding is performed as the major classification process.

\subsection{Comprehensive Frameworks}

To the best of our knowledge, the only existing comprehensive framework for automatic quantification of LGE CMR images is \cite{A_Comprehensive_Approach}, which conducted a joint analysis of the LGE and perfusion CMR images.
However, its segmentation of the myocardium in LGE images was done semi-interactively.
Although a solution for the misalignment artifacts of LGE slices was proposed, the navigator-gated whole-heart coronary angiography volume data, which was used as the reference for the realignment of LGE slices, is not widely available in usual clinical settings.
Finally, this work did not consider intensity inconsistencies across LGE slices and it is unclear whether the detection of infarcts was done in 2D or 3D.

\section{Method}
\label{sec:method}
The proposed framework comprises four steps:
(i) correction of misalignment artifacts;
(ii) 3D segmentation of the myocardium;
(iii) intensity normalization of the LV regions in a stack of SA slices; and
(iv) 3D classification of the myocardium into infarcted and normal regions.

\subsection{Correction of Misalignment Artifacts}

We propose an effective and robust method for misalignment correction of LGE CMR slices.
Unlike previous methods which merely rely on similarity measurement at intersections of slices, we propose to also utilize the innate physical continuity of the heart throughout the SA stack to establish interlinks within the SA slices and group them together.
This is achieved by introducing a novel cost term, namely, the contiguous cost, between adjacent SA slices.
The contiguous cost not only prevents unrealistic movement of the slices, but also improves accuracy of the misalignment correction.

In our method, we incorporate only translational but no rotational correction based on two research findings:
(i) a study on motion and deformation of the heart due to respiration found only very small rotations (typically a couple of degrees) \cite{heartMotionStudy}.
(ii) in another study \cite{alignEqualSpacing2010}, a method capable of full 3D translational and rotational corrections found no significant improvement on final results with the rotational correction.

\subsubsection{Intersecting Cost}
The intersecting cost refers to the dissimilarity measure along the intersecting line between a pair of intersecting slices.
Assume there are $m$ SA and $n$ LA slices (denoted by $I^{k}_{\mathrm{SA}}$ and $I^{j}_{\mathrm{LA}}$, where $k\in[1,m]$ and $j\in[1,n]$).
The position and orientation of a slice in the patient coordinate system are uniquely defined by two fields of the standard DICOM header -- ImagePositionPatient (IPP) and ImageOrientationPatient (IOP), respectively.
Given IPP and IOP of two slices, their intersection line can be computed.
Then, two line segments are sampled along the intersection line on both slices (Fig.~\ref{fig:misalignment_n_correction_n_intersection}, third column).
The smaller PixelSpacing (PS, from the DICOM header) of the two slices' is used as the step size during sampling with linear interpolation.
Finally, the intersecting cost between the two slices is defined as the dissimilarity measure of the two sampled line segments:
\begin{equation}\label{eq:intersectingErr}
    E_{\mathrm{int}}=\begin{cases}
        &E_{\mathrm{int}}(I^{k}_{\mathrm{SA}}\bigcap I^{j}_{\mathrm{LA}}) = \mathcal{C}_{\boldsymbol{s}}(\boldsymbol{s}^{k}_{\mathrm{SA}}, \boldsymbol{s}^{j}_{\mathrm{LA}}),\\
        &E_{\mathrm{int}}(I^{j}_{\mathrm{LA}}\bigcap I^{j'}_{\mathrm{LA}}) = \mathcal{C}_{\boldsymbol{s}}(\boldsymbol{s}^{j}_{\mathrm{LA}}, \boldsymbol{s}^{j'}_{\mathrm{LA}}),
    \end{cases}
\end{equation}
where $I^{k}_{\mathrm{SA}}\bigcap I^{j}_{\mathrm{LA}}$ denotes the intersection between the $k^\mathrm{th}$ SA and $j^\mathrm{th}$ LA slices, and $I^{j}_{\mathrm{LA}}\bigcap I^{j'}_{\mathrm{LA}}$ the intersection between the $j^\mathrm{th}$ and ${j'}^\mathrm{th}$ LA slices ($j\neq j'$), $\boldsymbol{s}^{k}_{\mathrm{SA}}$ and $\boldsymbol{s}^{j}_{\mathrm{LA}}$ the sampled line segments on SA and LA slices respectively, and $\mathcal{C}_{\boldsymbol{s}}$ a function measuring dissimilarity between two line segments.

$\boldsymbol{s}^{k}_{\mathrm{SA}}$ and $\boldsymbol{s}^{j}_{\mathrm{LA}}$ depend on IPP, IOP and PS of the two intersecting slices.
Of the three, IOP and PS can be considered constants as neither is changed by our method.
Consequently, $E_{\mathrm{int}}$ is actually a function of IPP.
Particularly, IPP of a slice stores the $x$, $y$ and $z$ coordinates of the slice in the patient coordinate system.
Thus, translation of a slice can be achieved conveniently by changing its IPP.
Because we are only interested in the LV, only a fraction of the intersection line between SA and LA slices is sampled to form $\boldsymbol{s}^{k}_{\mathrm{SA}}$ and $\boldsymbol{s}^{j}_{\mathrm{LA}}$.
The sampling range is confined within a region of interest (ROI) on each SA slice, which is predefined by identifying a bounding box that well contains the LV (Fig.~\ref{fig:ROI}(a)).
However, for the intersection line between different LA slices, the entire length is sampled to form $\boldsymbol{s}^{j}_{\mathrm{LA}}$ and $\boldsymbol{s}^{j'}_{\mathrm{LA}}$.
This is because if only the portion within the LV is sampled, there would be too little discriminative information (mostly from the BP) to determine relative positions of the LA slices.

\begin{figure*}
  \centering
  \includegraphics[width=.9\textwidth]{ROI}\\
  \caption{Illustration of the sampling ranges for
  (a) the intersecting cost: $\boldsymbol{s}^{k}_{\mathrm{SA}}$ and $\boldsymbol{s}^{j}_{\mathrm{LA}}$ are sampled along the intersection (dashed white) line only from the portion (yellow segment) lying within the SA slice's ROI (green rectangle);
  (b) the contiguous cost: the regions (white rectangles) in which $R^{k}_{\mathrm{SA}}$ and $R^{k+1}_{\mathrm{SA}}$ are sampled are defined by finding the smallest cuboid whose top and bottom faces are in the involved SA slices (planes) and can completely contain the corresponding ROI (dashed blue and green) in each slice.
  Distance between contiguous SA slices is exaggerated here for better visualization.}\label{fig:ROI}
\end{figure*}

\subsubsection{Contiguous Cost}
The contiguous cost refers to the discontinuity measure between two adjacent SA slices.
The LV is naturally smooth, and hence its representation throughout the stack of SA slices should be continuous\footnote{
Although there could sometimes be a gap between adjacent SA slices in clinical CMR data, this gap is always small
enough to obtain a correct 3D representation of the LV, and thus the smoothness of the LV is retained.}.
Although the continuity of the LV throughout the SA stack is an important physical property, it has been ignored by previous works.
We define the discontinuity measure between contiguous SA slices as the contiguous cost (as opposed to the intersecting  cost).
The contiguous cost is combined with the already defined intersecting cost to form the total cost to be minimized.
On one hand, the contiguous cost acts as a regularization to the intersecting cost, preventing it from pushing slices to unrealistic locations;
on the other hand, it can also improve the correction accuracy since it is derived from a strong physical property of the LV anatomy.

The discontinuity is measured by the dissimilarity between two regions of the same area sampled from adjacent SA slices.
Although the LV appears different in adjacent slices due to different imaging locations, the change is gradual.
Therefore, by minimizing the dissimilarity via translation of slices, the original relative positions of the SA slices can be recovered.
Again the ROI is predefined for each SA slice.
However, since the ROI is defined as the bounding box which well contains the LV, its size is varied for most SA slices as the LV size varies in most SA slices;
furthermore, during the translation, SA slices are deviated from each other, and so are their ROIs.
To solve this problem, we project the two ROIs of adjacent SA slices onto a middle plane and find the smallest rectangle which can completely contain the two projections on that plane, and then project this rectangle back onto the SA slices under consideration to determine the regions used for the discontinuity measure (Fig.~\ref{fig:ROI}(b)).
Noting that the projections are all along the normal to the SA planes, these procedures are actually finding the smallest cuboid whose top and bottom faces are in the SA slices (planes) under consideration and can completely contain the corresponding ROI in each slice.
Denoting the two regions sampled from the two adjacent SA slices by $R^{k}_{\mathrm{SA}}$ and $R^{k+1}_{\mathrm{SA}}$, the contiguous cost is defined as:
\begin{equation}\label{eq:contiguousCost}
    E_{\mathrm{cnt}}(I^{k}_{\mathrm{SA}}\parallel I^{k+1}_{\mathrm{SA}})=\mathcal{C}_{R}(R^{k}_{\mathrm{SA}}, R^{k+1}_{\mathrm{SA}}),\, k\in[1,m-1],
\end{equation}
where $I^{k}_{\mathrm{SA}}\parallel I^{k+1}_{\mathrm{SA}}$ denotes the adjacency of two consecutive SA slices, and $\mathcal{C}_{R}$ is a function measuring dissimilarity between two regions.
$E_{\mathrm{cnt}}$ is also a function of IPP, as the regions sampled on adjacent SA slices are determined by their IPPs (plus the constants IOP and PS).
Furthermore, because the translation of slices takes place in the patient coordinate system (instead of any of those coinciding with the SA plane), $E_{\mathrm{cnt}}$ is fully dependent on slices' translations in all $x$-, $y$-, and $z$- directions.

\subsubsection{Total Cost}
The total cost is a summation of all the intersecting and contiguous costs.
The intersecting cost is computed for all intersecting slices, including intersections between SA and LA slices, and between different LA slices;
each intersection is counted only once.
The contiguous cost is computed for all pairs of adjacent SA slices;
each adjacency is counted once too.
Now we can write the total cost as:
\begin{align}\label{eq:finalCost}
    &E_{\mathrm{align}}(\textit{IPP}_{\mathrm{all}}) =
    \sum^{m}_{k=1}\sum^{n}_{j=1} E_{\mathrm{int}}(I^{k}_{\mathrm{SA}}\bigcap I^{j}_{\mathrm{LA}})\ + \\
    &\sum^{n-1}_{j=1}\sum^{n}_{j'=j+1} E_{\mathrm{int}}(I^{j}_{\mathrm{LA}}\bigcap I^{j'}_{\mathrm{LA}})
    +\gamma\sum^{m-1}_{k=1} E_{\mathrm{cnt}}(I^{k}_{\mathrm{SA}}\parallel I^{k+1}_{\mathrm{SA}}), \nonumber
\end{align}
where $\textit{IPP}_{\mathrm{all}}=\{\textit{IPP}^{k}_{\mathrm{SA}},\ \textit{IPP}^{j}_{\mathrm{LA}} \mid k\in[1,m], j\in[1,n]\}$ represents $x$, $y$, $z$ coordinates of all slices under the patient coordinate system, and $\gamma$ a weighting factor.
$E_{\mathrm{align}}$ is minimized to find the optimal $\textit{IPP}^{\mathrm{opt}}_{\mathrm{all}}$ set, which is the corrected positions of all the slices.

For both $\mathcal{C}_{\boldsymbol{s}}$  and $\mathcal{C}_{R}$ we use the mean of the squared differences.
Noting that the intensities of the same structures may vary in different slices, the sampled $\boldsymbol{s}^{k}_{\mathrm{SA}}$, $\boldsymbol{s}^{j}_{\mathrm{LA}}$ and $R^{k}_{\mathrm{SA}}$ are pre-normalized to zero mean and unit standard deviation before being fed to $\mathcal{C}_{\boldsymbol{s}}$  and $\mathcal{C}_{R}$.
The ROI can be automatically estimated by a variant of the Hough transform tailored to the detection of annular shapes \cite{relatedworks01_Ciofolo2008} and our method is insensitive to small variations of the ROI.
The weighting factor $\gamma$ is set to 0.01 empirically.
The second row of Fig.~\ref{fig:misalignment_n_correction_n_intersection} shows two examples of the misalignment correction results.
We can see that the artifacts are largely reduced or even eliminated.

\subsection{3D Segmentation of the Myocardium}
\label{sec:method:subsec:myoSeg}
For high-quality automatic segmentation of the myocardium, we proposed a sophisticated 3D method.
Given myocardial contours in cine images as \emph{a priori} knowledge, either (semi-)automatically or manually obtained, the method initially propagates the \emph{a priori} segmentation from cine to LGE images via 2D translational registration.
Two \emph{simplex} meshes \cite{simplexmeshIJCV1999} representing respectively endocardial and epicardial surfaces are then constructed with the propagated segmentation contours.
After construction, the two meshes are nonrigidly deformed towards the myocardial edge points detected in both SA and LA LGE images in a unified 3D coordinate system.
Based on intensity characteristics of the LV in LGE images, we propose a novel parametric model of the LV for consistent myocardial edge point detection regardless of pathological status of the myocardium (infarcted or healthy) and of the type of the LGE images (SA or LA).
The final meshes after the nonrigid deformation are themselves a 3D segmentation of the myocardium.

Our method is distinct from the existing ones \cite{relatedworks01_Ciofolo2008,relatedworks02_Dikici2004,myownMICCAI2011} in two aspects.
(i) We integrate LA images (standard 4-chamber (4C) and 2-chamber (2C) views) with SA images for the 3D segmentation.
Besides providing complementary information about the LV between the largely spaced SA images, the LA images are also used for the correction of misalignment artifacts among slices.
(ii) We propose a novel parametric model of the LV for LGE images based on 1D intensity profiles.
This model is flexible to be applied to both SA and LA images, and self-adaptive to detect edge points of both infarcted and healthy myocardium.

We have evaluated the proposed method with the real patient data which is also used in this work.
The mean distance error between the automatic contours by our method and manual contours by the expert is $0.97\pm0.69$~mm, and the mean Dice coefficient (DC) \cite{Dice1945} for the myocardium is $87.60\pm8.63\%$.
The results confirm the observation that the method can generate accurate and reliable results for myocardial segmentation of LGE images.
We have also tested the robustness of the method with respect to varied \textit{a priori} segmentation.
In one of the experiments, we use manual and automatic cine segmentations as \textit{a priori} knowledge for a set of LGE data, and find that the segmentation results of the LGE images are similarity good: $1.02\pm0.18$ mm versus $1.12\pm0.29$~mm and $87.84\pm2.73\%$ versus $86.18\pm3.76\%$.
These results demonstrate that the method can compensate variations in the given \textit{a priori} knowledge and consistently produce accurate segmentations.

Due to the page limit, we will present the detailed description and comprehensive evaluation of this method in another paper.

\subsection{Correction of Intensity Inconsistencies}

Signal intensities may vary from one LGE image to another even in the same dataset due to different time delays after injection of the contrast agent.
A usual pattern is that the intensities tend to become brighter from the mitral valve to the apex (see the first row of Fig.~\ref{fig:intensityCorrectionComparison}).
Therefore, for any infarct classification method that combines multiple slices, normalization of the intensity across slices is necessary.
We propose a novel method which only considers local regional intensities of the LV and thus can effectively eliminate or reduce the intensity inconsistencies in the LV region.
The rationale underlying our normalization is that the intensity of the BP should be roughly the same across slices.
The BP region in each image can be simply obtained by masking the region enclosed by the endocardial contour.
However, dark papillary muscles are also considered as BP pixels in this way.
In order to normalize the BP and hence LV intensities more accurately, we have implemented an iterative algorithm based on the intensity distribution of the LV.

\begin{figure*}[!t]
  \centering
  \includegraphics[width=\textwidth]{intensityCorrectionComparison}\\
  \caption{Illustration of the intensity normalization.
  First row: a stack of original SA images in one dataset; images become brighter from the mitral valve to the apex.
  Second row: the same image stack after the correction of intensity inconsistencies; the image intensities, especially in the LV regions, are more consistent.}\label{fig:intensityCorrectionComparison}
\end{figure*}

\subsubsection{Rician Distribution of the LV in LGE CMR Images}
In the presence of noise, the intensity of MR images is shown to be governed by a Rician distribution \cite{ricianMRI1995}.
In the case of the LV in LGE images, the Rayleigh distribution models normal myocardium (dark regions) while the Gaussian distribution models infarcts plus the BP (bright regions), and the entire LV is modeled with the mixture model:
\begin{equation}\label{eq:RicianDistribution}
\begin{split}
    &Rcn(x) = Pr_\mathrm{Ray}(x) + Pr_\mathrm{Gau}(x)\\
    &=\underbrace{\alpha_{\mathrm{R}}\frac{x}{\sigma_{\mathrm{R}}^{2}}\exp(-\frac{x^{2}}{2\sigma_{\mathrm{R}}^{2}})}_{\text{\footnotesize{Rayleigh distribution}}} +\underbrace{\alpha_{\mathrm{G}}\frac{1}{\sqrt{2\pi}\sigma_{\mathrm{G}}}\exp[-\frac{1}{2}(\frac{x-\mu}{\sigma_{\mathrm{G}}})^2]}_{\text{\footnotesize{Gaussian distribution}}},
\end{split}
\end{equation}
where $Pr_\mathrm{Ray}(x)$ and $Pr_\mathrm{Gau}(x)$ denote the Rayleigh and Gaussian distributions respectively, and $\alpha_\mathrm{R}$, $\delta_\mathrm{R}$, $\alpha_\mathrm{G}$, $\delta_\mathrm{G}$, $\mu$ their parameters.
The combination of the BP and infarcts as one class is not only reasonable because there is no significant difference between their intensity distributions \cite{infarctSeg2005}, but also helpful to make the assumed bimodal distribution more prominent \cite{infarctSeg2010b}.
During the acquisition of LGE CMR images, technologists try to null signals from the normal myocardium as much as possible in order to give prominence to the hyper-enhanced regions.
Consequently, the intensities of the normal myocardium is often shifted to the dark side with an unknown offset.
It has been suggested that introducing an extra offset parameter $a$ to the Rayleigh distribution would lead to more accurate estimation of (\ref{eq:RicianDistribution}) \cite{infarctSeg2010a}.
We follow this approach and thus $Pr_\mathrm{Ray}(x)$ becomes:
\begin{equation}
    Pr_\mathrm{Ray}(x) = \alpha_{\mathrm{R}}\frac{(x+a)}{\sigma_{\mathrm{R}}^{2}}\exp\left[-\frac{(x+a)^{2}}{2\sigma_{\mathrm{R}}^{2}}\right].
\end{equation}

We extract all the LV voxels from the regions enclosed by the epicardial contours in the stack of SA slices and obtain the intensity histogram.
We then obtain a \emph{relative probability} distribution by normalizing the histogram with the largest frequency count and fit the relative probability distribution to (\ref{eq:RicianDistribution}) to solve for $Rcn(x)$.
Figure~\ref{fig:RicianFit} shows an example of the fitted Rician distribution overlaid on the relative probability distribution.
The intersection of $Pr_\mathrm{Ray}(x)$ and $Pr_\mathrm{Gau}(x)$ is found and denoted as $i_\mathrm{thrh}$.
Thereafter BP voxels can be extracted from the regions enclosed by the endocardial contours but excluding those with intensities smaller than $i_\mathrm{thrh}$ (considered papillary muscles).
We take the middle slice in the SA stack as the reference and calculate the mean intensity of its BP pixels $i_\mathrm{m\_ref}$.
Then the rest of the slices are normalized by multiplying the factor $i_\mathrm{m\_ref} / i_\mathrm{m\_k}$, where $i_\mathrm{m\_k}$ is the mean intensity of the BP pixels in the $k^\mathrm{th}$ slice.
Finally, the whole stack of slices are normalized \emph{together} to the range $[0,1]$.

\begin{figure}
  \centering
  \includegraphics[width=.95\columnwidth]{RicianFit}\\
  \caption{The fitted Rician distribution overlaid on the relative probability distribution.
  Also overlaid are vertical position lines of $\sigma_\mathrm{R}-a$, $i_\mathrm{thrh}$ and $\mu$.}\label{fig:RicianFit}
\end{figure}

\subsubsection{Iterative Normalization}

Since the above described procedure may have changed the shape of the underlying intensity histogram of the stack of LGE slices, it is likely that if we repeat we would still find $i_\mathrm{m\_ref} / i_\mathrm{m\_k}$ unequal to 1 (as $i_\mathrm{thrh}$ is changed along with the histogram shape).
Therefore, we iterate the normalizing procedure until the histogram stabilizes, which is determined by examining how close the ratio $i_\mathrm{m\_ref} / i_\mathrm{m\_k}$ is to 1 for every slice.
Figure~\ref{fig:intensityCorrectionComparison} shows slices of a dataset before and after the intensity normalization as a comparison.
The final $Pr_\mathrm{Ray}(x)$ and $Pr_\mathrm{Gau}(x)$ obtained when the iteration stops will be used further for infarct classification by the graph-cut algorithm.

\subsection{Classification of Infarcts by Graph-Cut}
\label{sec:subsec:graphCutClassify}

\subsubsection{3D Graph-Cut}
We employ the ND graph-cut algorithm presented in \cite{graphCutJolly2001} for our 3D classification of infarcts.
Let $p$, $q$ be any two different pixels in the image $I$ and $l_p$ the label assigned to $p$ in a two-label system $l_p\in\{0,1\}$, then a labeling set $L$ of every pixel in $I$: $L=\{l_p\ |\ p\in I\}$ represents a segmentation of $I$.
Based on the concept of maximum \textit{a posteriori} estimate in a locally dependent Markov random field \cite{MRFgraphcut1998}, the energy functional used by the graph-cut algorithm can be written as:
\begin{equation}
  E_\mathrm{seg}(L) = \lambda \sum_{p\in I} -\ln Pr(i_p\ |\ l_p) + \sum_{p\in I} \sum_{q\in\mathcal{N}_p} V_{(p,q)}(l_p,l_q),
\end{equation}
where $Pr(i_p\ |\ l_p)$ is the \emph{relative} likelihood, $V_{(p,q)}$ a \emph{clique potential}, $\mathcal{N}_p$ the set of all neighboring pixels of $p$, and $\lambda$ a positive weighting factor.
The term $-\ln Pr(i_p\ |\ l_p)$ denotes the cost of assigning label $l_p$ to the pixel $p$ of intensity $i_p$, and $V_{(p,q)}(l_p,l_q)$ denotes the cost of assigning labels $l_p$ and $l_q$, respectively, to the two neighboring pixels $p$ and $q$.

The input to the graph-cut algorithm is just the segmented myocardium, that is, regions other than the myocardium are masked out to eliminate irrelevant information.
Note that the input is a stack of 2D images, i.e., a 3D image.
Since there are only two allowable values for $l_p$, i.e., 1 for infarct and 0 for normal myocardium, $Pr(i_p\ |\ l_p)$ is defined as:
\begin{align}
    Pr(i_p\ |\ 1) &= \begin{cases}
        Pr_\mathrm{Gau}(i_p),\quad &\text{if}\quad i_p\leq \mu \\
        1,\quad                    &\text{if}\quad i_p>\mu
    \end{cases}\\
    Pr(i_p\ |\ 0) &= \begin{cases}
        Pr_\mathrm{Ray}(i_p),\quad &\text{if}\quad i_p\geq(\sigma_\mathrm{R}-a) \\
        1,\quad                    &\text{if}\quad i_p<(\sigma_\mathrm{R}-a)
    \end{cases}
\end{align}
The two special cases -- $Pr(i_p\ |\ 1)=1$ for $i_p$ larger than $\mu$ and $Pr(i_p\ |\ 0)=1$ for $i_p$ smaller than $(\sigma_\mathrm{R}-a)$ -- conforms to the clinical experience that a pixel brighter than the mode of the Gaussian distribution must be given the infarct label 1 and one darker than the mode of the Rayleigh distribution must be given the myocardium label 0.
In these cases, there should be no cost on the corresponding labelings, which is embodied by the fact that the logarithm of 1 is 0.

We also use the interaction potential $V_{(p,q)}$ proposed in \cite{graphCutJolly2001}.
However, instead of setting $\sigma$ -- the constant that controls the effective range of the potential well -- to be the estimated image noise, we set $\sigma=\mu-(\sigma_\mathrm{R}-a)$, which is the distance between the modes of $Pr_\mathrm{Ray}(x)$ and $Pr_\mathrm{Gau}(x)$.
$E_\mathrm{seg}(L)$ is minimized with the efficient $\alpha$-expansion algorithm \cite{alphaExpansion1999}.
Because only two labels are involved in our specific problem, the exact global minimum of $E_\mathrm{seg}(L)$ can be reached \cite{firstGraphCut}.
As for the definition of $\mathcal{N}_p$, we use the 6-connected neighborhood.
$\lambda$ is set to 1 at all times.

\subsubsection{Post-Processing}

False positive of infarcts is usually caused by inclusion of thin layers of epicardial fats or endocardial BP inside the myocardial contours.
This kind of contour tracing flaws are inevitable for both manual and \mbox{(semi-)} automatic delineations.
Sometimes patches of artifacts within the myocardium are another source of false positive.
Meanwhile, false negative may happen at regions with intermediate intensities, which are caused by the partial volume effect, i.e., the mixing of infarcted and normal myocardium.
Therefore, post-processing is necessary to remove these misclassifications.
Since there are already extensive works on the post-processing techniques, we simply follow those proposed in \cite{infarctSeg2010b} for the implementation, which have been proven effective and also work readily in our experiments.
To include MVO, sub-endocardial dark blobs surrounded by hyper-enhanced infarcts are located and classified as infarcts as well.

\subsection{Implementation}

Our framework is implemented under the \textsc{Matlab} R2011a environment, on a PC with Intel Core2 Duo T9400 processer, 3 GB RAM and running a 32-bit Microsoft Windows 7 OS.
Optimization of $E_{\mathrm{align}}$ and fitting of $Rcn(x)$ are achieved with the Optimization Toolbox and Curve Fitting Toolbox of \textsc{Matlab} R2011a, respectively.
The average running times for each step in the framework are: 26 s for the misalignment correction, 2.1 min for the 3D myocardium segmentation and 1.5 s for the intensity normalization plus 3D infarct classification.
The running times can be significantly reduced if the proposed framework is implemented in a more efficient programming language with multi-threading enabled.
The framework is found insensitive to the exact values of the parameters $\gamma$ and $\lambda$ by our preliminary experiments.

\section{Experimental Results and Discussion}
\label{sec:results}

\subsection{Data Specification and Experimental Settings}
\label{sec:subsec:dataSpecs}

20 sets of real patient data (denoted by DAT$_{1-20}$) with varying amounts of infarcts are used for validation of the proposed framework.
The data were acquired with ECG gating by a 1.5T Siemens Symphony MRI scanner from 20 patients (18 males and 2 females, 38-81 years old) three months after myocardial infarction, following a bolus injection of gadolinium-based contrast agent.
The sequence parameters used for the LGE image acquisition are as follows: flip angle 25$^\circ$, repetition time 650-1000 ms, echo time 4.18 ms, bandwidth 130 Hz/pixel, image width and height 176-256 pixels, in-plane isotropic pixel size 1.17-1.56 mm, slice thickness 7 mm and gap 3 mm.
Depending on the individual heart size, there could be 6-8 usable SA slices within the LV in each dataset;
standard 4C and 2C LA views were acquired along with the SA views.
In total 149 SA slices were included in the study.
The data were analyzed by medical experts (including segmentation of the myocardium and infarcts) manually and the results were used as the reference standard in our experiments.
The percentage of the infarcts with respect to the entire myocardium is commonly used by cardiologists to reflect the extent and severity of infarction.
According to the reference standard, volumetric infarct percentages for the 20 patients range from 0 to $32.31\%$ (two cases of non-infarct) with the mean and standard deviation of 18.60$\pm$10.43\%, with five patients having MVO.

For myocardium segmentation, the general quantitative evaluation results are already provided in Section \ref{sec:method:subsec:myoSeg}.
In this paper we focus on infarct classification and the final quantification results.
Apart from the infarct percentage (denoted by I/M\% afterwards), the DC is employed to measure the overlapping rates of infarcts classified automatically and manually.
I/M\% calculated from the manual and automatic results are compared by the Bland-Altman (BA) analysis \cite{Bland-Altman1986}.

\subsection{Quantification Results}

Figure \ref{fig:finalResults} shows the automatic and manual results on 3 LGE images with different extents of infarction and different intensities.
Qualitatively we find that our framework consistently produces myocardial and infarct segmentations close to those by the experts, though several exceptions of failed cases are found in few very apical slices (e.g., the slice shown in Fig.~\ref{fig:failure}).

\begin{figure}
    \centering
    \includegraphics[width=.95\columnwidth]{001B_SA04}\\
    \includegraphics[width=.95\columnwidth]{003B_SA06}\\
    \includegraphics[width=.95\columnwidth]{055B_SA06}\\
    \begin{minipage}{.95\columnwidth}
        \centering\footnotesize
        \makebox[.32\linewidth][c]{(a)} \hfill \makebox[.32\linewidth][c]{(b)} \hfill \makebox[.32\linewidth][c]{(c)}
    \end{minipage}\\
    \caption{Exemplary results (each row shows a slice from a different subject): (a) the original image; (b) the proposed framework; (c) the reference standard.
    The green and yellow contours delineate the endocardium and epicardium respectively, and the red contours delineate the infarcts.}\label{fig:finalResults}
\end{figure}

\subsubsection{Volumetric Analysis}
Of the 20 sets of data, 2 sets present no infarct after three months' recovery since the heart attack and our
framework correctly reported zero infarct.
The average DC for the other 18 sets of data is $80.90\pm5.04\%$, with the minimum and maximum values of $67.23\%$ and $88.25\%$, respectively.
Figure \ref{fig:Bland-Altman_plots} top shows the BA plot of the volumetric I/M\% (automatic minus manual results).
The BA analysis indicates that the automatic results report slightly higher percentages of infarcts than the manual results with a mean difference of $0.14\pm3.07\%$.
Furthermore, a Wilcoxon rank sum test indicates no significant difference between respective I/M\% ($P=0.95$).

\begin{figure}
    \rotatebox{90}{\hspace{6mm} \footnotesize Automatic I/M\%$-$Manual I/M\%}
    \begin{minipage}[b]{.45\textwidth}
      \centering
      \includegraphics[width=\linewidth]{Bland-Altman_plot}\\
      \footnotesize (Automatic I/M\%$+$Manual I/M\%)$/$2\vspace{3mm}
      \label{fig_first_case}
    \end{minipage}
    \\
    \rotatebox{90}{\hspace{2.5mm} \footnotesize Automatic I/M\%$-$Manual I/M\%}
    \begin{minipage}[b]{.45\textwidth}
      \centering
      \includegraphics[width=\linewidth]{Bland-Altman_allSeg}\\
      \footnotesize (Automatic I/M\%$+$Manual I/M\%)$/$2
      \label{fig_first_case}
    \end{minipage}

    \caption{Bland-Altman plots of volumetric (top) and AHA segment-wise (bottom) I/M\%.}\label{fig:Bland-Altman_plots}
\end{figure}

\subsubsection{Segment-wise Analysis}
The American Heart Association (AHA) has recommended to divide the LV myocardium into 17 segments for standardization \cite{AHA17segments};
since the 17$^\mathrm{th}$ segment, which is extracted from the apex in 4C and 2C LA views, is rarely used for infarct quantification by cardiologists in LGE CMR studies, only the 1$^\mathrm{st}$-16$^\mathrm{th}$ segments (denoted by Seg$_{1-16}$) extracted from SA views are included in our study.
As the segments are standardized with assignments to coronary arterial territories, examining segment-wise infarct presence can help identify the correlated vessel that is potentially occluded.
For easy visual interpretation, I/M\%'s of the 16 segments are plotted together in color-coded Bull's eye plots (Fig. \ref{fig:BullseyePlot}).

\begin{figure}
  \centering
  \rotatebox{90}{\hspace{5mm} \footnotesize{Infarct percentage}}
  \includegraphics[width=.9\columnwidth]{BullseyePlot}\\
  \leftline{\footnotesize
    \hspace{3.1cm}(a) \hspace{3.2cm}(b)
  }
  \caption{Bull's eye plots of I/M\%'s in standardized 16 segments for a set of LGE CMR data: (a) the automatic results; (b) the manual results. The extent of infarction is coded by the hot map. }\label{fig:BullseyePlot}
\end{figure}

Compared to the LV volume-wise analysis, segment-wise analysis is more sensitive and expected to reveal more insights due to the employment of location-specific smaller volumes.
The BA analysis with all the $20 \times 16 = 320$ segments (Fig.~\ref{fig:Bland-Altman_plots} bottom) indicates that overall the automatic results report slightly higher I/M\% than the manual results, with a difference of $0.66\pm6.68\%$ ($P = 0.99$).
While the data points look compact around the zero level line at the smaller end of average I/M\%, the automatic results tend to report considerably more infarcts than the manual results when the average I/M\% is larger than 75\% (see the red dashed line in Fig.~\ref{fig:Bland-Altman_plots} bottom).

Of the 320 segments, 142 segments are deemed infarcted and 173 segments deemed free of infarct by both results.
Two false positives (i.e., the manual results reported zero while the automatic results reported positive I/M\%) are found with very small values: $3.06\%$ and $1.56\%$, respectively.
Meanwhile, three false negatives (i.e., the manual results reported positive while the automatic results reported zero I/M\%) are found and they are: $-3.03\%$, $-28.49\%$ and $-6.82\%$.
We have also compared absolute differences between I/M\% calculated from the automatic and manual results at different slice levels: basal (Seg$_{1-6}$), mid (Seg$_{7-12}$) and apical (Seg$_{13-16}$).
The multiple comparison following a Kruskal-Wallis test indicates that while there is no significant difference between basal ($2.14\pm5.21\%$) and mid ($3.12\pm5.47\%$) levels, the absolute differences at the apical level ($5.12\pm6.94\%$) are significantly larger than both of the basal and mid levels ($P\ll0.05$).

\subsection{3D Versus 2D Quantifications}

For 2D segmentation of the myocardium, we employed the method proposed in \cite{myownMICCAI2011}.
On the 20 sets of data, the mean distance error is $1.09\pm0.80$ mm, and the DC for the myocardium is $85.84\pm9.86\%$.
For 2D classification of infarcts, we applied the method presented in Section \ref{sec:subsec:graphCutClassify} to one slice at a time without the volumetric correction of intensity inconsistencies.
Also, the post processing took place in 2D.
As the 2D quantification works slice-wisely, naturally we compare the quantification results at the slice level.
The comparison is shown in Table \ref{tab:2d_vs_3d}.
The quantification results of our 3D framework excel those of the 2D methods in three of the four compared aspects, and are equally good in the last.
\begin{table}[h]
  \centering
  \caption{3D versus 2D quantifications. Note: the DCs were calculated with slices deemed infarcted by all manual, 2D and 3D results~($N=102$).}\label{tab:2d_vs_3d}
  \begin{tabular}{c|c|c|c}
    \hline

     & 2D & 3D & p-value \\
    \hline
    DC (\%) & $73.09\pm16.17$ & $79.11\pm15.50$ & $\ll0.05$ \\
    Absolute diff. in I/M\% & $5.64\pm6.54$ & $3.56\pm5.04$ & 0.01 \\
    No. false positives & 13 & 3 & -- \\
    No. false negatives & 1 & 1 & -- \\
    \hline
  \end{tabular}
\end{table}

\subsection{Discussion}

\subsubsection{Accuracy and Applicability of the Framework}

Experimental results on the basis of LV volumes indicate that our framework is capable of producing a reliable evaluation of the overall extent and severity of infarction for a patient.
The volumetric BA analysis reports a nearly zero discrepancy ($0.14\%$) between the percentages of infarcts derived from the automatic and manual results with a very narrow limit of agreement ($SD = 3.07\%$).
The volumetric DCs also suggest that the automatic and manual results largely agree with each other (mean DC $= 80.90\pm5.04\%$).
However, we should be cautious when applying the framework to evaluate the extent of infarction in each AHA segment and identify potential supplying arteries of stenosis.
Special attention should be paid to quantification results of segments in the apical level, i.e., Seg$_{13-16}$, since the Kruskal-Wallis test with the slice levels as the factor indicates that discrepancies between the automatic and manual results are significantly larger at the apical level than at the basal and mid levels ($5.12\pm6.94\%$ versus $2.14\pm5.21\%$ and $3.12\pm5.47\%$).
Further, the only significant false negative ($-28.49\%$) out of the five cases of false errors is Seg$_{15}$ of DAT$_5$.
Therefore, the framework can still be applied to segment-wise quantification of infarcts if careful attention is paid to apical segments and proper correction schemes (where necessary) are provided.

\subsubsection{Implications}

As mentioned, we find that discrepancies between the automatic and manual results are significantly larger at the apical level than at the basal and mid levels.
This may be due to two reasons:
(i) Our framework tends to report considerably higher infarct percentages than the human observers when there are a large amount of infarcts.
Naturally there are significantly more infarcts at the apical than at the other two levels in real patients \cite{infarctSeg2005}, so our framework tends to report more infarcts for apical segments.
(ii) Automatic segmentation of the myocardium in apical slices is difficult, and myocardial contours of very poor quality cause failures anyway despite the good performance of the subsequent infarct classification method.
Such cases are illustrated in Fig.~\ref{fig:failure} with the most apical slice in DAT$_5$.
This failure occurs because our method utilizes myocardial contours in corresponding cine CMR images as the \textit{a priori} knowledge to guide the segmentation in LGE images.
When there is a large change in myocardium shape between corresponding cine and LGE images, the method may not be able to compensate and the segmentation fails.
Due to the use of ECG gating, such severe changes only happen at the extremely apical slices, where the LV is small and the apex's motion is large.
In fact the segmentation of the myocardium at apical slices is always so difficult that no method can consistently produce reliable segmentations.
Therefore, a user-correction scheme should be provided to handle such cases.

\begin{figure}
  \centering
  \includegraphics[width=.95\columnwidth]{006B_SA08}\\
  \begin{minipage}{.95\columnwidth}
        \centering\footnotesize
        \makebox[.32\linewidth][c]{(a)} \hfill \makebox[.32\linewidth][c]{(b)} \hfill \makebox[.32\linewidth][c]{(c)}
    \end{minipage}\\
  \caption{The most apical slice in DAT$_5$: (a) the original LGE image; (b) the automatic result; (c) the manual result.
  The incorrect classification of infarcts is caused by the failed automatic segmentation of the myocardium.}\label{fig:failure}
\end{figure}

\subsubsection{Advantages of 3D Quantification}
The experimental results show that our 3D quantification framework outperform the 2D methods in the control group in three of the four examined aspects and has the same good performance in the last.
There are two reasons for this.
First, the myocardial contours produced by our 3D segmentation method is more accurate than those by the 2D method.
Second, due to the utilization of 3D spatial and intensity information in a stack of SA slices, the 3D classification of infarcts is also better than a 2D application of the same method to individual slices.
The latter has been confirmed by a subsequent comparative evaluation of merely the 3D and 2D \textit{classification} methods, in which the same manual myocardial contours drawn by the experts were input to both the 3D and 2D classification methods.
The results show that the 3D classification still performs better in the first three examined aspects in Table \ref{tab:2d_vs_3d} and equally well in the last.

\subsubsection{Limitations and Future Works}
One potential limitation of our framework is that the contiguous cost in the misalignment correction step does not work for sparse SA slices, e.g., when there is only one slice for each of the basal, mid and apical levels of the LV.
However, our misalignment correction method would still work for such cases by setting the weight for the contiguous cost $\gamma=0$, in which case the method is purely intersection based and thus quite similar to \cite{alignManyLA2004}.
Fortunately, most LGE CMR data nowadays have SA slices with slice spacings smaller than 1 cm (unless otherwise intentionally scanned in a sparse way) and our method utilizes this advantageous feature when available.
For future work, it would be interesting to study the impacts of different slice spacings and pixel sizes on the entire framework.

The evaluation in this work was conducted with limited (i.e., 20 sets) LGE CMR data.
If more datasets are involved, the framework could be tested by more potential variations of the human heart anatomy.
Moreover, reliable Kruskal-Wallis tests with the AHA segment order as the factor would be possible with a considerably larger database, other than the slice level used in this work.
This would provide insights specific to every AHA segment.
Another drawback of this work is that we did not explore the potential variations between different human observers.
As a study reported volumetric DCs of $80\pm8\%$ between the infarcts segmented by two observers given the same myocardial contours \cite{infarctSeg2010b}, we expect that, if the observers draw and use their respective delineations of the myocardium, the DC would drop further.

It has recently been reported that the infarcts may not always be homogeneous and that tissue heterogeneity (core and peri-infarct zones) has great diagnostic and predictive potentials~\cite{yan2006characterization,schmidt2007infarct}.
In future work, we plan to extend our framework to also support the quantification of the infarct core and peri-infarct zones within the identified infarcts.
Finally, though the 17$^{th}$ segment in the AHA nomenclature, which stands for the apex in the LA views, is rarely used for the quantification of infarcts, it can provide valuable information on the presence and transmurality of the infarcts right at the apex.
Hence we are currently extending our 3D segmentation method to the apex in 4C and 2C LA views to support the analysis of infarcts there.

\section{Conclusion}
\label{sec:conclusion}

This paper presents a complete 3D framework for automatic quantification of LGE CMR images.
This framework achieves 3D segmentation of the myocardium as well as 3D classification of infarcts within the segmented myocardium, with robust and effective pre-processing measures overcoming misalignment artifacts and intensity inconsistencies across slices.
Qualitative and quantitative evaluations using real patient data demonstrated that the proposed framework is able to produce accurate segmentation and classification results and has the potential to be developed further as a clinical tool to generate objective quantification of LGE CMR images.

\appendices

\section*{Acknowledgment}

The authors would like to thank the Academic Research Fund, National University of Singapore, Ministry of Education, Singapore for funding the CMR studies. We are also grateful to the radiographers and staff at the Department of Diagnostic Imaging, National University Hospital, Singapore, for helping with the CMR scans.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtranN}