\section{Introduction}

Most of the events that we experience in our life consist of visual and acoustic stimuli. Recordings of such events may lack the acoustic component, for example due to limitations of the recording equipment or technical issues in the transmission of the information. Since acoustic and visual modalities are often correlated, methods to reconstruct audio signals using videos have been proposed \cite{davis2014visual, owens2016visually, zhou2018visual}.

In this paper, we focus on one particular case of the aforementioned problem: \textit{speech reconstruction (or synthesis) from a silent video}. Solving this task might be useful to automatically generate speech for surveillance videos and for extremely challenging speech enhancement applications, e.g. hearing assistive devices, where noise completely dominates the target speech, making the acoustic signal worth less than its video counterpart.

A possible way to tackle the problem is to decompose it into two steps: first, a \textit{visual speech recognition} (VSR) system \cite{assael2016lipnet, Stafylakis2017, chung2017lip} predicts the spoken sentences from the video; then, a \textit{text-to-speech} (TTS) model \cite{sotelo2017char2wav, ping2018deep, shen2018natural} synthesises speech based on the output of the VSR system. However, at least two drawbacks can be identified when such an approach is used. In order to generate speech from text, each word should be spoken in its entirety to be processed by the VSR and the TTS systems, imposing great limitations for real-time applications. Furthermore, when the TTS method is applied, useful information that should be captured by the system, such as emotion and prosody, gets lost, making the synthesised speech unnatural. For these reasons, approaches that estimate speech from a video, without using text as an intermediate step, have been proposed.

Le Cornu and Miller \cite{cornu2015reconstructing, le2017generating} developed a video-to-speech method with a focus on speech intelligibility rather than quality. This is achieved by estimating spectral envelope (SP) audio features from visual features and then reconstructing the time-domain signal with the STRAIGHT vocoder \cite{kawahara1999restructuring}. Since the vocoder also requires other audio features, i.e. the fundamental frequency (F0) and the aperiodic parameter (AP), these are artificially created independently of the visual features.

Ephrat and Peleg \cite{ephrat2017vid2speech} treated speech reconstruction as a regression problem using a neural network which takes as input raw visual data and predicts a line spectrum pairs (LSP) representation of linear predictive coding (LPC) coefficients computed from the audio signal. The waveform is reconstructed from the estimated audio features using Gaussian white noise as excitation, producing unnatural speech. This issue is tackled in a subsequent work \cite{ephrat2017improved}, where a neural network estimates the mel-scale spectrogram of the audio from video frames and optical flow information derived from the visual input. The time-domain speech signal is reconstructed using either example-based synthesis, in which estimated audio features are replaced with their closest match in the training set, or speech synthesis from predicted linear-scale spectrograms.

Akbari et. al. \cite{akbari2018lip2audspec} tried to reconstruct natural sounding speech using a neural network that takes as input the face region of the talker and estimates bottleneck features extracted from the auditory spectrogram by a pre-trained autoencoder. The time-domain signal is obtained with the algorithm in \cite{chi2005multiresolution}. This approach shows its effectiveness when compared to \cite{ephrat2017vid2speech}.

All the methods reported until now have a major limitation: they estimate either a magnitude spectrogram, SPs or LSPs, which do not contain all the information of a speech signal. Vougioukas et al. \cite{Vougioukas2019} addressed this issue and proposed an end-to-end model that can directly synthesise audio waveforms from videos using a generative adversarial network (GAN). However, their direct estimation of a time-domain signal causes artefacts in the reconstructed speech.

In this work, we propose an approach, \textit{vid2voc}, to estimate WORLD vocoder \cite{morise2016world} features from the silent video of a speaker\footnote{Although this paper aims at synthesising speech from frontal-view silent videos, it is worth mentioning that some methods using multi-view video feeds have also been developed \cite{kumar2018harnessing, kumar2018mylipper, kumar2019lipper, Uttam2019}.}. We trained the systems using either the whole face or the mouth region only, since previous work \cite{ephrat2017vid2speech} shows a benefit in using the entire face. Our method differs from the work in \cite{cornu2015reconstructing, le2017generating}, because we predict all the vocoder features (not only SP) directly from raw video frames. The estimation of F0 and AP, alongside with SP, allows to have a framework with a focus on speech intelligibility (as in \cite{cornu2015reconstructing, le2017generating}) and speech quality, able to outperform even the recently proposed GAN-based approach in \cite{Vougioukas2019} in several conditions.
In addition, we train a system that can simultaneously perform speech reconstruction (our main goal) and VSR, in a multi-task learning fashion. This can be useful in all the applications that require video captioning without adding considerable extra complexity to the system. Although Kumar et al. \cite{kumar2019lipper} incorporate a text-prediction model in their multi-view speech reconstruction pipeline, this model is trained separately from the main system and it is quite simple: it classifies encoded audio features estimated with a pre-trained network into 10 text classes. This makes the method dependent on the number of different sentences of the specific database used for training and not suitable for real-time applications. Instead, we make use of the more flexible connectionist temporal classification (CTC) \cite{graves2006connectionist} sequence modelling which has already shown its success in VSR \cite{assael2016lipnet}.

Additional material, including samples of reconstructed speech that the reader is encouraged to listen to for a better understanding of the effectiveness of our approach, can be found in \url{https://danmic.github.io/vid2voc/}.

\section{Methodology and Experimental Setup}

\subsection{Audio-Visual Speech Corpus}

Experiments are conducted on the GRID corpus \cite{cooke2006audio}, which consists of audio and video recordings from 34 speakers (s1$-$34), 18 males and 16 females, each of them uttering 1000 six-word sentences with the following structure: $<$command$>$ $<$color$>$ $<$preposition$>$ $<$letter$>$ $<$digit$>$ $<$adverb$>$. Each video has a resolution of 720$\times$576 pixels, a duration of 3 s and a frame rate of 25 frames per second. The audio tracks have the same duration as the videos and a sample frequency of 50 kHz. In addition, text transcription for every utterance is provided.

As in \cite{Vougioukas2019}, we evaluate our systems in speaker dependent and speaker independent settings. Regarding the speaker dependent scenario, the data from 4 speakers (s1, s2, s4, s29) is pooled together, then 90\% of the data is used for training, 5\% for validation and 5\% for testing. Regarding the speaker independent scenario, the data from 15 speakers (s1, s3, s5$-$8, s10, s12, s14, s16, s17, s22, s26, s28, s32) is used for training, the data from 7 speakers (s9, s20, s23, s27, s29, s30, s34) for validation and the data from 10 speakers (s2, s4, s11, s13, s15, s18, s19, s25, s31, s33) for testing.

\subsection{Audio and Video Preprocessing}

The acoustic model used in this work is based on the WORLD vocoder \cite{morise2016world}, with a sample frequency of 50 kHz and a hop size of 250 samples\footnote{The window length is automatically determined by the WORLD algorithm.}. WORLD consists of three analysis algorithms to determine SP, F0 and AP features, and a synthesis algorithm which incorporates these three features. Here, we use SWIPE \cite{camacho2007swipe} and D4C \cite{morise2016d4c} to estimate F0 and AP, respectively. As done in \cite{blaauw2017neural}, a dimensionality reduction of the features is applied: SP is reduced to 60 log mel-frequency spectral coefficients (MFSCs) and
AP is reduced to 5 coefficients according to the D4C band-aperiodicity estimation. In addition, a voiced-unvoiced (VUV) state is obtained by thresholding the F0 obtained with SWIPE. All the acoustic features are min-max normalised using the statistics of the training set as in \cite{chandna2019vocoder}.

As in \cite{Vougioukas2019}, videos are preprocessed as follows: first, the faces are aligned to the canonical face\footnote{We use the face processor library in \url{https://github.com/DinoMan/face-processor}, which makes use of \cite{bulat2017far}.}; then, the video frames are normalised in the range $[-1, 1]$, resized to 128$\times$96 pixels and, for the models that use only the mouth region as input, cropped preserving the bottom half; finally, the videos are mirrored with a probability of 0.5 during training.

\subsection{Architecture and Training Procedure}

As shown in Figure \ref{fig:architecture}, our network maps video frames of a speaker to vocoder features and consists of a \textit{video encoder}, a \textit{recursive module} and five decoders: \textit{SP decoder}, \textit{AP decoder}, \textit{VUV decoder}, \textit{F0 decoder} and \textit{VSR decoder}. We also tried not to use the VSR decoder, to see whether it has any impact on the performance.

The video encoder is inspired by \cite{Vougioukas2019}: it takes as input one video frame concatenated with the three previous and the three next frames and applies five 3-D convolutions (conv3D). Each of the first four convolutional layers is followed by batch normalisation (BN) \cite{ioffe2015batch}, ReLU activation and dropout \cite{srivastava2014dropout}, while the last one is followed by Tanh activation.

To model the sequential nature of video data, a recursive module is used: it consists of a single-layer gated recurrent unit (GRU) \cite{cho2014learning}, BN, ReLU activation and dropout.

Each decoder takes the GRU features as input. For every video frame the SP decoder produces an eight-frame-long estimate $\widehat{W}_{se} \in {\rm I\!R}^{60\times8}$ of the normalised dimensionality-reduced SP, $W_{se}$, through three 2-D transposed convolutions (convT2D), each followed by BN, ReLU activation and dropout, and another convT2D followed by ReLU activation.

The VUV decoder consists of a linear layer followed by ReLU activation. A threshold of 0.2 is applied to the output obtaining $\widehat{{W}}_{vuv} \in {\rm I\!R}^{8}$, an estimate of the VUV state, $W_{vuv}$.

The AP decoder has a structure similar to the SP decoder, with a total of three convT2D in this case. Its output, $O_{nap} \in {\rm I\!R}^{5\times8}$, together with $\widehat{{W}}_{vuv}$ is used to get $\widehat{W}_{nap}$, an estimate of $W_{nap} = I_{5,8} - W_{ap}$, where $I_{5,8}$ indicates an all-ones matrix with 5 rows and 8 columns, and $W_{ap}$ is the normalised dimensionality-reduced AP:
\begin{align}
\begin{split}\label{eqn:nap}
   (\widehat{W}_{nap})_i ={}& (O_{nap})_i \odot  \widehat{W}_{vuv} \;\;\; \text{for} \; i \in \{1, \dots, 5\}
\end{split}
\end{align}
where $(A)_i$ indicates the $i$-th row of $A$ and $\odot$ denotes the element-wise product.

The F0 decoder has a linear layer followed by a sigmoid activation function. Its output, $O_{f0} \in {\rm I\!R}^{8}$, is point-wise multiplied with $\widehat{{W}}_{vuv}$ to obtain $\widehat{{W}}_{f0}$, an estimate of the normalised F0, $W_{f0}$:
\begin{align}
\begin{split}\label{eq:f0}
   \widehat{{W}}_{f0} ={}& O_{f0} \odot  \widehat{W}_{vuv}.
\end{split}
\end{align}

Finally, the VSR decoder, consisting of a linear and a softmax layers, outputs a CTC character that will be used to predict the text transcription of the utterance.

The system is trained to minimise the following loss:
\begin{multline}J =
\frac{{\lambda}_{1}}{{\lambda}} J_{se} +
\frac{{\lambda}_{2}}{{\lambda}}  J_{nap} +
\frac{{\lambda}_{3}}{{\lambda}}  J_{f0} +
\frac{{\lambda}_{4}}{{\lambda}}  J_{vuv} +
\frac{{\lambda}_{5}}{{\lambda}}  J_{vsr}
\label{eqn:loss}\end{multline}
where ${\lambda}_{1}=600$, ${\lambda}_{2}=50$, ${\lambda}_{3}=10$, ${\lambda}_{4}=10$, ${\lambda}_{5}=1$, ${\lambda}=\sum_{i=1}^{5}{{\lambda}_{i}}$ and:
\begin{itemize}
	\item $J_{se}$: mean squared error (MSE) between $W_{se}$ and $\widehat{W}_{se}$.
	\item $J_{nap}$: MSE between $W_{nap}$ and $\widehat{W}_{nap}$.
	\item $J_{f0}$: MSE between ${W}_{f0}$ and $\widehat{W}_{f0}$.
	\item $J_{vuv}$: MSE between ${W}_{vuv}$ and $\widehat{W}_{vuv}$.
	\item $J_{vsr}$: CTC loss \cite{graves2006connectionist} between the target text transcription and the estimated one.
\end{itemize}

Details regarding architecture and training hyperparameters can be found in Table \ref{tab:architecture}.

\begin{table}
\caption{Architecture and training hyperparameters. Activation, batch normalisation and dropout omitted for brevity.}
\centering
\resizebox{0.45\textwidth}{!}{
\begin{tabular}{c | c c c c c }
\toprule
\multicolumn{6}{c}{\textbf{Input Size}}\\
\bottomrule
\multicolumn{6}{c}{$B{\times}S{\times}C{\times}F{\times}H{\times}W$}\\
\toprule
\multicolumn{6}{c}{\textbf{Video Encoder}}\\
\bottomrule
{Layer} & {\begin{tabular}{@{}c@{}}Input \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Output \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Kernel \\ Size\end{tabular}}   & Stride   & Padding    \\
\midrule
Conv3D &        3        &      64      &   (7,4,4)      &     (1,2,2)          &   (0,1,1)      \\
Conv3D &      64        &    128      &  (1,4,4)       &     (1,2,2)          &   (0,1,1)      \\
Conv3D &    128        &    256      &   (1,4,4)       &    (1,$d_1$,2)  &   (0,1,1)      \\
Conv3D &    256        &    512      &   (1,4,4)       &     (1,2,2)         &   (0,1,1)      \\
Conv3D &    512        &    128      &   (1,$d_2$,6) &     (1,1,1)       &   (0,0,0)      \\
\toprule
\multicolumn{6}{c}{\textbf{Recursive Module}}\\
\bottomrule
{Layer} & \multicolumn{2}{c}{Input Size}   & \multicolumn{2}{c}{Hidden Size}      \\
\midrule
GRU &        \multicolumn{2}{c}{128}        &      \multicolumn{2}{c}{128}               \\
\toprule
\multicolumn{6}{c}{\textbf{Spectral Envelope (SP) Decoder}}\\
\bottomrule
{Layer} & {\begin{tabular}{@{}c@{}}Input \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Output \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Kernel \\ Size\end{tabular}}   & Stride   & Padding   \\
\midrule
ConvT2D &    128        &    256      &   (1,6)       &     (1,1) &     (0,0)   \\
ConvT2D &    256        &    128      &   (2,4)       &     (1,2) &     (0,0)     \\
ConvT2D &    128        &      64      &   (4,4)       &      (1,2) &     (0,0)     \\
ConvT2D &      64        &        1      &   (4,2)       &     (1,2)  &     (0,0)    \\
\toprule
\multicolumn{6}{c}{\textbf{Aperiodic Parameter (AP) Decoder}}\\
\bottomrule
{Layer} & {\begin{tabular}{@{}c@{}}Input \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Output \\ Channels\end{tabular}}   & {\begin{tabular}{@{}c@{}}Kernel \\ Size\end{tabular}}   &  Stride   & Padding  \\
\midrule
ConvT2D &    128        &    128      &   (4,1)  &     (1,1) &     (0,0)         \\
ConvT2D &    128        &      64      &   (3,3)  &     (1,1) &     (0,0)         \\
ConvT2D &      64        &        1      &   (3,3)  &     (1,1) &     (0,0)          \\
\toprule
\multicolumn{6}{c}{\textbf{Voiced-Unvoiced (VUV) Decoder}}\\
\bottomrule
{Layer} & \multicolumn{2}{c}{Input Size}   & \multicolumn{2}{c}{Output Size}      \\
\midrule
Linear &        \multicolumn{2}{c}{128}        &      \multicolumn{2}{c}{8$^{\mathrm{a}}$}               \\
\toprule
\multicolumn{6}{c}{\textbf{Fundamental Frequency (F0) Decoder}}\\
\bottomrule
{Layer} & \multicolumn{2}{c}{Input Size}   & \multicolumn{2}{c}{Output Size}      \\
\midrule
Linear &        \multicolumn{2}{c}{128}        &      \multicolumn{2}{c}{8$^{\mathrm{a}}$}               \\
\toprule
\multicolumn{6}{c}{\textbf{Visual Speech Recognition (VSR) Decoder}}\\
\bottomrule
{Layer} & \multicolumn{2}{c}{Input Size}   & \multicolumn{2}{c}{Output Size}      \\
\midrule
Linear &        \multicolumn{2}{c}{128}        &      \multicolumn{2}{c}{28$^{\mathrm{b}}$}               \\
\toprule
\multicolumn{6}{c}{\textbf{Extra Information}}\\
\bottomrule
\multicolumn{6}{l}{The system is implemented in Pytorch \cite{NEURIPS2019_9015} and trained for $N$}\\
\multicolumn{6}{l}{iterations using the Adam optimizer \cite{kingma2014adam} with a learning rate}\\
\multicolumn{6}{l}{of 0.0001, $\beta_1$=0.5 and $\beta_2$=0.9. The model that performs the}\\
\multicolumn{6}{l}{best in terms of PESQ on the validation set is used for testing.}\\
\multicolumn{6}{l}{$S$=75 (sequence length). $C$=3 (image channels).}\\
\multicolumn{6}{l}{$F$=7 (consecutive video frames). $W$=96 (video frame width).}\\
\multicolumn{6}{l}{If the full face is used as input:}\\
\multicolumn{6}{l}{$B$=16 (batch size). $H$=128 (video frame height). $d_1$=3. $d_2$=5.}\\
\multicolumn{6}{l}{If only the mouth is used as input:}\\
\multicolumn{6}{l}{$B$=24 (batch size). $H$=64 (video frame height). $d_1$=2. $d_2$=4.}\\
\multicolumn{6}{l}{In the speaker dependent case, the dropout probability of each}\\
\multicolumn{6}{l}{dropout layer is $p_d$=0.2. $N$=300000.}\\
\multicolumn{6}{l}{In the speaker independent case, $p_d$=0.5 for the video encoder}\\
\multicolumn{6}{l}{and the GRU, and $p_d$=0.2 for the rest. $N$=185000.}\\
\bottomrule
\multicolumn{6}{l}{$^{\mathrm{a}}$Eight is the number of the output audio frames corresponding}\\
\multicolumn{6}{l}{to the video frame used as input (together with its context).}\\
\multicolumn{6}{l}{$^{\mathrm{b}}$The 28 CTC characters consist of the 26 letters of the English}\\
\multicolumn{6}{l}{alphabet, one space character and one blank token.}\\
\end{tabular}}
\label{tab:architecture}
\end{table}

\begin{table}
\caption{Systems used in this study.}
\centering
\resizebox{0.45\textwidth}{!}{%
\begin{tabular}{l c c}
\toprule
\multicolumn{1}{c}{}   & \multicolumn{2}{c}{Input}\\
\midrule
 & Mouth & Face \\
\midrule
w/o VSR Decoder & vid2voc-M & vid2voc-F  \\
w/ VSR Decoder & vid2voc-M-VSR & vid2voc-F-VSR  \\
\bottomrule
\end{tabular}}
\label{tab:systems}
\end{table}

\subsection{Waveform Reconstruction and Lipreading}

The network outputs are used to reconstruct the speech waveform with the WORLD synthesis algorithm \cite{morise2016world} and to get a text transcription adopting the best path CTC decoding scheme~\cite{graves2006connectionist}.

\subsection{Evaluation Metrics}

The system is evaluated in terms of perceptual evaluation of speech quality (PESQ) \cite{rix2001perceptual} and extended short-time objective intelligibility (ESTOI) \cite{jensen2016algorithm}, two of the most used measures that provide estimates of speech quality and speech intelligibility, respectively. PESQ scores are in the range from $-0.5$ to $4.5$ and ESTOI scores practically lie between $0$ and $1$. In both cases, higher values correspond to better performance.

For the systems having the VSR decoder, we also provide the word error rate (WER), a standard metric for automatic speech recognition systems. In this case, lower values correspond to better performance.

\section{Results and Discussion}

As shown in Table \ref{tab:systems}, four systems are trained based on the input (mouth or full face) and the presence of the VSR decoder (only speech synthesis or speech synthesis and VSR).

The systems are compared with the recently proposed GAN-based approach in \cite{Vougioukas2019}.
As an additional baseline, we also report the PESQ score for \cite{akbari2018lip2audspec}, since this method, which makes use of bottleneck features extracted from auditory spectrograms,  outperforms \cite{Vougioukas2019} in terms of estimated speech quality for the speaker dependent case.

\begin{table}
\caption{Results for the speaker dependent and the speaker independent cases. Best performance (except WORLD) in bold.}
\centering
\resizebox{0.45\textwidth}{!}{%
\begin{tabular}{l | c c c | c c c}
\toprule
\multicolumn{1}{c}{}   & \multicolumn{3}{c}{Speaker Dependent} & \multicolumn{3}{c}{Speaker Independent}\\
\midrule
Mean Scores & PESQ $\uparrow$ & ESTOI $\uparrow$ & WER $\downarrow$ & PESQ $\uparrow$ & ESTOI $\uparrow$ & WER $\downarrow$  \\
\midrule
Approach in \cite{akbari2018lip2audspec}$^{\mathrm{a}}$  & 1.82 & - & - & - & - & - \\
Approach in \cite{Vougioukas2019} & 1.71 & 0.329 & - & 1.24 & 0.198 & - \\
vid2voc-M & 1.89 & 0.448 & - &  1.20 & 0.214 & - \rule{0pt}{3ex} \\
vid2voc-M-VSR & \textbf{1.90} & \textbf{0.455} & 15.1\% & 1.23 & \textbf{0.227} & \textbf{51.6\%} \\
vid2voc-F & 1.85 & 0.439 & - & 1.19 & 0.202 & -  \rule{0pt}{3ex}\\
vid2voc-F-VSR & 1.88 & 0.447 & \textbf{14.4\%} & \textbf{1.25} & 0.210 & 69.3\% \\
WORLD$^{\mathrm{b}}$ & 3.06 & 0.759 & - & 3.03 & 0.759 & - \rule{0pt}{3ex}\\
\bottomrule
\multicolumn{7}{l}{$^{\mathrm{a}}$Value taken from the experiments in \cite{Vougioukas2019}.}\\
\multicolumn{7}{l}{$^{\mathrm{b}}$WORLD indicates the reconstruction retrieved from the vocoder features of the}\\
\multicolumn{7}{l}{clean speech signals and it is a performance upper bound of our systems.}
\end{tabular}}
\label{tab:results}
\end{table}

\subsection{Speaker Dependent Case}

Table \ref{tab:results} (left part) shows the speaker dependent results. We observe that our models outperform the approach in \cite{Vougioukas2019} in terms of both PESQ and ESTOI by a considerable margin. Vougioukas et al. \cite{Vougioukas2019} mention that their  system produces low-power hum artefacts that affect the performance. They tried to solve the issue by applying average filtering to the output of their network, experiencing a rise of the PESQ score from 1.71 to 1.80 (not shown in Table \ref{tab:results}), comparable to \cite{akbari2018lip2audspec}, but still appreciably lower than the results we achieve. However, this filtering negatively affected the intelligibility of the produced speech signals, and was not used in the final system.

Among the systems we developed (cf. Table \ref{tab:systems}), we observe that including the VSR decoder in the pipeline is beneficial for the speech reconstruction task (see Table \ref{tab:results}). Moreover, the use of the mouth as input not only is sufficient to synthesise speech, but it also allows to achieve higher estimated speech quality and intelligibility if compared to the models that use the whole face of the speaker as input. This might be explained by the fact that handling an input with a larger dimensionality is harder if we want to keep roughly the same deep architecture with a similar number of parameters. However, when the whole face is used as input, the WER is slightly lower, indicating that there might be a performance trade-off between VSR and speech reconstruction that should be further investigated in future work in relation with other multi-task learning techniques.

\begin{figure}
	\centering
		\includegraphics[scale=.43]{plot.pdf}
	\caption{Results of the vid2voc-M-VSR models for the speaker dependent (SD) and the speaker independent (SI) cases. Each marker indicates the mean score of a speaker.}
	\label{fig:plot}
\end{figure}

\subsection{Speaker Independent Case}

Regarding the speaker independent scenario (cf. right part of Table \ref{tab:results}), we observe that the performance gap between the approach in \cite{Vougioukas2019} and our systems is not as large as for the speaker dependent case. Although our models appear to perform slightly better than  \cite{Vougioukas2019} in terms of ESTOI, the PESQ scores are similar. This can be explained by the fact that some speech characteristics, e.g. F0, cannot be easily estimated for unseen speakers. Since it is reasonable to think that people having similar facial characteristics (e.g. due to gender, age etc.) have similar speech characteristics (cf. \cite{oh2019speech2face}, where the face of a person was predicted from a speech signal), we expect that training a network with a dataset that includes more speakers might be beneficial: such a network can produce an average voice of speakers from the training set that share similar facial traits with an unseen talking face.

Among the systems we developed, the presence of the VSR decoder still gives an advantage for speech reconstruction. Unlike the speaker dependent case, the WER for the model that uses the whole face as input is higher than the system using only the mouth. This is due to the early stopping technique that we adopt, which tends to favour speech reconstruction over VSR, indicating again the trade-off between these two tasks.

Finally, Figure \ref{fig:plot} shows the results for the vid2voc-M-VSR models by speaker. We can see that the spread of the scores is much higher for the speaker independent case in particular for WER. This is in line with the observations reported in \cite{Vougioukas2019}, suggesting the different performance between the estimated speech of subjects whose facial traits substantially differ from the speakers in the training set and the others.

\section{Conclusion}

In this study, we reconstructed speech from silent videos using a deep model that estimates WORLD vocoder features. We tested our approach in both speaker dependent and speaker independent scenarios. In both cases, we were able to obtain speech signals with estimated speech quality and intelligibility generally higher if compared to a recently proposed GAN-based approach. In addition, we designed our system to simultaneously perform visual speech recognition by using a decoder that estimates CTC characters from a given video sequence.

Future work includes: (a) the adoption of self-paced multi-task learning techniques; (b) the improvement of the visual speech recognition performance, e.g. with a beam search decoding scheme; (c) the design of a system that can generalise well to unseen speakers in noncontrolled environments.

\section{Acknowledgment}

The authors would like to thank Konstantinos Vougioukas, Stavros Petridis, Pritish Chandna and Merlijn Blaauw.

This research is partially funded by: the William Demant Foundation; the TROMPA H2020 project (770376); the Spanish Ministry of Economy and Competitiveness under the Mar{\'i}a de Maeztu Units of Excellence Program (MDM-2015-0502) and the Social European Funds; the MICINN/FEDER UE project (PGC2018-098625-B-I00); the H2020-MSCA-RISE-2017 project (777826 NoMADS).

\bibliographystyle{IEEEtran}


\section{Introduction and statement of the Problem}

There are algebraic and variational approaches in the spline theory. In the
algebraic approach splines are considered as some smooth piecewise polynomial
functions. In the variational approach splines are elements of Hilbert or Banach
spaces minimizing certain functionals. The first spline functions were constructed
from pieces of cubic polynomials. After that, this construction was modified, the
degree of polynomials increased. The theory of splines based on variational methods studied and
developed, for example, by J.H.Ahlberg et al. \cite{Ahlb67},    C. de Boor \cite{deBoor63}, A.R.Hayotov, G.V.Milovanivi\'{c} and Kh.M.Shadimetov \cite{HayShadMil14},  I.J.Schoenberg \cite{Schoen64},
L.L.Schumaker \cite{Schum81}, S.L.Sobolev \cite{Sob61a}, V.A.Vasilenko \cite{Vas83} and others.

The present work is also devoted to the variational method of construction of optimal interpolation formulas.

Assume, we are given the table of the values $\varphi(x_\beta)$, $\beta=0,1,...,N$ of a
function $\varphi$ at the points $x_{\beta}\in [0,1]$. It is required to approximate the function $\varphi$ by
another more simple function $P_\varphi$, i.e.
\begin{equation}\label{(1)}
\varphi (x) \cong P_\varphi  (x) = \sum\limits_{\beta = 0}^N
{C_\beta(x) \cdot \varphi (x_\beta )}
\end{equation}
which satisfies the following interpolation conditions
$$
\varphi(x_{\beta})=P_{\varphi}(x_{\beta}), \ \beta=0,1,...,N.
$$
Here $C_\beta (x)$ and $x_\beta$ ($\in [0,1]$) are
\emph{the coefficients} and \emph{the nodes} of the interpolation formula (\ref{(1)}), respectively.
By ${K_2}(P_2)$ we denote the class of all functions $\varphi$ defined on [0,1] which posses an absolutely continuous first derivative on [0,1] and whose second derivative is in $L_2(0,1)$. The class $K_2(P_2)$ under the pseudo-inner product
\begin{equation}\nonumber
\langle {\varphi,\psi} \rangle  =
\int\limits_0^1 \left(\varphi ''(x)+\omega^2\varphi(x)\right)\left(\psi_{\ell}''(x)+\omega^2\psi_{\ell}(x)\right)dx.
\end{equation}
is a Hilbert space if we identify functions that differ by $\sin(\omega x)$ and $\cos(\omega x)$, where $\omega\neq 0$.
Here we consider the norm
\begin{equation}\nonumber
\left\| {\varphi|K_2(P_2)} \right\| = \left\{
{\int\limits_0^1 {\left( \varphi''(x)+\omega^2\varphi(x)
\right)} ^2 \d x} \right\}^{1/2}.
\end{equation}
For the fixed $z\in [0,1]$ the error
\begin{equation}\label{(3)}
(\ell,\varphi)=\int\limits_{-\infty}^{\infty}\ell(x)\varphi(x)dx=\varphi(z)-P_{\varphi}(z)=\varphi(z)-\sum\limits_{\beta = 0}^N
{C_\beta(z) \cdot \varphi (x_\beta )}
\end{equation}
of the interpolation formula (\ref{(1)}) is a linear functional.\\
Here, in fixed $z\in [0,1]$,
$$
\ell(x,z)=\delta(x-z)-\sum_{\beta=0}^NC_\beta(z)\delta(x-x_\beta),
$$
and is the error functional of the
interpolation formula (\ref{(1)}) belongs to the space $K_2^*(P_2)$. Here $K_2^*(P_2)$ is the
conjugate space to the space $K_2(P_2)$, $\delta$ is  Dirac's delta-function.

By the Cauchy-Schwarz inequality the absolute value of the error (\ref{(3)}) is
estimated as follows
$$
|(\ell,\varphi)|\leq
\|\varphi|{K_2(P_2)}\|\cdot
\|\ell|{K_2^{*}(P_2)}\|,
$$
where
$$
\left\| {\ell|{K_2^{*}(P_2)} } \right\| = \mathop {\sup }\limits_{\varphi,\ \|\varphi \| \neq 0}
\frac{|(\ell,\varphi)|}{\|\varphi|{K_2(P_2)}\|}.
$$

Therefore, in order to estimate the error of the interpolation formula (\ref{(1)}) on
functions of the space $K_2(P_2)$ it is required to find the norm of the error functional $\ell$
in the conjugate space $K_2^{*}(P_2)$. That is we get the following problem.

\medskip

\textbf{Problem 1.1.} \textit{Find the norm of the error functional $\ell$ of the interpolation
formula (\ref{(1)}) in the space $K_2^{*}(P_2)$.}

\medskip

It is clear that the norm of the error functional $\ell$ depends on the coefficients $C_\beta(z)$
and the nodes $x_\beta$. The problem of minimization of the quantity $\|\ell\|$ by
coefficients $C_\beta(z)$ is the linear problem and by nodes $x_\beta$ is, in general, nonlinear
and complicated problem. We consider the problem of minimization of the
quantity $\|\ell\|$ by coefficients $C_\beta(z)$ when nodes $x_\beta$ are fixed.

The coefficients $\mathring{C}_{\beta}(z)$ (if there exist) satisfying the  equality
\begin{equation}\label{(4)}
\left\| \mathring{\ell}|{K_2^{*}(P_2)}\right\| = \mathop {\inf }\limits_{C_\beta(z) }  \left\|\ell|{K_2^{*}(P_2)} \right\|
\end{equation}
are called \textit{the optimal coefficients} and corresponding interpolation formula
$
\mathring{P}_{\varphi}(z)=\sum\limits_{\beta=0}^N\mathring{C}_{\beta}(z)\varphi(x_\beta)
$
is called \textit{the optimal interpolation formula} in the space $K_2^{*}(P_2)$.
Therefore, for construction of the interpolation formula we should solve the next problem.

\medskip

\textbf{Problem 1.2.} \textit{Find the coefficients $\mathring{C}_{\beta}(z)$ which satisfy equality (\ref{(4)}) when the
nodes $x_\beta$ are fixed.}

\medskip

It should be noted that Problems 1.1 and 1.2 were solved in \cite{BabHay18} in the Hilbert space $W_2^{(m,m-1)}$.
There the optimal interpolation formulas, which are exact for any polynomial of degree $m-2$ and for the function $\exp(-x)$, were obtained.

The rest of the paper is organized as follows. In Section 2, using the extremal function, the norm of the error functional is found.
Existence and uniqueness of the optimal interpolation formula of the form (\ref{(1)}) is discussed in Section 3. Section 4 is devoted to construction of the optimal interpolation formula. Finally, in Section 5 some numerical results are presented.

\section{The extremal function and the norm of the error functional $\ell$ }

Here we find explicit form of
the norm of the error functional $\ell$.

For finding the explicit form of the norm of the error functional
$\ell$ in the space $K_2(P_2)$ we use its extremal function which was introduced by Sobolev
\cite{Sob61a,Sob74}. The function  $\psi_{\ell}$ from
${K_2}(P_2)$ space is called \emph{the extremal
function} for the error functional $\ell$ if the following
equality is fulfilled
$$
\left( {\ell,\psi_{\ell}} \right) = \left\| {\ell\left|
{{K_2^*(P_2)}} \right.} \right\| \cdot \left\|
{\psi_{\ell}\left|{{K_2(P_2)}} \right.} \right\|.
$$
According to the Riesz theorem any linear continuous functional
$\ell$ in a Hilbert space is represented in the form of a
inner product. So, in our case, for any function $\varphi$ from ${K_2}(P_2)$
space we have
\begin{equation} \label{(5)}
\left( {\ell,\varphi} \right) = \langle{\psi _\ell,\varphi} \rangle.
\end{equation}
Here $\psi_\ell$ is the function from
${K_2}(P_2)$ is defined  uniquely by functional $\ell$ and is the extremal function.

It is easy to see from (\ref{(5)}) that the error functional $\ell$, defined on the space $K_2(P_2)$, satisfies the following equalities
\begin{eqnarray}
&&(\ell,\sin(\omega x))=0, \label{(6)}\\
&&(\ell,\cos(\omega x))=0.    \label{(7)}
\end{eqnarray}
The equalities (\ref{(6)}) and (\ref{(7)}) mean that our interpolation formula is exact for the functions $\sin(\omega x)$ and $\cos(\omega x)$.

The equation (\ref{(5)}) was solved in \cite{HayShadMil14} and for the extremal function $\psi_{\ell}$ was obtained the following
expression
\begin{equation}\label{(8)}
\psi_{\ell}(x)=(\ell*G_2)(x)+d_1\sin(\omega x)+d_2\cos(\omega x),
\end{equation}
where
\begin{equation}\label{(9)}
G_2(x)=\frac{\mathrm{sgn}x}{4\omega^3}[\sin(\omega x)-\omega x\cos(\omega x)],
\end{equation}
$*$ is the operation of convolution which for the functions $f$ and $g$ is defined as follows
\begin{equation}\label{(10)}
(f*g)(x)=\int\limits_{-\infty}^{\infty}f(x-y)g(y)\d y=\int\limits_{-\infty}^{\infty}f(y)g(x-y)\d y.
\end{equation}
Now we obtain the norm of the error functional
$\ell$. Since the space ${K_2}(P_2)$ is the Hilbert
space then by the Riesz theorem we have
\begin{equation}\label{(11)}
 \left( {\ell,\psi_\ell} \right) = \|\ell\| \cdot \| \psi _\ell\| =\| \ell\|^2 .
\end{equation}
Hence, using (\ref{(8)}) and (\ref{(9)}), taking
into account (\ref{(10)}) and (\ref{(11)}),  we get
{\small
\begin{eqnarray*}
\|\ell\|^2&=&
(\ell,\psi_{\ell})=\int\limits_{-\infty}^{\infty}\ell(x,z)\psi_{\ell}(x)\d x\\
&=&\int\limits_{-\infty}^{\infty}\left(\delta(x-z)-\sum\limits_{\beta=0}^NC_{\beta}(z)\delta(x-x_{\beta})\right)\\
&&\times \left(G_2(x-z)-\sum\limits_{\beta=0}^NC_{\beta}(z)G_2(x-x_{\beta})\right)\d x.
\end{eqnarray*}}
Hence, keeping in mind that $G_2(x)$, defined by (\ref{(9)}), is the even function, we
have
\begin{equation}
\|\ell\|^2=(-1)^m\left(\sum\limits_{\beta=0}^N\sum\limits_{\gamma=0}^N
C_{\beta}(z)C_{\gamma}(z)G_2(x_\beta-x_\gamma)-2\sum\limits_{\beta=0}^NC_{\beta}(z)G_2(z-x_{\beta})
\right).\label{(12)}
\end{equation}
Thus, Problem 1.1 is solved.

Further, we solve Problem 1.2.

\section{Existence and uniqueness of the optimal interpolation formula}

Assume that the nodes $x_\beta$  of the interpolation formula
(\ref{(1)}) are fixed. The error functional (\ref{(3)})
satisfies the conditions (\ref{(6)}) and (\ref{(7)}). The
norm of the error functional $\ell$ is  a multivariable
function with respect to the coefficients $C_\beta(z)$ $(\beta  =
\overline {0,N} )$. For finding the point of the conditional
minimum of the expression (\ref{(12)}) under the conditions
(\ref{(6)}) and (\ref{(7)}) we apply the Lagrange
method.

Consider the function
\begin{eqnarray*}
&&\Psi (C_0(z),C_1(z),...,C_N(z),d_1,d_2)\\
&&\qquad=\left\| \ell \right\|^2  - 2\Bigg(d_1(z)(\ell,\sin(\omega x))+d_2(z)(\ell,\cos(\omega x))\Bigg).
\end{eqnarray*}
Equating to 0 the partial derivatives of the function $\Psi$
by $C_\beta(z)$ $(\beta =\overline{0,N})$, $d_1(z)$ and $d_2(z)$, we get the following
system of $N+3$ linear equations of $N+3$ unknowns
\begin{eqnarray}
&&\sum\limits_{\gamma=0}^N C_\gamma(z) G_2(x_\beta   - x_\gamma  )
+ d_1(z)\sin(\omega x_\beta)+d_2(z)\cos(\omega x_\beta) = G_2(z-x_{\beta}),\nonumber\\
&&\beta  =0,1,...,N,\    \label{(13)} \\
&& \sum\limits_{\gamma=0}^N C_\gamma(z)\sin(\omega x_\gamma)=\sin(\omega z),\ \label{(14)}\\
&&\sum\limits_{\gamma=0}^NC_\gamma(z)\cos(\omega x_\gamma)=\cos(\omega z),\label{(15)}
\end{eqnarray}
where $G_2(x)$ is defined by equality (\ref{(9)}).

The system (\ref{(13)})-(\ref{(15)}) has a unique solution
and this solution gives the minimum to $\left\|\ell\right\|^2
$ under the conditions (\ref{(14)}) and (\ref{(15)}).

The uniqueness of the solution of the system
(\ref{(13)})--(\ref{(15)}) is proved as the uniqueness of the solution of the system (24)--(26) of the work
\cite{SobVas}.

Therefore, in fixed values of the nodes $x_\beta$ the square of
the norm of the error functional $\ell$, being quadratic
function of the coefficients $C_\beta(z)$, has a unique minimum in
some concrete value $C_\beta(z)= \mathring{C}_{\beta }(z)$.

\begin{rem}\label{Rem(1)}
It should be noted that by integrating both sides of the system (\ref{(13)})-(\ref{(15)}) by $z$ from 0 to 1
we get the system (3.1)-(3.1) of the work \cite{HayShadMil14}. This means that by integrating the optimal interpolation formula (\ref{(1)}) in  the space $K_2(P_2)$ we get the optimal quadrature formula of the form (\ref{(1)}) in the same space (see \cite{HayShadMil14}).
\end{rem}

\begin{rem} \label{Rem(2)}
It is clear from the system (\ref{(13)})-(\ref{(15)}) that for the optimal coefficients the following are true
\begin{eqnarray*}
\mathring{C}_{\beta}(h\gamma)&=&\left\{
\begin{array}{ll}
1, & \gamma =\beta,\\
0,& \gamma\neq \beta,
\end{array}
\right. \gamma=0,1,...,N,\ \ \beta=0,1,...,N.
\end{eqnarray*}
\end{rem}

Below for convenience the optimal coefficients $ \mathring{C}_{\beta}(z)$ we remain as $C_\beta(z)$.
\section{The algorithm for computation of coefficients of the optimal interpolation formula}

In the present section we give the algorithm for solution of the
system (\ref{(13)})-(\ref{(15)}). Below mainly is used the concept of discrete
argument functions and operations on them. The theory of discrete
argument functions is given, for instance, in \cite{Sob74,SobVas}. For
completeness we give some definitions about functions of discrete
argument.

Assume that the nodes $x_\beta$ are equal spaced, i.e. $x_\beta=
h\beta,$ $h = {1 \over N}$, $N = 1,2,...$.

\begin{defn}\label{def1}
  The function $\varphi (h\beta )$ is a
\emph{discrete argument function} (or \emph{discrete function}) if it is given on some set of
integer values of $\beta$.
\end{defn}

\begin{defn}\label{def2}
\emph{The inner product} of two discrete
functions $\varphi(h\beta )$ and $\psi (h\beta )$ is given by
$$
\left[ {\varphi(h\beta),\psi(h\beta) } \right] =
\sum\limits_{\beta  =  - \infty }^\infty  {\varphi (h\beta ) \cdot
\psi (h\beta )},
$$
if the series on the right hand side converges absolutely.
\end{defn}

\begin{defn}\label{def3}
\textit{The convolution} of two functions
$\varphi(h\beta )$ and $\psi (h\beta )$ is the inner product
$$
\varphi (h\beta )*\psi (h\beta ) = \left[ {\varphi (h\gamma ),\psi
(h\beta  - h\gamma )} \right] = \sum\limits_{\gamma  =  - \infty
}^\infty  {\varphi (h\gamma ) \cdot \psi (h\beta  - h\gamma )}.
$$
\end{defn}

Now we turn on to our problem.

Suppose that $C_\beta(z)=0$  when $\beta  < 0$ and $\beta  > N$.
Thus we have the following problem.

\medskip

\textbf{Problem 4.1} {\it
Find the discrete functions $C_\beta(z)$, $\beta=0,1,...,N$,  $d_1(z)$ and $d_2(z)$ which
satisfy the system (\ref{(14)})-(\ref{(15)}).}

\medskip

Further we investigate Problem 4.1 which is equivalent to Problem 1.2.
Instead of $C_\beta(z)$ we introduce the following functions
\begin{eqnarray}
v_2(h\beta ) &=& G_2(h\beta )*C_\beta(z),\nonumber\\
u_2(h\beta)&=&v_2\left({h\beta}\right)+ d_1(z)\sin(\omega h\beta)+d_2\cos(\omega h\beta).\label{(17)}
\end{eqnarray}

Now we should express the coefficients
$C_\beta(z)$ by the function $u(h\beta)$. For this we use the operator $D_2(h\beta)$ which satisfies the equality
\begin{equation}
D_2(h\beta )*G_2(h\beta ) =\delta(h\beta) ,\label{(18)}
\end{equation}
where $\delta(h\beta)$ is equal to 0 when $\beta \ne 0$ and is
equal to 1 when $\beta  = 0$, i.e. $\delta(h\beta)$ is the
discrete delta-function.

In \cite{Hay04a} the
operator $D_2(h\beta )$  which satisfies
equation (\ref{(18)}) is constructed and its some properties
are studied.

The following theorems are proved in \cite{Hay04a}.

\begin{thm}\label{THM4.1}
The discrete analogue of the differential operator ${{d^{4}}\over {d x^{4} }}+2\omega^2{{d^{2}} \over {d x^{2} }}+\omega^4$ satisfying the
equation (\ref{(18)}) has the form
\begin{equation}
D_2(h\beta)=p\left\{
\begin{array}{ll}
A\lambda^{|\beta|-1}, &|\beta|\geq 2,\\
1+A,& |\beta|=1,\\
C+\frac{A}{\lambda}, & \beta=0.
\end{array}
\right. \label{(19)}
\end{equation}
where
\begin{equation}
\displaystyle
\begin{array}{l}
p=\frac{2\omega^3}{\sin(\omega h)-\omega h\cos(\omega h)},\nonumber\\
C=\frac{2h\omega\cos(2h\omega)-\sin(2h\omega)}{\sin(h\omega)-h\omega\cos(h\omega)},\nonumber\\
A=\frac{(2h\omega)^2\sin^4(h\omega)\cdot\lambda^2}{(\lambda^2-1)(\sin(h\omega)-h\omega\cos(h\omega))^2},\nonumber\\
\lambda=\frac{2h\omega-\sin(2h\omega)-2\sin(h\omega)\cdot\sqrt{(h\omega)^2-\sin^2(h\omega)}}{2(h\omega\cos(h\omega)-\sin(h\omega))},\  |\lambda|<1.
\end{array}
\end{equation}

\end{thm}
\begin{thm}\label{THM4.2}
The discrete analogue $D_2(h\beta)$ of the
differential operator ${{d^{4}} \over {d x^{4} }}+2\omega^2{{d^{2}}\over {d x^{2} }}+\omega^4$ satisfies the following equalities

1) $D_2(h\beta)*\sin(h\omega\beta)=0,$

2) $D_2(h\beta)*\cos(h\omega\beta)=0,$

3) $D_2(h\beta)*(h\omega\beta)\cos(h\omega\beta)=0,$

4) $D_2(h\beta)*(h\omega\beta)\sin(h\omega\beta)=0.$\\
\end{thm}

Then taking into account (\ref{(18)}),  (\ref{(19)}), using
Theorems \ref{THM4.1} and \ref{THM4.2}, for optimal coefficients we
have
\begin{equation}
C_\beta(z)=D_2(h\beta )*u_2(h\beta ).\label{(28)}
\end{equation}

Thus if we find the function $u_2(h\beta )$ then the optimal
coefficients $C_{\beta}(z)$ will be found from equality
(\ref{(28)}).

In order to calculate the convolution (\ref{(28)}) it is required to
find the representation of the function $u_2(h\beta )$ at all
integer values of $\beta$. From equality (\ref{(17)}) we get that
$u_2(h\beta ) = G_2(z-h\beta )$ when $h\beta \in [0,1]$. Now we find the representation of the function $u_2(h\beta )$ when
$\beta < 0$ and $\beta>N$.

Since $C_\beta(z)= 0$ when $h\beta \notin [0,1]$ then
$$
C_\beta(z)   = D_2(h\beta )*u_2(h\beta ) =
0,\,\,\,\, h\beta \notin [0,1].
$$

Now we calculate the convolution $v_2(h\beta ) =G_2(h\beta)*C_\beta(z)$ when $h\beta \notin [0,1]$.

Suppose $\beta  \leq 0$ and $\beta\geq N$ then taking into account equalities
(\ref{(9)}), (\ref{(14)}) and (\ref{(15)}), we have
\begin{equation}
u_2(h\beta)=\left\{
\begin{array}{ll}
d_1^-(z)\sin(\omega h\beta)+d_2^-(z)\cos(\omega h \beta)+\frac{h\beta}{4\omega^2}\cos(\omega (h\beta-z)),& \beta\leq 0,\\[2mm]
G_2(z-h\beta),& 0\leq \beta\leq N,\\[2mm]
d_1^+(z)\sin(\omega h\beta)+d_2^+(z)\cos(\omega h\beta)-\frac{h\beta}{4\omega^2}\cos(\omega(h\beta-z)),& \beta\geq N,\\
\end{array}
\right. \label{(29)}
\end{equation}
where $d_1^-(z),\ d_1^+(z),\ d_2^-(z)$ and $d_2^+(z)$ are unknowns.

From (\ref{(29)}) when $\beta=0$ and $\beta=N$ we get
\begin{eqnarray}
&&d_2^-(z)=G_2(z),  \ \label{(30)}\\
&&d_2^+(z)=\frac{1}{\cos(\omega)}\Bigg[G_2(z-1)+\frac{\cos(\omega(1-z))}{4\omega^2}-d_1^+(z)\sin(\omega)\Bigg].\ \label{(31)}
\end{eqnarray}
Thus, putting (\ref{(30)}) and (\ref{(31)}) to (\ref{(29)}) we have the following explicit form of the function $u_2(h\beta)$:
\begin{equation}
u_2(h\beta)=\left\{
\begin{array}{ll}
d_1^-(z)\sin(\omega h\beta)+G_2(z)\cos(\omega h \beta)+\frac{h\beta}{4\omega^2}\cos(\omega (h\beta-z)),& \beta\leq 0,\\[2mm]
G_2(z-h\beta),& 0\leq \beta\leq N,\\[2mm]
d_1^+(z)\frac{\sin(\omega (h\beta-1))}{\cos(\omega)}+\frac{\cos(\omega h\beta)}{\cos(\omega)}\Bigg[G_2(z-1)+\frac{\cos(\omega(1-z) )}{4\omega^2}\Bigg]\\
-\frac{h\beta}{4\omega^2}\cos(\omega (h\beta-z)),& \beta\geq N.\\
\end{array}
\right. \label{(32)}
\end{equation}
In the last expression of the function $u_2(h\beta)$ we have only two unknowns $d_1^-(z)$ and $d_1^+(z)$.

Hence using (\ref{(19)}) and (\ref{(32)}) we get the following problem.\\

\medskip

\textbf{Problem 4.2.} {\it
Find the solution of the equation
\begin{equation}
D_2(h\beta )*u_2(h\beta) = 0,\,\, h\beta  \notin[0,1],\nonumber
\end{equation}
having the form (\ref{(32)}).
Here $d_1^{-}(z)$ and $d_1^{+}(z)$ are unknowns.}

\medskip

Unknowns $d_1^{-}(z)$ and $d_1^{+}(z)$ we find from the equation
\begin{equation}\label{(34)}
D_2(h\beta )*u_2(h\beta)=0
\end{equation}
when $\beta=-1$ and $\beta=N+1$. From the last equation

From the system (\ref{(34)}) in the case $\beta=-1$ and after some simplifications we have the following system of equations
\begin{equation}
A_{11}d_1^-(z)+A_{12}d_1^+(z)=S_1 \ \label{(35)}
\end{equation}
where
\begin{equation}
\begin{array}{ll}
A_{11}=-\sin(2\omega h)-C\sin(\omega h)-\frac{A\sin(\omega h)}{\lambda(\lambda^2-2\lambda\cos(\omega h)+1)},\\ \nonumber
A_{12}=\frac{A\lambda^{N+1}\sin(\omega h)}{\cos(\omega)\left(\lambda^2-2\lambda\cos(\omega h)+1\right)},\\
S_1=\frac{1}{4\omega^2}\Bigg(Ch\cos(\omega(h+z)+2h\cos(\omega(2h+z))-\frac{A\lambda^{N}\cos(\omega(1-z))}{\cos(\omega)})K_1\\
\ \ +\frac{Ah}{\lambda^2}K_2+Ah\lambda^NK_3\Bigg)-AK_4-\frac{A\lambda^NG_2(z-1)}{\cos(\omega)}K_1\\
\ \ -G_2(z)\Bigg(C\cos(\omega h)+1+\cos(2\omega h)+\frac{A}{\lambda^2}K_5\Bigg),\\
K_1=\frac{\lambda(\cos(\omega(1+h))-\lambda\cos(\omega))}{\lambda^2-2\lambda\cos(\omega h)+1},\\
K_2=\frac{\lambda(\lambda^2\cos(\omega(h-z))-2\lambda\cos(\omega z)+\cos(\omega(h+z)))}{(\lambda^2-2\lambda\cos(\omega h)+1)^2},\\
\end{array}
\end{equation}
\begin{equation}
\begin{array}{ll}
K_3=\sum\limits_{\gamma = 1}^\infty\lambda^\gamma(N+\gamma)\cos(\omega (h(N+\gamma)-z)),\\ \nonumber
K_4=\sum\limits_{\gamma = 0}^N\lambda^\gamma G_2(z-h\gamma),\\
K_5=\frac{\lambda(\cos(\omega h)-\lambda)}{\lambda^2-2\lambda\cos(\omega h)+1}.\\
\end{array}
\end{equation}

Now, from (\ref{(34)}) in the case $\beta = N + 1$ doing some calculations we get the next equation
\begin{equation}
A_{21}d_1^-(z)+A_{22}d_1^+(z)=S_2,  \ \label{(36)}
\end{equation}
where
\begin{equation}
\begin{array}{ll}
A_{21}=-\frac{A\lambda^{N+1}\sin(\omega h)}{\lambda^2-2\lambda\cos(\omega h)+1},\\ \nonumber
A_{22}=\frac{1}{\cos(\omega)}\Bigg[C\sin(\omega h)+\sin(2\omega h)+\frac{A\sin(\omega h)}{\lambda(\lambda^2-2\lambda\cos(\omega h)+1)}\Bigg],\\
S_2=\frac{1}{4\omega^2}\Bigg(A\lambda^NhK_2+\cos(\omega(1-z))+Ch(N+1)\cos(\omega(h(N+1)-z))\\
\ \ +h(N+2)\cos(\omega(h(N+2)-z))+\frac{AhK_3}{\lambda^2}-\frac{\cos(w(1-z))K_6}{\cos(\omega)}\Bigg)\\
\ \ -A\lambda^NK_7-G_2(z)A\lambda^NK_5-\frac{G_2(z-1)}{\cos(\omega)}K_6,\\
K_6=\cos(\omega)+C\cos(\omega h(N+1))+\cos(\omega h(N+2))+\frac{AK_1}{\lambda^2},\\
K_7=\sum\limits_{\gamma = 0}^N\lambda^{-\gamma} G_2(z-h\gamma),\\
K_8=\frac{\lambda(\sin(\omega (h+1))-\lambda\sin(\omega))}{\lambda^2-2\lambda\cos(\omega h)+1}.\\ \nonumber
\end{array}
\end{equation}
Then solving the system  (\ref{(35)}), (\ref{(36)}) of equations we get $d_1^-(z)$ and $d_1^+(z)$.
Finally, from (\ref{(28)}) for $\beta=0,1,...,N$ we get the explicit formulas for optimal coefficients as claimed in the following theorem.

\begin{thm}\label{THM4.3}
Coefficients of the optimal interpolation formula (\ref{(1)})
with equal spaced nodes in the space $K_2(P_2)$ have the
following form
\begin{equation}
\begin{array}{ll}
C_{0}(z)=&
p\Bigg[-d_1^-(z)\sin(\omega h)+d_2^-(z)\cos(\omega h)-\frac{1}{4\omega^2}h\cos(\omega(h+z))+C G_2(z)\\\nonumber
&+G_2(z-h)+\frac{A}{\lambda}\Bigg[\sum\limits_{\gamma=0}^N\lambda^\gamma G_2(z-h\gamma)+M_1+\lambda^N N_1\Bigg]\Bigg],\\
C_{\beta}(z)=&
p\Bigg[G_2(z-h(\beta-1))+C G_2(z-h\beta)+G_2(z-h(\beta+1))\\ \nonumber
\end{array}
\end{equation}
\begin{equation}
\begin{array}{ll}
&+\frac{A}{\lambda}\Bigg[\sum\limits_{\gamma=0}^N\lambda^{|\beta-\gamma|} G_2(z-h\gamma)+\lambda^\beta M_1+\lambda^{N-\beta} N_1\Bigg]\Bigg],  \beta=1,2,..,N-1.\\ \nonumber
\end{array}
\end{equation}
\begin{equation}
\begin{array}{ll}\label{(37)}
C_{N}(z)=&p\Bigg[d_1^+(z)\sin(\omega h(N+1))+d_2^+(z)\cos(\omega h(N+1))\\ \nonumber
&-\frac{1}{4\omega^2}h(N+1)\cos(\omega(h(N+1)-z))+C G_2(z-1)\\
&+G_2(z-h(N-1))+\frac{A}{\lambda}\Bigg[\sum\limits_{\gamma=0}^N\lambda^{N-\gamma} G_2(z-h\gamma)+\lambda^N M_1+ N_1\Bigg]\Bigg],\\
\end{array}
\end{equation}
here
\begin{equation}
\begin{array}{ll}
M_1=-d_1^-(z)\frac{\lambda\sin(\omega h)}{\lambda^2-2\lambda\cos(\omega h)+1}+d_2^-(z)K_5-\frac{h}{4\omega^2}K_2, \\ \nonumber
N_1=d_1^+(z)K_8+d_2^+(z)K_1-\frac{h}{4\omega^2} K_3,\\
d_1^-(z)=\frac{S_1\cdot A_{22} - S_2 \cdot A_{12}}{A_{11} \cdot A_{22} - A_{21} \cdot A_{12} },\\
d_1^+(z)=\frac{S_2\cdot A_{11} - S_1 \cdot A_{21}}{A_{11} \cdot A_{22} - A_{21} \cdot A_{12} },\\
d_2^-(z)=G_2(z),\\
d_2^+(z)=\frac{1}{\cos(\omega)}\Bigg[G_2(z-1)+\frac{\cos(\omega(1-z))}{4\omega^2}-d_1^+(z)\sin(\omega)\Bigg].
\end{array}
\end{equation}
\end{thm}

\section{Numerical results}
In this section we give some numerical results.

First, when $N=5$ using Theorem \ref{THM4.3}, we get the graphs of the coefficients of the optimal interpolation formulas
$$
\varphi(z)\cong \mathring{P}_{\varphi}(z)=\sum\limits_{\beta=0}^5\mathring{C}_{\beta}(z)\varphi(h\beta),\ z\in [0,1].
$$ They are presented in Fig 1.
These graphical results confirm Remark \ref{Rem(2)} for the case $N=5$ , i.e. for the optimal coefficients the following hold
$$
\mathring{C}_{\beta}(h\gamma)=\delta_{\beta\gamma},\ \ \beta,\gamma=0,1,...,5,
$$
where $\delta_{\beta\gamma}$ is the Kronecker symbol.

Now, in numerical examples, we interpolate the functions
$$
\varphi_1(z)=z^2,\ \varphi_2(z)=\E^{z}\mbox{ and } \varphi_3(z)=\sin z
$$
by optimal interpolation formulas of the form (\ref{(1)}) in the cases $N=5,\ 10$, using Theorem \ref{THM4.3}. For the functions
$\varphi_i$, $i=1,2,3$ the graphs of absolute errors $|\varphi_i(z)-\mathring{P}_{\varphi_i}(z)|$, $i=1,2,3$, are given in Fig 2, Fig 3, Fig 4. In these Figures one can see that by increasing value of $N$ absolute errors between optimal interpolation formulas and given functions are decreasing.
 \vspace{-50pt}
\begin{center}
  \includegraphics[width=0.3\textwidth]{m=2_C_0}
  \includegraphics[width=0.3\textwidth]{m=2_C_1}
  \includegraphics[width=0.3\textwidth]{m=2_C_2}\\
  \vspace{-50pt}
  \includegraphics[width=0.3\textwidth]{m=2_C_3}
  \includegraphics[width=0.3\textwidth]{m=2_C_4}
  \includegraphics[width=0.3\textwidth]{m=2_C_5}\\
   \vspace{-50pt}
Fig. 1. {\small Graphs of coefficients of the optimal interpolation formulas (\ref{(1)}) in the case  $N=5$.}
\end{center}
\begin{center}
  \includegraphics[width=0.45\textwidth]{m=2_f1_x2_N=5}
  \includegraphics[width=0.45\textwidth]{m=2_f1_x2_N=10}\\
   \vspace{-60pt}
Fig. 2. { Graphs of absolute errors for  $N=5$ and $N=10$:  $|z^2-P_{z^2}(z)|$.}
\end{center}
 \vspace{-40pt}
 \begin{center}
\includegraphics[width=0.45\textwidth]{m=2_f2_exp2x_N=5}
  \includegraphics[width=0.45\textwidth]{m=2_f2_exp2x_N=10}\\
   \vspace{-60pt}
Fig. 3. {\small Graphs of absolute errors for  $N=5$ and $N=10$:  $|\exp(z)-P_{\exp(z)}(z)|$.}
\end{center}
 \vspace{-30pt}
 \begin{center}
  \includegraphics[width=0.45\textwidth]{m=2_f3_sinx_N=5}
  \includegraphics[,width=0.45\textwidth]{m=2_f3_sinx_N=10}\\
 \vspace{-60pt}
Fig. 4. {\small Graphs of absolute errors for  $N=5$ and $N=10$: $|\sin(z)-P_{\sin(z)}(z)|$.}
\end{center}

The Figure 4 shows exactness our optimal interpolation formula  for the function $\sin(x)$.
\newpage
\begin{center}
\textbf{References}
\end{center}
{\small
\begin{enumerate}
\bibitem{Ahlb67} Ahlberg.J.H, Nilson.E.N, Walsh.J.L, The theory of splines and their applications, Mathematics in Science and Engineering, New York: Academic Press, 1967.
\bibitem{BabHay18}  Babaev.S.S, Hayotov.A.R, Optimal interpolation formulas in $W_2^{(m,m-1)}$ space. arXiv:1802.00562.
\bibitem{deBoor63} de Boor.C, Best approximation properties of spline functions of odd degree, J. Math. Mech. 12, (1963), pp.747-749.
\bibitem{Hay04a} Hayotov.A.R, The discrete analogue of a differential operator and its applications. Lithuanian Mathematical Journal. 2014.-vol54. No2, pp.290-307.
\bibitem{HayShadMil14} Hayotov.A.R , Milovanovic.G.V , Shadimetov.Kh.M,  Optimal quadrature formulas and interpolation splines minimizing the semi-norm in $K_2(P_2)$ space // Milovanovic.G.V and Rassias.M.Th (eds.), Analytic Number Theory, Approximation Theory, and Special Functions, Springer, 2014, -pp.573-611.
\bibitem{Schoen64} Schoenberg.I.J, On trigonometric spline interpolation, J. Math. Mech. 13, (1964), pp.795-825.
\bibitem{Schum81} Schumaker.L, Spline functions: basic theory, Cambridge university press, 2007, 600 p.
\bibitem{Sob61a} Sobolev.S.L, On Interpolation of Functions of $n$ Variables, in: Selected Works of Sobolev.S.L, Springer, 2006, pp. 451-456.
\bibitem{Sob74} Sobolev.S.L, Introduction to the Theory of Cubature Formulas, Nauka, Moscow, 1974, 808 p.
\bibitem{Sob77} Sobolev.S.L, The coefficients of optimal quadrature formulas, in: Selected Works of Sobolev.S.L. Springer, 2006, pp.561-566.
\bibitem{SobVas} Sobolev.S.L, Vaskevich.V.L, The Theory of Cubature Formulas. Kluwer Academic Publishers Group, Dordrecht (1997).
\bibitem{Vas83} Vasilenko.V.A, Spline-fucntions: Theory, Algorithms, Programs, Nauka, Novosibirsk, 1983, 216 p. (in Russian)
\end{enumerate}
}
\medskip

\medskip
{\small
\begin{tabular}{p{9cm}}
Bukhara state university, Bukhara, Uzbekistan\\
\end{tabular}
}

\label{Babaev2}
\end{document}


\section{Introduction}\label{Sec-Int}
The celebrated theorem of Bernstein \cite{Be1915} (as for a simple proof, see \cite{Ni1957}), stating that the only entire minimal graph in the Euclidean $3$-space ${\R}^3$ is a plane,
has played a seminal role in the evolution of surface theory.
For instance, by using this theorem,  we can show the uniqueness theorem for an entire graph of constant mean curvature in ${\R}^3$.
Indeed Heinz \cite{He1955} proved that if $\varphi (u^1, u^2)$ is a $C^2$-differentiable function defined on the open disk $B^2 (R)$ with center at the origin and radius $R (>0)$ in the $(u^{1}, u^{2})$-plane, and the mean curvature $H$ of the graph
$\Gamma_{\varphi}:= \{ (u^{1}, u^{2}, \varphi (u^{1}, u^{2})) \in \R^3 \, ; \,  (u^{1}, u^{2}) \in B^2 (R) \}$ of $\varphi$ satisfies $|H|\geq \alpha >0$, where $\alpha$ is constant, then $R\leq 1/\alpha$. As a corollary, we have a vanishing theorem of mean curvature for entire graphs of constant mean curvature, that is, if  $\Gamma_{\psi}$ is an entire graph of
constant mean curvature, then $H\equiv 0$.
Combining these results, we obtain that the only entire graph of constant mean curvature in ${\R}^3$
is a plane  (see \cite[Section 2.1]{Ke2003}). The Heinz result was extended to the case of graphic hypersurfaces in the Euclidean $(n+1)$-space ${\R}^{n+1}$
by Chern \cite{Ch1965} and Flanders \cite{Fl1966}.

We will study this subject for graphic hypersurfaces in the Lorentz-Minkowski space. Here we recall some basic definitions and fundamental facts. We denote by ${\R}^{n+1}_{1} = (\R^{n+1}, \langle \,  \,  ,  \, \rangle_{L})$ the Lorentz-Minkowski $(n+1)$-space with the Lorentz metric
\[
\langle (x^{1}, \ldots , x^{n}, x^{n+1}), (y^{1}, \ldots , y^{n}, y^{n+1}) \rangle_{L} := x^{1}y^{1} + \cdots + x^{n}y^{n} -x^{n+1}y^{n+1},
\]
where $ (x^{1}, \ldots , x^{n}, x^{n+1})$, $(y^{1}, \ldots , y^{n}, y^{n+1}) \in \R^{n+1}$. Let $M^{n}$ be a smooth orientable $n$-manifold.
An immersion $f\colon M^n \to \R^{n+1}_{1}$ is called {\it space-like}
if the induced metric $g:= f^{\ast} \langle \,  \,  ,  \, \rangle_{L}$ on $M^n$ is Riemannian, and is called {\it time-like} if $g$ on $M^n$ is Lorentzian.
A space-like (resp. time-like) immersion of constant mean curvature is called a {\it space-like (resp. time-like) constant mean curvature immersion}.
In particular, a space-like immersion with vanishing mean curvature is called a {\it space-like maximal immersion} and
a time-like immersion with vanishing mean curvature is called a {\it time-like minimal immersion}.

Let $\psi \colon \Omega \to \mathbf{R}$ be a $C^2$-differentiable function defined on a domain $\Omega$ in $\R^{n}$. We consider the graph $\Gamma_{\psi}:= \{ (u^{1}, \ldots, u^{n}, \psi (u^{1}, \ldots , u^{n}))\in {\R}^{n+1}_{1} \, ; \,
(u^{1}, \ldots , u^{n})\in \Omega \}$ of $\psi$. If $\Omega = \R^{n}$, then $\Gamma_{\psi}$ is called an {\it entire graph} of $\psi$ in ${\R}^{n+1}_{1}$.
If the graph $\Gamma_{\psi}$ of $\psi$ is space-like, the gradient $\nabla \psi := (\psi_{u^1}, \ldots, \psi_{u^n})$ of $\psi$ satisfies $|\nabla \psi|:= \sqrt{(\psi_{u^{1}})^2 +
\cdots + (\psi_{u^{n}})^2} <1$ on $\Omega$, and
\begin{equation}\label{1-1}
\displaystyle nH = \sum_{i=1}^{n} \dfrac{\partial }{\partial u^i}\left( \dfrac{\psi_{u^i}}{\sqrt{1-|\nabla \psi|^2}} \right) = \text{div} \left( \dfrac{\nabla \psi}{\sqrt{1-|\nabla \psi|^2}}\right) ,
\end{equation}
that is,
\[
\displaystyle (1-|\nabla \psi|^{2})\sum_{i=1}^{n} \dfrac{\partial^{2} \psi}{(\partial u^{i})^{2}} +\sum_{i, j=1}^{n}\dfrac{\partial \psi}{\partial u^{i}} \dfrac{\partial \psi}{\partial u^{j}} \dfrac{\partial^{2} \psi}{\partial u^{i} \partial u^{j}}= nH(1-|\nabla \psi|^{2})^{3/2}
\]
holds. Here $H$ is the mean curvature of $\Gamma_{\psi}$ and $\psi_{u^i}:= \partial \psi / \partial u^{i}$.
If the graph $\Gamma_{\psi}$ of $\psi$ is time-like, its gradient $\nabla \psi$ satisfies $|\nabla \psi |>1$ on $\Omega$, and
\begin{equation}\label{1-1-2}
\displaystyle nH = \sum_{i=1}^{n} \dfrac{\partial }{\partial u^i}\left( \dfrac{\psi_{u^i}}{\sqrt{|\nabla \psi|^2 -1}} \right) = \text{div} \left( \dfrac{\nabla \psi}{\sqrt{|\nabla \psi|^2 -1}}\right)
\end{equation}
holds.

For entire space-like maximal graphs in $\R^{n+1}_{1}$, the following uniqueness theorem, called the Calabi-Bernstein theorem, is well-known.
This result was proved by Calabi \cite{Ca1968} for $n\leq 4$ and Cheng-Yau \cite{CY1976} for all $n$.
\begin{fact}[The Calabi-Bernstein theorem]\label{fact1}
Any entire space-like maximal graph is a space-like hyperplane in  $\R^{n+1}_{1}$.
\end{fact}
We remark that an improvement of this theorem and a fluid mechanical interpretation of the duality between minimal graphs in ${\R}^3$ and space-like maximal graphs in $\R^{3}_{1}$ are discussed in \cite{AUY2019} (cf. \cite{Be1958}).

On the other hand, for entire space-like constant mean curvature graphs, the uniqueness theorem does not hold in general. In fact, the graph of the function
\begin{equation}\label{eq-hyperboloid}
\psi (u^1, \ldots, u^n) = \left( (u^1)^2 + \cdots + (u^n)^2 + \dfrac{1}{H^2} \right)^{1/2} \quad (H>0)
\end{equation}
is an entire space-like constant mean curvature $H$ graph which is not a hyperplane. This graph is a hyperboloid.
Many other entire space-like constant mean curvature graphs are constructed by Treibergs \cite{Tr1982}.
Hence the following questions naturally arise. What are the situations where Heinz-type mean curvature estimate and Bernstein-type uniqueness theorem hold for space-like graphs in ${\R}^{n+1}_{1}$?  What are these situations for time-like graphs in ${\R}^{n+1}_{1}$?

In this paper, we perform a systematic study for these questions. In particular, we give a unified description of Heinz-type mean curvature estimates under an assumption on
the gradient bound for space-like graphs and time-like graphs
in ${\R}^{n+1}_{1}$ (Theorem \ref{thm1}). As a corollary, we obtain a unified vanishing theorem of mean curvature for entire graphs of constant mean curvature
with a gradient bound (Corollary \ref{cor1}). In Section \ref{Sec-App}, by applying these results, we answer the questions mentioned above.
Moreover we give a short commentary about the recent study by the authors for entire constant mean curvature graphs in $\R^{n+1}_{1}$.

Finally, the authors gratefully acknowledge the useful comments from Shintaro Akamine and Atsushi Kasue during the preparation of this paper.

\section{Main results}\label{Sec-Main}
The main theorem is stated as follows:
\begin{theorem}\label{thm1}
Let $B^{n}(R)$ be the open $n$-ball with center at the origin and radius $R (> 0)$ in $\mathbf{R}^{n}$,
and $\psi (u^{1}, \ldots , u^{n})$ a $C^2$-differentiable function on  $B^{n}(R)$. Suppose that there exist constants $M>0$ and $k\in \R$ such that
\begin{equation}\label{eq2-0}
\dfrac{|\nabla \psi|}{\sqrt{|1-|\nabla \psi|^2|}}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{k}
\end{equation}
on $B^{n}(R)$. Assume that $\Gamma_{\psi}$ is a space-like or time-like graph of $\psi$ in $\R^{n+1}_{1}$.
If the mean curvature $H$ of $\Gamma_{\psi}$ satisfies the inequality
\[
|H|\geq \alpha >0,
\]
where $\alpha$ is constant, then the following inequality holds:
\begin{equation}\label{eq2-1}
\alpha \leq MR^{2k-1}.
\end{equation}
\end{theorem}

\begin{proof}
Let $\omega$ be the $(n-1)$-form given by
\begin{equation}\label{eq2-a}
\displaystyle \omega = \sum_{i=1}^{n} (-1)^{i-1} \left( \dfrac{\psi_{u^i}}{\sqrt{|1-|\nabla \psi|^2}|} \right)\, du^{1}\wedge \cdots \wedge \widehat{du^i}\wedge \cdots \wedge du^n,
\end{equation}
where $\widehat{du^i}$ means that $du^i$ is omitted from the wedge product. Then we can obtain $d\omega = nH\, du^{1}\wedge \cdots \wedge du^n$
by \eqref{1-1} and \eqref{1-1-2}.
Take any positive number $R'$ satisfying $0< R' <R$. The Stokes theorem yields the following equality:
\begin{equation}\label{eq2-2}
\displaystyle \int_{\overline{B}^{n}(R')} d\omega = \int_{\partial B^{n}(R')} \omega,
\end{equation}
where $\overline{B}^{n}(R')$ is the closure and $\partial B^{n}(R')$ is the boundary of  $B^{n}(R')$.

We may assume $H\geq \alpha >0$ by changing the direction of the normal vector if necessary. Let $V_{n}$ be the volume of the unit closed $n$-ball $\overline{B}^n (1)$ and
$A_{n-1}$ the volume of the unit $(n-1)$-sphere $\partial B^n (1)$. Then, by radial integration, $nV_{n}=A_{n-1}$ holds. The absolute value of the left-hand side of \eqref{eq2-2} becomes
\begin{eqnarray*}
\left| \displaystyle \int_{\overline{B}^{n}(R')} d\omega  \right| & =  & n\int_{\overline{B}^{n}(R')} H\, du^{1}\cdots du^{n} \\
                                                                                & \geq & n\alpha \int_{\overline{B}^{n}(R')} du^{1}\cdots du^{n} = n\alpha (R')^n V_{n}.  \\
\end{eqnarray*}
On the other hand, by using the Cauchy-Schwarz inequality, the absolute value of the right-hand side of \eqref{eq2-2} becomes
\begin{eqnarray*}
\left| \int_{\partial B^{n}(R')} \omega \right| &=& \left| \int_{\partial B^{n}(R')}  \sum_{i=1}^{n} (-1)^{i-1} \left( \dfrac{\psi_{u^i}}{\sqrt{|1-|\nabla \psi|^2|}} \right)\, du^{1}\cdots \widehat{du^i}\cdots du^n \right| \\
                                                                &\leq & \displaystyle \int_{\partial B^{n}(R')} \dfrac{|\nabla \psi|}{\sqrt{|1- |\nabla \psi|^2|}} \left( \sum_{i=1}^{n} \left( du^{1}\cdots \widehat{du^i}\cdots du^n \right)^{2} \right)^{1/2} . \\
\end{eqnarray*}
By \eqref{eq2-0}, we have
\begin{eqnarray*}
\displaystyle \left| \int_{\partial B^{n}(R')} \omega \right| &\leq & M (R')^{2k} \int_{\partial B^{n}(R')} \left( \sum_{i=1}^{n} \left( du^{1}\cdots \widehat{du^i}\cdots du^n \right)^{2} \right)^{1/2} \\
                                                                                &=& M (R')^{n+2k-1} A_{n-1}. \\
\end{eqnarray*}
We thus obtain $n\alpha (R')^{n}V_{n} \leq M(R')^{n+2k-1}A_{n-1}$, that is, $\alpha\leq M(R')^{2k-1}$. The proof is completed by letting $R' \to R$.
\end{proof}

From the above argument, we give the following generalization of the result which was obtained by Salavessa \cite[Theorem 1.5]{Sa2008} for space-like
graphs in ${\R}^{n+1}_{1}$.

\begin{proposition}\label{thm-sala}
Let $D$ be a relatively compact domain (i.e. its closure $\overline{D}$ is compact) of ${\R}^n$ with smooth boundary $\partial D$.
Assume that $\Gamma_{\psi} \subset {\R}^{n+1}_{1}$ is a space-like or time-like graph of  a $C^2$-differentiable function $\psi$.
Set $m_{D}:= \max_{\overline{D}} |\nabla \psi|$ if $\Gamma_{\psi}$ is a space-like graph in $\R^{n+1}_{1}$
and $m_{D}:= \min_{\overline{D}} |\nabla \psi|$ if $\Gamma_{\psi}$ is a time-like graph
in $\R^{n+1}_{1}$. Then, for the mean curvature $H$ of $\Gamma_{\psi}$, we have
\begin{equation}\label{eq-sala-1}
\displaystyle \min_{\overline{D}} |H| \leq \dfrac{1}{n} \dfrac{m_{D}}{\sqrt{|1- (m_{D})^2|}} \dfrac{A (\partial D)}{V (\overline D)}.
\end{equation}
Here $V(\overline{D})$ (resp. $A (\partial D)$\,) is the volume of $\overline{D}$ (resp. $\partial D$).
\end{proposition}

\begin{remark}
To be precise, in \cite[Theorem 1.5]{Sa2008} Salavessa obtained this inequality for space-like graphs of $M\times N$, where $M$
is a Riemannian $n$-manifold and $N$ is an oriented Lorentzian $1$-manifold.
\end{remark}

Indeed, from the Stokes theorem for the $(n-1)$-form $\omega$ given by \eqref{eq2-a}, we obtain
\begin{equation}\label{eq2-2-2}
\displaystyle \int_{\overline{D}} d\omega = \int_{\partial D} \omega.
\end{equation}
Then the absolute value of the left-hand side of \eqref{eq2-2-2} becomes
\[
\left| \displaystyle \int_{\overline{D}}\, d\omega \,  \right| \geq n \, \min_{\overline{D}} |H|\, V (\overline{D}).
\]
Since the function $t/\sqrt{1-t^2}$ is monotone increasing on $0< t< 1$ and $t/\sqrt{t^2 -1}$ is monotone decreasing on $t> 1$, we have
\[
\left| \int_{\partial D} \, \omega \, \right| \leq \dfrac{m_{D}}{\sqrt{|1- (m_{D})^2|}}\, A (\partial D).
\]
We have thus proved the inequality \eqref{eq-sala-1}.

As a corollary of Theorem \ref{thm1}, we give the following unified vanishing theorem of mean curvature for entire space-like graphs and time-like graphs
of constant mean curvature in ${\R}^{n+1}_{1}$.

\begin{corollary}\label{cor1}
Assume that $\Gamma_{\psi}$ is an entire space-like or time-like graph of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ in $\R^{n+1}_{1}$.
If the entire graph $\Gamma_{\psi}$ has constant mean curvature and there exist constants $M>0$ and $\varepsilon >0$ such that
\begin{equation}\label{eq-Be}
\dfrac{|\nabla \psi|}{\sqrt{|1-|\nabla \psi|^2|}}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{(1/2) - \varepsilon}
\end{equation}
on ${\R}^n$, then its mean curvature must vanish everywhere.
\end{corollary}
\begin{proof}
By Theorem \ref{thm1}, the mean curvature $H$ of ${\Gamma}_{\psi}$ satisfies
\[
|H|\leq \dfrac{M}{R^{2\varepsilon}}
\]
on $B^{n}(R)$. We obtain $H\equiv 0$ by letting $R\rightarrow +\infty$.
\end{proof}

\section{Applications}\label{Sec-App}
\subsection{Space-like case}\label{subsec-sp}
Then $|\nabla \psi|< 1$ holds. We give a geometric interpretation for $|\nabla \psi|/ \sqrt{1-|\nabla \psi|^2}$.
When the graph $\Gamma_{\psi}$ is space-like,
\begin{equation}\label{eq2-nu}
\nu (=  \nu (u^1, \cdots, u^n)) =\dfrac{1}{\sqrt{1- |\nabla \psi|^2}} (\psi_{u^1}, \cdots, \psi_{u^n}, 1)
\end{equation}
is the time-like unit normal vector field $\nu$ of $\Gamma_{\psi}$. Since $e_{n+1}=(0, \cdots , 0, 1) \in \R^{n+1}_{1}$ is also time-like,
there exists a unique real-valued function $\theta (= \theta (u^1, \cdots , u^n ) )\geq 0$ such that $\langle \nu , e_{n+1} \rangle_{L} = -\cosh{\theta}$.
This function $\theta$ is called the {\it hyperbolic angle} between $\nu$ and $e_{n+1}$ (see \cite{On1983}). By simple calculation, we have
\begin{equation}\label{eq2-sinh}
\sinh{\theta} = \dfrac{|\nabla \psi|}{\sqrt{1-|\nabla \psi|^2}}.
\end{equation}
We remark that L\'opez \cite{Lo2020} studies the function $|\nabla \psi|/ \sqrt{1-|\nabla \psi|^2}$ from the viewpoint of the
Dirichlet problem.

From Theorem \ref{thm1}, we obtain the following Heinz-type mean curvature estimate.

\begin{corollary}\label{corA}
Let $\Gamma_{\psi}$ be a space-like graph of  a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ on $B^{n}(R)$ in ${\R}^{n+1}_{1}$ and
$\theta$ the hyperbolic angle between $\nu$ and $e_{n+1}$. Suppose that there exist $M >0$ and $k\in \R$ such that
\begin{equation}\label{eq-corA}
\sinh{\theta} \leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{k}
\end{equation}
on $B^{n}(R)$. If the mean curvature $H$ of $\Gamma_{\psi}$ satisfies the inequality $|H|\geq \alpha >0$, where $\alpha$ is constant, then $\alpha \leq MR^{2k-1}$ holds.
\end{corollary}

By virtue of Fact \ref{fact1} and Corollary \ref{cor1}, we obtain the following Bernstein-type theorem for entire space-like constant mean curvature
graphs in ${\R}^{n+1}_{1}$.

\begin{corollary}\label{corB}
If an entire space-like graph $\Gamma_{\psi}$ of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ in ${\R}^{n+1}_{1}$ has constant mean curvature and there exist constants $M>0$ and $\varepsilon >0$
such that
\begin{equation}\label{eq-corB}
\sinh{\theta}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{(1/2) - \varepsilon}
\end{equation}
on $\R^n$, then it must be a space-like hyperplane. Here $\theta$ is the hyperbolic angle between $\nu$ and $e_{n+1}$.
\end{corollary}

Corollary \ref{corB} is optimal because there exists an example which is not congruent to a space-like hyperplane and satisfies \eqref{eq-corB} for $\varepsilon =0$. In fact,
the function $\psi$ given by \eqref{eq-hyperboloid} is not linear and its graph has constant mean curvature. Moreover it satisfies
\[
\sinh{\theta} = H \sqrt{(u^1)^2+\cdots +(u^n)^2},
\]
which is the equality in \eqref{eq-corB} for $\varepsilon =0$.

From the case where $\varepsilon =1/2$ in Corollary \ref{corB}, we obtain the following  Liouville-type result.

\begin{corollary}\label{cor2}
Let $\Gamma_{\psi}$ be an entire space-like constant mean curvature graph of  a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$
in $\R^{n+1}_{1}$ and $\theta$ the hyperbolic angle between $\nu$ and $e_{n+1}$.
If $\sinh{\theta}$ is bounded on ${\R}^n$, then ${\Gamma}_{\psi}$ must be a space-like hyperplane.
\end{corollary}

The following result is equivalent to Corollary \ref{cor2}.

\begin{corollary}{\cite[Corollary 9.5]{AMR2016}}\label{cor3}
If an entire space-like graph $\Gamma_{\psi}$ of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ in $\R^{n+1}_{1}$ has constant mean curvature and there exists some constant $C(<1)$ such that
$|\nabla \psi|\leq C$ on ${\R}^n$, then it must be a space-like hyperplane.
\end{corollary}

Indeed, if $|\nabla \psi | \leq C $ holds on $\R^n$, then we obtain
\[
\dfrac{|\nabla \psi|}{\sqrt{1-|\nabla \psi|^2}}\leq \dfrac{C}{\sqrt{1-C^2}}
\]
on $\R^n$, that is, $\sinh{\theta}$ is bounded because the function $t/\sqrt{1-t^2}$ is monotone increasing on $0<t<1$.
By the same argument, the converse also holds, that is, if $\sinh{\theta}$ is bounded on ${\R}^n$, then there exists some constant $C (<1)$ such that
$|\nabla \psi|\leq C$ on ${\R}^n$.

As described in \cite[Section 9.3]{AMR2016},  Corollary \ref{cor3} is a case of nonparametric form of the following  value-distribution-theoretical property
of the Gauss map.
\begin{fact}[\cite{Ai1992, Xi1991}, \cite{Pa1990} for a first weaker version]\label{fact2}
Let $f\colon M^n \to {\R}^{n+1}_{1}$ be a complete space-like constant mean curvature immersion.
If the image of the time-like unit normal vector field (i.e., Gauss map) $\nu (M^{n})$ is contained in a geodesic ball in the hyperbolic $n$-space ${\H}^{n}(-1)$,
then $f(M^{n})$ is a space-like hyperplane in ${\R}^{n+1}_{1}$.
\end{fact}
In fact, every entire space-like constant mean curvature graph in ${\R}^{n+1}_{1}$ is complete (see \cite{CY1976}, note that the converse statement is also true, that is,
every complete space-like hypersurface in ${\R}^{n+1}_{1}$ is an entire graph (\cite{ARS1995}))
and the image $\nu ({\R}^n)$ is contained in a geodesic ball in ${\H}^{n}(-1)$ if and only if $\sqrt{1- |\nabla \psi|^2}$ is bounded away from zero on ${\R}^n$. However the proof of Fact \ref{fact2} is not as simple as our argument in this paper because it is obtained via an application of the Omori-Yau maximum principle (\cite{CX1992}, \cite{Om1967}, \cite{Ya1975}). See the book \cite{AMR2016} for details.

For a space-like graph in ${\R}^{n+1}_{1}$, $\sinh{\theta}$ can be replaced with $1/\sqrt{1-|\nabla \psi|^2}$ in \eqref{eq-corA} and \eqref{eq-corB}.
In fact,
\[
\sinh{\theta}= \dfrac{|\nabla \psi|}{\sqrt{1-|\nabla \psi|^2}} < \dfrac{1}{\sqrt{1-|\nabla \psi|^2}}
\]
holds. Thus we have the following result.

\begin{corollary}\label{cor-Do}
If an entire space-like graph $\Gamma_{\psi}$ of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ in ${\R}^{n+1}_{1}$ has constant mean curvature and there exist constants $M>0$ and $\varepsilon >0$
such that
\[
\dfrac{1}{\sqrt{1-|\nabla \psi|^2}}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{(1/2) - \varepsilon}
\]
on $\R^n$, then it must be a space-like hyperplane.
\end{corollary}

From Corollary \ref{cor-Do}, we can show the following proposition which was given by Dong.

\begin{proposition}[\cite{Do2007}]
Let $\Gamma_{\psi}$ be an entire space-like constant mean curvature graph of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$ in ${\R}^{n+1}_{1}$.
If the function $\psi$ satisfies
\[
\dfrac{1}{\sqrt{1-|\nabla \psi|^2}} = o (r)  \quad {as} \quad r\to +\infty,
\]
where $r=\sqrt{(u^1)^2+\cdots +(u^n)^2}$, then $\Gamma_{\psi}$ must be a space-like hyperplane.
\end{proposition}

\begin{remark}
In \cite{Do2007}, Dong showed a generalization of this result for an entire space-like graph of parallel  mean curvature in the pseudo-Euclidean $(m+n)$-space ${\R}^{m+n}_{n}$ of index
$n$.
\end{remark}

\subsection{Time-like case}\label{Subsec-time}
Then $|\nabla \psi|> 1$ holds. By Theorem \ref{thm1}, we obtain the following Heinz-type estimate for the mean curvature of
time-like graphs in $\R^{n+1}_{1}$.

\begin{corollary}\label{thm2}
Let $\Gamma_{\psi}$ be a time-like graph of  a $C^2$-differentiable function $\psi (u^{1}, \ldots , u^{n})$ on $B^{n}(R)$ in $\R^{n+1}_{1}$.
Suppose that there exist constants $M>0$ and $k\in \R$ such that
\begin{equation}\label{eq3-0}
\dfrac{|\nabla \psi|}{\sqrt{|\nabla \psi|^2 -1}}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{k}
\end{equation}
on $B^{n}(R)$. If the mean curvature $H$ of $\Gamma_{\psi}$ satisfies the inequality
\[
|H|\geq \alpha >0,
\]
where $\alpha$ is constant, then the following inequality holds:
\begin{equation}\label{eq3-1}
\alpha \leq MR^{2k-1}.
\end{equation}
\end{corollary}

Moreover, we obtain the following result from Corollary \ref{cor1}.

\begin{corollary}\label{cor4}
If an entire time-like graph $\Gamma_{\psi}$ of a $C^2$-differentiable function $\psi (u^{1}, \ldots, u^{n})$
in $\R^{n+1}_{1}$ has constant mean curvature $H$ and
there exist constants $M>0$ and $\varepsilon >0$ such that
\begin{equation}\label{eq-Be2}
\dfrac{|\nabla \psi|}{\sqrt{|\nabla \psi|^2 -1}}\leq M \biggl{(}(u^1)^2+\cdots +(u^n)^2  \biggr{)}^{(1/2) - \varepsilon}
\end{equation}
on ${\R}^n$, then it must be minimal (i.e. $H\equiv 0$).
\end{corollary}

We note that there exist many entire time-like minimal graphs in $\R^{n+1}_{1}$. In fact, for a $C^2$-differentiable function $h(t)$ of one variable which satisfies $h'(t) >0$ for
all $t\in \R$, $\psi (u^1 , \cdots, u^n):= u^n +h(u^1)$ gives an entire time-like minimal graph ${\Gamma}_{\psi}$ in ${\R}^{n+1}_{1}$.

We do not know whether Corollary \ref{cor4} is optimal or not. However, for example, from the case where $\varepsilon =1/2$ and $n=2$ in Corollary \ref{cor4},
we provide the following result.

\begin{corollary}\label{cor5}
Let $\Gamma_{\psi}$ be an entire time-like constant mean curvature graph of a $C^2$-differentiable function $\psi (u^1, u^2)$ in $\R^{3}_{1}$.
If $|\nabla \psi|/\sqrt{|\nabla \psi|^2 -1}$ is bounded on $\R^2$, then $\Gamma_{\psi}$ is a translation surface,
that is, $\psi (u^1, u^2) = h (u^1) + k(u^2) $, where $h (u^1)$ is a $C^2$-differentiable function of $u^1$ and $k (u^2)$ is a $C^2$-differentiable function of $u^2$.
\end{corollary}

The latter half of the statement of Corollary \ref{cor5} follows from the following fact.

\begin{fact}[\cite{Ma1989, Mi1987}]\label{fact3}
Every entire time-like minimal graph in ${\R}^{3}_{1}$ is a translation surface.
\end{fact}

\subsection{Closing Remark}\label{subsec-CR}
This paper is focused on the study of entire space-like graphs and time-like graphs in ${\R}^{n+1}_{1}$.
Here we explain the recent development of the study of entire constant mean curvature graphs in $\R^{n+1}_{1}$.
The following definitions are given in \cite{AHUY2019}.
Let $\psi \colon \Omega \to \R$ be a $C^2$-differentiable function defined on a domain $\Omega$ in ${\R}^n$.
A point where $|\nabla \psi | <1$ (resp. $|\nabla \psi|> 1$, $|\nabla \psi| = 1$) is called a {\it space-like} (resp. {\it time-like}, {\it light-like}) point of
${\Gamma}_{\psi}$. If $\Omega$ contains only space-like (resp. time-like) points, then $\Gamma_{\psi}$ is space-like (resp. time-like).
We remark that the first fundamental form $g$ of the graph ${\Gamma}_{\psi}$ in ${\R}^{n+1}_{1}$ is degenerate at each light-like point.
If $\Omega$ contains both space-like points and time-like points, then ${\Gamma}_{\psi}$ is called of {\it mixed type}.
Akamine, the first author, Umehara and Yamada proved the following fact.

\begin{fact}{\cite[Corollary C]{AHUY2019}}\label{fact-light}
An entire real analytic constant mean curvature graph in ${\R}^{n+1}_{1}$ which has no time-like points and does have a light-like point must be a light-like hyperplane.
\end{fact}

See \cite{AHUY2019} (also \cite{UY2019}) for the precise definition of constant mean curvature graph which does not depend on its causal type.
By Fact \ref{fact-light}, we know that there exists no entire real analytic constant mean curvature graph in ${\R}^{n+1}_{1}$ which has both space-like points and light-like points although there exist many nontrivial examples of entire space-like constant mean curvature graphs.
Moreover, the first and third authors, Kokubu, Umehara and Yamada \cite[Remark 2.2]{HKKUY2017} showed that there does not exist any real analytic nonzero constant mean curvature graph of mixed type in $\R^{n+1}_{1}$.
\begin{remark}
Several entire zero mean curvature graphs of mixed type in $\R^{3}_{1}$ were found in recent years. A systematic construction of this class is given in \cite{FKKRUYY2016}.
See \cite{HST2020} for details about mixed type surfaces in Lorentzian 3-manifolds.
\end{remark}


\section{Introduction}
An estimated 17.5 million people died from cardiovascular diseases in 2005, accounting 30\% of deaths around the world~\cite{who}. Over the last fifteen years, magnetic resonance imaging (MRI)  has become a routine modality for the determination of patient cardiac morphology. An accurate  extraction of the anatomical information is crucial to obtain reproducible quantitative measurements to support the diagnosis and follow-up of cardiac pathologies as well as to develop new clinical applications. Manual delineation of the anatomical structures is a long and tedious task, which is prone to inter- and intraobserver variations. Therefore, it is highly desirable to develop techniques for automatic segmentation.

For years, left ventricle (LV) segmentation in cardiac MRI has been the main focus of several research groups. For instance, it has been the subject of two different segmentation challenges\footnote{MICCAI 2009 - Cardiac MRI Left Ventricle Segmentation Challenge and MICCAI 2011 - The STACOM'11 Cardiac Left Ventricular Segmentation Challenge}. We refer the reader to \cite{suri,petitjean2011} for detailed reviews on the topic. Due to its larger anatomical complexity, right ventricle (RV) segmentation has not been as extensively explored. However, recent work~\cite{caudron} has highlighted the importance of RV structure and function assessment in the management of cardiac diagnosis and the determinant role of a high quality RV segmentation.  The aim of this paper is to participate in the RV Segmentation Challenge in Cardiac MRI, an initiative to gather up and review state of the art methods in the field.

\begin{figure}
\begin{center}
\begin{tabular}{ccc}
\includegraphics[scale=0.15]{imgs/basal_ed.png} & \includegraphics[scale=0.15]{imgs/basal_ed2.png} &
\includegraphics[scale=0.15]{imgs/basal_es.png} \\
\includegraphics[scale=0.15]{imgs/apical_ed.png} &
\includegraphics[scale=0.15]{imgs/apical_ed2.png} &
\includegraphics[scale=0.15]{imgs/apical_es.png}
\end{tabular}
\end{center}
\caption{The shape of the ventricles can vary significantly between basal (top row) and apical slices (bottom row), end-diastole (left column) and systole (center column) and different patients (left and center column vs. right column). The use of atlases exhibiting a great variability can have a negative effect on the label fusion process.  }
\label{fig:examples}
\end{figure}

In the remaining sections of the paper, we first describe the proposed segmentation approach, followed by the results obtained with the segmentation challenge data. Finally, we discuss on the obtained results and some conclusions are presented.

\section{Methodology}
The RV segmentation method presented here is inspired on a previous approach for 3D whole heart segmentation using an atlas-based propagation segmentation propagation framework~\cite{zhuang} using a single atlas, combining an intensity image and a label image. Following this idea, we make use of multiple atlases for segmentation propagation in order to extract the endocardium and epicardium of the RV.

\begin{figure}
\begin{center}
\includegraphics[width=\columnwidth]{imgs/diagram.png}
\end{center}
\caption{Automatic RV coarse-to-fine segmentation framework.  }
\label{fig:framework}
\end{figure}

Multi-atlas segmentation propagation is typically performed  in two stages. First,  the unseen image is registered to all the intensity images of the atlases, followed by the mapping of the label images to the unseen image space using the obtained transformations. Then, at the second stage, the transformed labeled images are combined into a single segmentation of the unseen image through a label fusion method.

As our segmentation is performed on 2D and cardiac images can show large variability (Fig.~\ref{fig:examples}), it is necessary to perform an atlas selection that chooses the best suited atlases for a particular unseen image slice. For this matter, we use a multi-label fusion strategy based on the STAPLE algorithm
that incorporates a fast locally normalised cross correlation (LNCC) based ranking combined with
a consensus based region-of-interest selection and an iterative Markov Random Field into the a multi-label STAPLE algorithm~\cite{warfield}. For further details on this multi-label fusion strategy, here denoted as STEPS, we refer the reader to~\cite{cardoso}.

Instead of applying only once the multi-atlas propagation framework to an unseen image, we use a coarse-to-fine strategy that incrementally refines the segmentation by means of the multi-atlas propagation. The coarser segmentation obtained at each propagation level is used as a mask to start the registration gradually up to the finest level. The different steps of our coarse-to-fine strategy, illustrated in Fig.~\ref{fig:framework}, are described in the following.

\paragraph{Phase I. Heart localization.} At a first stage, we seek to define a region of interest (ROI) that encloses the heart in the unseen image. At this phase, the unseen image is globally registered to the atlases using a block matching approach~\cite{ourselin}. The obtained transformations  are applied to the atlas' label images, which are all fused using  majority voting. Figure~\ref{fig:framework} shows the typical output ROI for the whole heart.

\paragraph{Phase II. Rough segmentation.} The unseen image ROI obtained from phase I is used as a mask in the registration step of phase II. The use of a mask allows the removal of structures surrounding the heart, e.g. the liver and the ribs, that can affect the registration procedure.

At this phase, two different registration steps are perfomed. The atlas are rigidly registered~\cite{ourselin} to the masked unseen image, followed by a non-rigid alignment  using a fast free form deformation algorithm~\cite{modat}. After non-rigidly transforming the atlases, i.e. both intensity and label images, to the unseen image space, we apply an NCC-based ranking strategy~\cite{cardoso} in order to determine for each unseen slice  which are the most suitable registered  atlases slices to use. We selected only the 10\% best ranked atlases as several of the registrations fail to provide a good alignement of the unseen image and the atlases. The selected atlases slices were then fused using the STEPS algorithm~\cite{cardoso} in order to obtain an initial RV rough segmentation (Fig.~\ref{fig:framework}). As at this stage only the top 10\% atlases are used, the results from thie fusion step still don't have a high quality.

\paragraph{Phase III. Segmentation refinement.} In order to accurately initialize an affine  registration for all the atlases, all the label images are affinely aligned~\cite{ourselin} to the estimated rough segmentation. The transformed label images are fused through majority voting. As in phase II, this mask is used to remove surrounding structures in the final image-to-image non-rigid registration. With the non-rigid transformations, the label images are propagated and resampled using a nearest neighbor interpolation to preserve the binary nature of the segmentations. Finally, we make use of the multi-label fusion algorithm~\cite{cardoso} that first ranks every atlas slice based on the NCC and then fuses locally the labels of the N best ranked atlases $\left(N=\frac{\left|{atlas}\right|}{3}\right)$.

\begin{table}
\begin{center}
\caption{Summary of the registration and fusion methods involved at each phase.}
\label{tab:coarsetofinesummary}
\begin{tabular}{llllll}
\hline\noalign{\smallskip}
\multicolumn{1}{c}{Phase} & \quad & \multicolumn{1}{c}{Registration} & \quad & \quad & \multicolumn{1}{c}{Fusion} \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
I & & Image-to-image global rigid registration & & \quad & Majority Voting  \\
& & & &\\
II & & 1) Image-to-image affine registration  & & \quad & STEPS  \\
   & & 2) Image-to-image non-rigid registration & & \quad & \\
   & & & &\\
III & & 1) Mask-to-mask affine registration & & \quad & 1) Majority voting \\
    & & 2) Image-to-image non-rigid registration & & \quad & 2) STEPS \\
\hline
\end{tabular}
\end{center}
\end{table}

Table~\ref{tab:coarsetofinesummary} summarizes the different registration and fusion algorithms involved at each phase of the coarse-to-fine strategy.

\section{Experiments and Results}
The RV Segmentation Challenge in Cardiac MRI provided the participants with an initial set of 16 MR images containing manual annotations of the RV epicardium and endocardium 2D contours. We used these annotations to build up label images. These, in combination with the intensity MR images, made up our atlas set. As each dataset contained end-diastole (ED) and end-systole slices, separate atlases were built for ED and ES.
\begin{table}
\begin{center}
\caption{Image-based Dice metric (DM) and Hausforff distance (HD) of RV segmentation for the test data set.}
\label{tab:results}
\begin{tabular}{lllll}
\hline\noalign{\smallskip}
 & \quad  & \multicolumn{1}{c}{DM}  & \quad &\multicolumn{1}{c}{HD}  \\
 & \quad & \multicolumn{1}{c}{mean(std)} & \quad & \multicolumn{1}{c}{mean(std)}  \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
\textbf{ED} & & & \\
Endocardium & & 0.83 (0.17) & &9.77 (7.88)  \\
Epicardium & &0.86 (0.13) & &10.23 (7.22) \\
\textbf{ES} & & & &\\
Endocardium & & 0.72 (0.27)& & 11.41 (10.49)  \\
Epicardium & & 0.77 (0.23) & &11.81 (9.46) \\
\hline
\end{tabular}
\end{center}
\end{table}

Training and parameter optimization of the algorithm was performed using a leave-one-out strategy using the training datasets. On a subsequent stage, the challenge organizers provided 16 additional cases (test set) to which our method was applied. In such scenario, the complete training set was used as an atlas. In average, the segmentation of a complete case (i.e. segmentation of both ED and ES) took 12 minutes. Evaluation of the results was performed in terms of both image- and patient-based criteria. Image-based criteria included the Dice metric (DM) and the Hausdorff distace (HD), whereas patient-based criteria, assessed endocardial volume at ED and ES, ejection fraction (EF) and ventricular mass (VM). Further details on the evaluation metrics can be found on the challenge website\footnote{http://www.litislab.eu/rvsc/}.

Table~\ref{tab:results} summarizes the results of our method for RV segmentation in terms of image-based metrics. The results show that the algorithm performs better at ED than ES. This can be explained by the fact that the image quality is higher at ED than at ES. By analyzing the DM and HD on a slice-basis, it can be seen that the proposed method has a high performance on basal slices, with an average $DM=0.93$ and $HD=4.31$. However, the segmentations on apical slices are of lower quality affecting the overall DM and HD scores.  Figure~\ref{fig:segmentations} shows segmentation results in six different cases obtained from the challenge that illustrate the differences between basal and apical segmentations.

\begin{figure}
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=0.2\columnwidth]{imgs/P21_7_ED.png} & \includegraphics[width=0.2\columnwidth]{imgs/P22_6_ED.png} &
\includegraphics[width=0.204\columnwidth]{imgs/P01_1_ED.png} \\
\includegraphics[width=0.2\columnwidth]{imgs/P21_6_ES.png} &
\includegraphics[width=0.201\columnwidth]{imgs/P22_6_ES.png} &
\includegraphics[width=0.201\columnwidth]{imgs/P01_0_ES.png}
\end{tabular}
\end{center}
\caption{Example segmentations on cases 21 (left), 22 (center) and 01 (right) at both ED (top row) and ES (bottom row). Basal slices achieve high quality segmentations, whereas it fails  on apical slices, specially, in the epicardium extraction.}
\label{fig:segmentations}
\end{figure}

Figure~\ref{fig:patientresults} shows the regression plots for the endocardial volume at ED and ES.  For both cases, the regression coeffictient is good (0.96 and 0.97 respectively) and the spread of the values is relatively low, demonstrating that the influence of the poor segmentation in the apical slices doesn't have much influence in the endocardial volume estimation. Errors obtained in the computation of the ejection fraction (EF) and ventricular mass (VM), which are based on volume estimations, are reported in table~\ref{tab:resultspatient}.

\section{Discussion and Conclusions}
In this paper, we have presented a multi-atlas propagation and segmentation approach for the segmentation of the RV epicardium and endocardium in ED and ES. As the RV segmentation is done currently in 2D and RV contours exhibit large variation along the heart, we make use of a multi-label ranking criterion, based on the local normalized cross correlation, to select the best atlases for label fusion.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.45\columnwidth]{imgs/EDvolume} &
\includegraphics[width=0.45\columnwidth]{imgs/ESvolume}
\end{tabular}
\end{center}
\caption{Regression curve for the endocardial volume at ED (left) and ES (right).  }
\label{fig:patientresults}
\end{figure}
\begin{table}
\begin{center}
\caption{Ejection fraction (EF) and ventricular mass (VM) errors.}
\label{tab:resultspatient}
\begin{tabular}{lllll}
\hline\noalign{\smallskip}
Error & \quad & EF & \quad & VM \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
Mean & & 0.15 & & 0.37  \\
Standard deviation & & 0.09 & & 0.33 \\
\hline
\end{tabular}
\end{center}
\end{table}

The results show that the method perform well on average but there are some cases in which it fails. In particular, our method has a very high performance on the basal slices, whereas the scores obtained for apical slices are lower (Figure~\ref{fig:segmentations}). These can be explained by two factors: 1) The image quality at the apical slices is rather low. As the registration, the atlas ranking and the label fusion are intensity-based, poor image quality can affect the results, and 2) the regions to be segmented at the apical slices are rather small, which implies that the atlases' masks are also small at these slices. When the masks are used in the registration process to supress undesired structures, the remaining information is insufficient, causing the intensity-based registration, rigid or non-rigid, to fail.

Nonetheless, these problems can be solved if the registrations are performed in 3D, where more image information can be used for both registration and fusion. Due mainly to slice thickness of this type of images, our current registrations are done in 2D. We will focus our efforts in the extension of the framework so that it can be applied to 3D.

In conclusion, we have presented a fully automatic segmentation method applied to the extraction of RV epicardium and endocardium. The results show that the method is in average succesful in the segmentation results. The segmentations obtained in the basal slices are of very high quality, whereas the results in the apical ones are less good mainly because the registration stage fails to properly align the atlases with the unseen image.

\section*{Acknowledgements}
M.A. Zuluaga and S. Ourselin have been founded by the EPRSC Grant 'Grand Challenges: Translating biomedical modelling into the heart of the clinic'   (EP/ H02025X/1). M.J. Cardoso was supported by a scholarship from the Fundacao para a Ciencia e a Tecnologia, Portugal (
SFRH/BD/43894/ 2008). M.J. Cardoso and S. Ourselin receive funding from
the EPSRC grant ' Intelligent Imaging: Motion, Form and Function Across Scale' (EP/H046410/1) and the
CBRC (Ref. 168).