Wael Saideni
email: wael.saideni@xlim.fr

Fabien Courreges
email: fabien.courreges@unilim.fr

David Helbert
email: david.helbert@univ-poitiers.fr

Jean Pierre Cances
email: cances@ensil.unilim.fr

Robust Spatiotemporal Convolutional Long Short-Term Memory Algorithm for Video Prediction

Keywords: Video prediction, deep learning, neural networks, computer vision, ConvLSTM, memory flow, hidden states

The use of recurrent neural networks in several applications has allowed to capture impressive results, especially in various applications such as video prediction and it has become a promising direction of scientific research. In this paper, we introduce a novel algorithm for video prediction called "Robust Spatiotemporal Convolutional Long Short-Term Memory (Robust-ST-ConvLSTM) Algorithm" that outperforms the stateof-the-art approaches. Robust-ST-ConvLSTM is a memory flow algorithm based on higher order ConvLSTM. This memory flow algorithm is holding the spatiotemporal information to optimize and control the prediction abilities of the ConvLSTM cell. Our approach is developed in the specific context of predicting future frames based on historical observations, and we experimentally validate the ability of the proposed algorithm on two spatiotemporal datasets, including a moving variant of MNIST dataset of handwritten digits, and KTH which is a human motion dataset.

I. INTRODUCTION

Video prediction or predicting what happens in the next frames is the key component of intelligent decision making systems. It is also, an emerging field of computer vision and deep learning that is facing many challenges. Actually, these predictive systems have many real-world applications such as video surveillance or human and buildings security which is one of the most frequently debated issues nowadays.

Video prediction networks are based on historical information gathered from continuous and unlabeled video frames. These networks aim to forecast future frames in a video after having some previous images. Formally, we suppose X t ∈ R w×h×c is the t-th frame of a dynamic scene X = (X t-n , ..., X t ) with n frames, where w, h, and c denote width, height and number of channels, respectively. The main target from this project is to predict the next m frames Y = (Y t+1 , ..., Y t+m ) from the input X.

Since deep learning has demonstrated its performances and effectiveness in video processing and because feature extraction from images using traditional supervised machine learning methods is challenging, deep learning for video prediction could be a promising research direction and a very powerful tool for better performances. However, despite the huge advancement is deep learning models, video prediction remains a big challenge in terms of image deblurring and longterm prediction. Therefore, we propose in this paper a new algorithm called Robust Spatiotemporal Convolutional Long Short-Term Memory (Robust-ST-ConvLSTM) algorithm for video prediction as a long-term prediction algorithm that can forecast more than 20 frames. Also, it outperforms the stateof-the-art architectures. Thus, we empirically demonstrate the performance of our model on the main two challenges mentioned above. In fact, the main cell of our architecture is a modified ConvLSTM cell. Actually, ConvLSTM architecture is commonly used for spatiotemporal predictive systems with a traditional roadway for the memory state. However, this approach is outdated and not optimal for many reasons and that's why some changes on its main structure have been done. Firstly, the predictive algorithms based on ConvLSTM focus on the long-term stochastic properties of the data rather than its spatial distortion. Nevertheless, for video prediction, spatial and temporal data structure are important and should be taken into account. Secondly, the spatiotemporal cell state designed will handle low-level and semantic aspects of the spatiotemporal video data which are significantly important to generate future frames. Thirdly, the memory flow will propose new cell state and hidden state transition functions different from the ordinary ones which have not fully used the homogeneity of the input and the output space. Finally, the ConvLSTM model proposes an explicit temporal information encoding in each cell [START_REF] Xingjian | Convolutional LSTM network: A machine learning approach for precipitation nowcasting[END_REF]. This first-order Markovian architecture updates the hidden states with information from the previous time step only which harden the capturing of longrange temporal correlations. In addition, most first-order RNN, i.e the current hidden state is updated based on one previous hidden state, suffer from the gradient vanishing problem in back-propagation resulting in a difficulty in learning RNN to model long-term dependency in data [START_REF] Soltani | Higher order recurrent neural networks[END_REF].

The remaining paper is organized as follows: Section II discusses related works. In Section III, we describe the main idea behind our proposed algorithm and its key components. In Section IV, we evaluate the capability of Robust-ST-ConvLSTM for multi-step video prediction on two spatiotemporal datasets, including a synthetic dataset of handwritten digits and a human motion dataset and report its performance by comparing it against the state-of-the-art algorithms. Finally, Section V provides conclusion and the future research directions.

II. RELATED WORKS

A. Optical flow based methods

Many research projects have proposed video prediction solutions based on optical flow or dense trajectory [START_REF] Walker | Dense optical flow prediction from a static image[END_REF]- [START_REF] Sedaghat | Hybrid Learning of Optical Flow and Next Frame Prediction to Boost Optical Flow in the Wild[END_REF]. In fact, optical flow is applied to report motion information about objects of successive frames. Technically, these approaches take the given dynamic scene as input to forecast the optical flow of the future frame. The obtained result is then merged with the last input frame to generate the future predicted video frame. However, those approaches that necessitate supervised training, use training datasets that contain optical flow information which is not obviously provided in the commonly used video datasets.

B. Deep Learning based methods

While the optical flow based models use the motion information to predict the frames, neural network approaches analyze the frames and extract their features in order to exploit the spatiotemporal representation to forecast the next frames. In this section, recent deep learning models for video prediction will be discussed after being classified into three categories: recurrent neural networks, convolutional networks and generative networks. Although these neural networks based methods are better than the traditional optical-flow-based solutions in terms of performances, they are challenging and produce sometimes blurry results. Obviously, it is a promising research area.

1) Recurrent models

Recurrent networks are commonly used for video sequences related problems since they are considered as sequential data with spatio-temporal representation. Recurrent neural networks (RNN) have demonstrated considerable success in video prediction research works that are detailed in [START_REF] Terwilliger | Recurrent flow-guided semantic forecasting[END_REF]- [START_REF] Vora | Future segmentation using 3d structure[END_REF] . In fact, along with the advancements in neural networks architectures, video prediction has been studied extensively in recent years. Zhang et al. proposed a ConvLSTM-based architecture where hidden states are updated along a z-order curve [START_REF] Zhang | Z-order recurrent neural networks for video prediction[END_REF]. The model presents a novel training approach based on two Z-Order Recurrent Networks (Znet): Znet-Predictor and Znet-Probe. Since the majority of video prediction algorithms based on ConvLSTM have duplicated features with same functionality in both cell state and hidden state of the LSTM unit, Znet came up with a novel routes updating to enhance the hidden states. Technically, to trick the neural network, the model is set to choose inputs that minimize the loss function instead of updating weights and biases that minimize the cost. W. Lotter et al. [START_REF] Lotter | Deep Predictive Coding Networks for Video Prediction and Unsupervised Learning[END_REF] presented a predictive neural network (PredNet) architecture. This network aims to forecast future video frames in dynamic scenes. Technically, every layer in the network makes local predictions and only sends the deviations from those predictions to the following layers.

The PredNet model is a series of recurrent blocks that make local predictions, then subtracted from the input before being forwarded to the subsequent network layer. C. Lu et al. [START_REF] Lu | Flexible Spatio-Temporal Networks for Video Prediction[END_REF] propose a Flexible Spatio-Temporal Network (FSTN). This model enables the generation of the frames lying between the observed frames in order to output slow-motion video sequences. Also, it proposes a novel loss function to optimize the training phase of the model. The architecture described above is based on two main models: extrapolation model and interpolation model. Both of them are considered as spatio-temporal autoencoders. However, the extrapolation model has a guided training phase by the the ground truth frames feeding each layer by the supervised information needed, while the interpolation model does not need the ground truth images. Another difference of the two models lies in their definition. The interpolation is the estimation of a value between given data points but the extrapolation is useful when looking for a value that is either higher or lower than the values in the dataset. A recent RNN architecture was proposed by Wang et al. in [START_REF] Wang | PredRNN: Recurrent neural networks for predictive learning using spatiotemporal lstms[END_REF]. The idea behind this research work remain behind the new spatiotemporal LSTM (ST-LSTM) unit that take out and memorize spatial appearances and temporal variations simultaneously since for video prediction we need to consider both the spatial and the temporal structures. In fact, the Predictive Recurrent Neural Network (PredRNN) is based on spatiotemporal memory flow which allows the memory cells to move vertically across stacked RNN layers and horizontally through all RNN states. This approach is different from stacked LSTM. Actually, in stacked LSTM, memory states are updated independently from the visual features which means that the first layer of the present time step could ignore the information memorized by the last layer at the previous time step. However, in PredRNN, a memory cell is introduced to handle the information between different time steps. Another problem is solved in [START_REF] Wang | PredRNN: A Recurrent Neural Network for Spatiotemporal Predictive Learning[END_REF]. The new memory cell can handle long-term and short-term information at the same time which can limit the predictive performances of the model. So, a pair of memory cells is used and and explicitly decoupled in order to satisfy the different variations. This model reduce the loss of visual information from the very first layer to the top of the recurrent network. Furthermore, another learning strategy was proposed called reverse scheduled sampling. This strategy enables to learn temporal dynamics from longer periods of the input video and reduce the training discrepancy between the encoding network and the prediction network.

2) Convolutional models

Different from recurrent neural networks, convolutional networks are feed-forward neural networks that are commonly used for computer vision challenges such as visual prediction. Many models are based on convolutions for video prediction. One of these architectures is a multi-model combining temporal and spatial sub-networks which is proposed in [START_REF] Yan | Mixpred: Video Prediction Beyond Optical Flow[END_REF] and called MixPred. The future frame prediction approach described is divided into two parts: a temporal model for modeling the time series of the input video and a spatial sub-network to model the spatial texture on the content. Then, the authors tested an information fusion method for feature map interaction between the two parts. This approach allows to copy the unchanging pixels from the last frame thanks to the temporal mask which means that the predicted frame has the same clearness as the original frame. Also, synchronously exchanging temporal and spatial information enables to fill the changed pixels in order to have a complete predicted image. This model uses only convolutional layers but it could be theoretically enhanced by using another models like the generative networks. However, the architecture of this approach, in particular the choice of the hyperparameters, differs with the model which their performances are compared to. The model described above could be used not only in future frames prediction but also in several applications such as object tracking, action recognition and video compression. In [START_REF] Liu | Video frame synthesis using deep voxel Flow[END_REF], the model trains a deep neural network to generate video frames by flowing pixel values from existing ones instead of initializing them from scratch. The model, called Deep Voxel Flow (DVF), takes usually 3 frames from the video scene without pre-processing: two frames are taken as input and the third frame is used as the generated target. This approach is based on the idea of borrowing voxels (3D-pixels) from the adjacent frames to generate more realistic results. The architecture is composed of a convolutional encoderdecoder to forecast the voxel flow and a volume sampling layer to generate the target image. As in [START_REF] Yan | Mixpred: Video Prediction Beyond Optical Flow[END_REF], the model can predict the in-between frames (interpolation) and the future frames (extrapolation) of the input dynamic scene. The voxel flow, used to sample the input frames with the volume sampling function to synthesize the target frame, has two main components: the spatial component and the temporal one. The spatial element is the optical flow for the predicted frame and the temporal part is used to form a color in that frame. The framework described above aims to predict one frame but it can naturally be extended to a multi-frame prediction framework with a fairly simple manipulation. In fact, the target becomes a 3D volume and not 2D image and the learning rate will be reduced to maintain stability in the training phase. In addition, the spatiotemporal coherence is maintained because of the preserve of local correlations due to the convolutions across the temporal layers. The strength of this model is that it combines the advantages of the optical-flow-based approach and the newer neuralnetwork-based models. Also, it can be trained and tested on any real-world video with any resolution. However, it fails in scenes with repetitive patterns. Also, it generates some blurry scenes, like most of neural-network-based implementations.

3) Generative models

Generative models are used to generate new samples from the same distribution as the input data. The target behind training generative models is to learn a probability distribution that is similar to the data's probability distribution. In video prediction, the models described above aim to output a single eventual outcome. However, generative approaches generate a wide spectrum of feasible predictions. The most common network structure in the field of video prediction and image generation in general is Generative Adversarial Networks (GAN). These networks are composed of two sub-networks jointly trained, the discriminator and the generator, to create fake samples that look like real data. Technically, the generator fools the discriminator by generating new samples from a random noise (e.g. Gaussian noise). Then the discriminator features the probability distribution function that describes real data. Nevertheless, in video prediction, some conditions could be added to the general implementation of GAN in order to forecast the future frames. In [START_REF] Kwon | Predicting future frames using retrospective cycle gan[END_REF], a generative approach was proposed to prediction frames based on cycle GAN. The main model is composed of one generator and two discriminators. In fact, the generator uses the retrospective cycle to predict both future and past frames and we train it with reversed input sequences. Moreover, one discriminator is dedicated to identify fake frames while the other is implemented to distinguish the sequences that contain fake frames which is crucial in forecasting temporally consistent frames. Technically, the loss function and the network architecture make this approach special when we compare it with the general formulation of GAN networks. Since this model enables to predict a limited number of frames before generating blurry images, a multi-frame prediction strategy is employed. The model starts by forecasting the next frame from an input video. Then, it constructs a new input video by concatenating the last frames of the input video and the predicted frame. Finally, the new input video will enable the prediction of the next frame. This strategy is repeated until we get the desired number of predicted frames. In [START_REF] Oprea | A review on deep learning techniques for video prediction[END_REF], the authors insisted on the fact that conditional Generative Adversarial Networks (cGAN) are suitable for video frames prediction because it can guarantee the spatiotemporal coherence between the predicted frames and the input video. Another approach is discussed in [START_REF] Tulyakov | MoCoGAN: Decomposing motion and content for video generation[END_REF] and is based on the idea of dividing the video signal into two parts: content and motion. Content to specify the objects in the sequence and motion to describe their movements. The model is based on mapping a sequence of random vectors to a sequence of frames in order to generate the predicted videos. These random vectors are composed of two parts: one for the content and the other for the motion. Since this framework is based on GAN, discriminators are used to learn motion and content decomposition in an unsupervised way by introducing a new adversarial learning scheme.

III. OVERVIEW OF THE PROPOSED ROBUST SPATIOTEMPORAL CONVLSTM ALGORITHM A. From LSTM to ConvLSTM

The idea behind the proposed algorithm is based on Convolutional LSTM (ConvLSTM) which is Long Short Term Memory (LSTM) network applied on high dimensional data.

1) LSTM Long short Term Memory Network is considered as an advanced type of RNN that was designed and developed by Hochhreiter and Schmidhhuber (1997) [START_REF] Hochreiter | Long Short-Term Memory[END_REF] to solve the vanishing gradient problem of standard RNNs. Theoretically, RNNs are designed to learn long term dependencies. However, in practice, many issues appear such as vanishing gradient that prevents those neural network to learn long term dependencies. Therefore, it has been proven that LSTM is a powerful tool to remember information for longer period of time. Indeed,the main idea behind LSTM consists of connecting the previous information to the future task. The main structure of LSTM based neural networks is the same: it consists of a chain of LSTM modules. However, the structure of those modules depends on the application. When we look at the LSTM cell in Figure 1, we notice that the information coming from the cell state c t passes through the forget gate that decide which information is going to be forgotten thanks to the sigmoid layer that outputs a number between 0 and 1. Then, the input gate uses the input x t and the hidden state h t-1 to update the cell state. Then, a tanh layer outputs new candidate values ĉt that have the possibility to be added to the cell state. The cell state update is created from the combination of ĉt and i t . Now, everything is ready to update the cell state. Firstly, the old cell state c t-1 is multiplied by f t then i t * ĉt is added. Finally, the output of the LSTM unit will be based on the cell state c t , the input x t and the hidden state h t-1 . Indeed, a sigmoid gate is applied to decide the parts of the cell state that will be involved in the output process. Then, the cell state c t is put through tanh and then multiplied by the output of the sigmoid layer. The main target of this last step is to output the new hidden state h t . To sum up the mechanism of LSTM: This neural network unit has 3 inputs: the input x t , the cell state c t-1 and the hidden state h t-1 that will be passed through 3 different gates in order to output 2 structures: the cell state c t and the new hidden state h t . The mechanism described above is explained by the following equations:

i t = σ(w i × x t + s i × h t-1 ) f t = σ(w f × x t + s f × h t-1 ) o t = σ(w o × x t + s o × h t-1 ) ĉt = tanh(w ĉ × x t + s ĉ × h t-1 ) c t = f t • c t-1 + i t • ĉt h t = o t • tanh( ĉt ) (1)
Where σ is the sigmoid function, × is a simple multiplication and • denotes the Hadamard product.

2) ConvLSTM

Although LSTM is considered as a powerful network for dealing with temporal relationship, its main drawback is that it is unable to handle spatial information because we need to flatten high dimensional data to 1D vectors to be compatible to the input common structure. However, Spatiotemporal data are commonly used is many applications such as video surveillance. So, we were forced to look for a new structure where we take advantage of LSTM by integrating spatiotemporal data. Convolutional LSTM is used to capture the spatial dimension for the prediction mode. The special feature of ConvLSTM is that the inputs x t , the cell states c t , the hidden states h t and the 3 gates are 3D tensors. In addition, the convolution operation is used instead of simple matrix multiplication as shown in the following equations:

I t = σ(W i * X t + S i * H t-1 ) F t = σ(W f * X t + S f * H t-1 ) O t = σ(W o * X t + S o * H t-1 ) Ĉt = tanh(W ĉ * X t + S ĉ * H t-1 ) C t = F t • C t-1 + I t • Ĉt H t = O t • tanh( Ĉt ) (2)
Where * denotes the convolution operation and • denotes the Hadamard product.

B. Robust Spatiotemporal ConvLSTM proposed algorithm

The Proposed Robust Spatiotemporal ConvLSTM (Robust-ST-ConvLSTM) algorithm is a memory flow algorithm based on higher order ConvLSTM. To make it simple, the novel algorithm aim to decide the cell state C t not only from the previous hidden state H t-1 but also from N previous hidden states (H t-2 , ..., H t-N ) (N will be fixed by the user and it can only affect the computational time). The second part of the algorithm is to implement a memory flow to hold spatiotemporal information to optimize and control the prediction capacities of ConvLSTM. In fact, the memory flow will be a second cell state for spatiotemporal data. However, the cell state will not be removed and will handle temporal Indeed, the novel algorithm uses a stack of ConvLSTM units to learn the spatial correlations and the temporal dynamics from the input video. These feature will be used later to forecast the future frames. So, a novel transition function is introduced based on spatiotempral memory flow and is able to leverage a deterministic number of previous hidden states. In the original implementation of ConvLSTM, the temporal memory states C t are updated only from one time step to another. However, in video prediction, the consecutive frames are having a close data distributions in the spatial dimensions and also many temporal correlations. That's why, we need to exploit these properties to make better predictions in terms of quality performances. Therefore, we believe that this higher order ConvLSTM based on memory flow will take advantage from the global motion changes of the consecutive frames and the information of the spatiotemporal memory to predict future frames. It is obvious that the memory state is updated horizontally in the original stacked ConvLSTM model. However, the previous model is enhanced by updating the memory state horizontally (the cell state) and also vertically (the spatiotemporal memory state) as illustrated in Figure 3. This approach ameliorates the way we handle the spatiotemporal information from the input to the output and enables to connect all the recurrent units of the entire network. The main equations of the new robust spatiotemporal unit represented in Figure 4 are:

I t = σ(W i * X t + f (H l t-1 , ..., H l t-N )) F t = σ(W f * X t + f (H l t-1 , ..., H l t-N )) Ĉt = tanh(W ĉ * X t + f (H l t-1 , ..., H l t-N )) C l t = F t • C l t-1 + I t • Ĉt I t = σ(W i * X t + M i * ST M l-1 t ) F t = σ(W f * X t + M f * ST M l-1 t ) Ĉ t = tanh(W ĉ * X t + M ĉ * ST M l-1 t ) ST M l t = F t • ST M l-1 t + I t • Ĉ t O t = σ(W ox * X t + f (H l t-1 , ..., H l t-N ) + W oc * C l t + W ostm * ST M l t ) H l t = O t * tanh(W 1×1 * [C l t , ST M l t ]) (3)
Where σ is the sigmoid activation function, * and • denote the convolution operator and the Hadamard product respectively. Like the original ConvLSTM, I t and I t e: the input gates, F t and F t : the forget gates, Ĉt and Ĉ t : the potential candidates for the cell states,O t : the output gate. X t denotes the input at the time step t. H l t denotes the hidden state of the lth layer at the time step t. C l t is the memory state of the lth layer at the time step t. ST M l t denotes the spatiotemporal memory of the lth layer at the time step t. f is the function that should be designed to combine N previous hidden states. The design of the f is quite difficult since it must satisfy the following conditions: the spatial structure of the hidden states must be preserved, the size of the filters that control the previous hidden states must increase with the time steps in order to capture the context of these structures and finally, the algorithm's complexity must not explode. To implement f we tested two main approaches. The first approach aims to return the mean value of all elements in the input tensor that handle the previous hidden states. In Robust-ST-ConvLSTM, the feedback signal is generated by combining multiple preceding hidden states. Therefore, the state of the Norder Robust-ST-ConvLSTM is recursively updated with the following f function:

f (H l t-1 , ..., H l t-N ) = 1 N N n=1 W hn H l t-n (4)
Analogous to the filter structures used in signal processing, the second approach in designing the f function in inspired from recursive least squares filters [START_REF] Cances | Joint channel estimation and data demodulation algorithms for fast time varying band limited frequency selective Rayleigh fading channels: A comparison study[END_REF]. It is now based on the weighted sum of the previous hidden states. Consequently, f is straightforward:

f (H l t-1 , ..., H l t-N ) = 1 N N n=1 α n n W hn H l t-n ( 5
)
Where α is the forgetting factor. The parameter α (0 < α < 1) gives more weight to recent hidden states. The gates of the Robust Spatiotemporal unit are no longer dependent on the the hidden state and the temporal memory state from the previous time step of the same layer. However, they depend on the previous hidden states from previous time steps at the same layer and the spatiotemporal memory state.

To be clear, the first layer in a stacked ConvLSTM model at time step t receives the spatiotempral memory of the last layer in the stacked model of the previous time step as illustrated in Figure 3 (ST M 1 t = ST M L t-1 with L is the number of stacked layers). So, we adopt the original structure of ConvLSTM and we added a second gated structure for the spatiotemporal memory ST M l t . However, the final hidden state H l t depends on the combination of the temporal memory state C l t and the spatiotemporal memory state ST M l t . The spatiotemporal memory parameter is dedicated to reduce the loss of spatiotempral information in the video sequences from the first layer to the last layer of the network. Besides, the previous hidden states used as input for the ConvLSTM blocks are implemented to expand the visibility of the neural units about the context of the current events at different time steps. It is clear that the proposed model increases the number of parameters when we compare it with the standard ConvLSTM but it will prevent as from unnecessarily expanding the Con-vLSTM model to obtain the same performances.

IV. PERFORMANCE EVALUATION, COMPARISON AND DISCUSSION

A. Datasets

As far as we are concerned, there are currently no datasets for video prediction because it is an emerging area of research. However, researchers basically use motion video datasets such as KTH and MovingMNIST used to compare the performances of our proposed algorithm with the state-of-the-art approaches.

1) KTH This dataset [START_REF] Schuldt | Recognizing human actions: a local svm approach[END_REF] has 2391 video sequences of 6 human actions (Walking , Jogging, Running, Boxing, Hand waving, Hand clapping) performed by 25 people in 4 different scenarios (outdoors s1, outdoors with scale variation s2, outdoors with different clothes s3 and indoors s4 as illustrated in Figure 5). Static cameras were used to capture the video scenes with 25fps as a frame rate. The sequences have a length of 4 seconds in average with a frame size of 160 × 120. The videos are stored in 600 video files for each combination of 25 subjects, 6 actions and 4 scenarios.

2) Moving MNIST The commonly used Moving MNIST dataset [START_REF] Srivastava | Unsupervised Learning of Video Representations using LSTMs[END_REF] contains 10,000 sequences of moving handwritten digits with random velocities. Each video of the dataset is a 20-frame dynamic scene and its frame size is 64 × 64.

B. Compared methods and performance metrics 1) Compared methods

To evaluate the performance of our proposed Robust-ST-ConvLSTM, we compare it with the performance of some advanced video prediction models:

• ConvLSTM: is commonly used for spatiotemporal predictive systems with a traditional roadway for the memory state. This algorithm is mentioned in almost every research work as the least efficient approach. However, it is the source of inspiration for video prediction algorithms based on recurrent neural networks. • PredRNN 2017: based on the spatiotemporal LSTM (ST-LSTM) unit that take out and memorize spatial appearances and temporal variations simultaneously. • PredRNN 2021: In this algorithm, a pair of memory cells is used and explicitly decoupled in order to enhance the performances of the previous algorithm and surpass its limitations. In addition, another learning strategy was proposed called reverse scheduled sampling.

2) Performance metrics

Because the results are video frames, we will use the most commonly used metrics to evaluate the quality of images between the ground truth and the prediction. Those metrics are Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM) [START_REF] Horé | Image Quality Metrics: PSNR vs[END_REF]. The PSNR measures, in decibels, the quality ratio between the original frame and the predicted one. The higher the PSNR, the better the quality of the predicted image. The PSNR is calculated by:

P SN R(Y, Ŷ ) = 10log 10 max 2 Ŷ 1 N N i=1 (Y -Ŷ ) 2 (6)
Where Y is the ground truth, Ŷ is the generated prediction, N is the number of pixels and max Ŷ is the maximum value of the frame intensities. The SSIM measures the similarity between two images in terms of luminance, contrast and structure. It is calculated as follows:

SSIM (Y, Ŷ ) = (2µ Y µ Ŷ + C 1 ) + (2σ Y Ŷ + C 2 ) (µ 2 Y + µ 2 Ŷ + C 1 )(σ 2 Y + σ 2 Ŷ + C 2 ) (7)
Where µ Y and µ Ŷ are the average of Y and Ŷ , respectively, σ Y and σ Ŷ are the variance of Y and Ŷ , respectively, σ Y Ŷ is the covariance of Y and Ŷ . C 1 and C 2 are constants. The higher the SSIM, the greater similarity between two images.

C. Implementation details

The proposed algorithm is implemented with Python 3.6 and Pytorch 1.4.0 as a deep learning framework. Pytorch is used because it offers an effective way to manipulate tensors or multi-dimensional matrices needed to store and process multidimensional data. We use Adam optimizer to train our model which is an optimization algorithm that combine the properties of AdaGrad and RMSProp algorithms to provide an optimization algorithm that is faster than the commonly used Stochastic Gradient Descent (SGD) algorithm especially with sparse data. A minibatch of 2 sequences is chosen at each training iteration and it is reduced to the maximum to handle the out of memory problem of our GPU. We choose a learning rate of 0.0001

D. Simulation results

1) on KTH dataset

Table I presents quantitative results of the proposed algorithm and state-of-the-art networks and the corresponding frame-wise comparisons are shown in Figure 6 and Figure 7. We adopt PSNR and SSIM as evaluation metrics. We can obviously confirm that our proposed algorithm show significant improvements in terms of short-term and long-term forecasting over the commonly used ConvLSTM approach. In fact, it increases the average PSNR and SSIM over the same number of predicted frames by 26% and 21.31%, respectively, by comparing it with the algorithm mentioned above. Also, it performs favorably against the PredRNN-v2017 and the PredRNN-v2021 algorithms of Wang et al. Our Robust-ST-ConvLSTM (with α = 0.9) performs better than PredRNN-v2021 by 1.72% and 2.77% in terms of PSNR and SSIM, respectively. These empirical results demonstrate the effectiveness and the efficiency of the Robust Spatiotemporal Convolutional LongShort-Term Memory algorithm in predicting future frames. In accordance with these results, Figure 8 that compares representative generated frames, proves that our algorithm outperforms the state-of-the-art approaches in terms of future movement and frames details. Robust ST-ConvLSTM predicts more accurate motion trajectories into the future because of the memory flow component that strengthen the long term prediction ability of the ConvLSTM cell and also because of updating the ConvLSTM cell using information from some previous time steps. we can notice also that the second approach in designing f which is inspired from recursive least squares filters slightly outperforms the first approach in terms of PSNR and SSIM. This means that further research work could be done in order to determine the optimal value of α that gives the best PSNR and SSIM performances. In this work, various values of α have been tested randomly (0 < α < 1) and the optimum one among them was the selected value 0.9. The presented results and the computational cost depend on the number of memory units used for feedback. In our implementation, we used 3 hidden states which means that we have 3r-order Robust-ST-ConvLSTM. Furthermore, the number of hidden states can affect the performances of our model in terms of the quality of its output and also in terms of the computational process. From the previous observations about the value of α and the number of hidden states, we can confirm that a trade-off should be done between quality performances and computational costs, in future research work, to have the best performances without training a computationally very expensive algorithm.

2) on Moving MNIST Table II presents the performance of the evaluated models on the Moving MNIST dataset by predicting the next 10 frames from the previous 10 input frames. We use the similarity index measure (SSIM) and the Peak signal-to-noise ratio (PSNR) for evaluation. As shown from table II and Figures 9 and 10, our architecture performs well against the state-of-the-Fig. 6. Frame-wise PSNR comparisons of different models on KTH dataset after 100 000 iterations Fig. 7. Frame-wise SSIM comparisons of different models on KTH dataset after 100 000 iterations art approaches in both metrics. Figure 9 and Figure 10 show the frame-wise PSNR and SSIM comparisons of different approaches on MNIST dataset. The results of these figures prove the ability of the Robust-ST-ConvLSTM in predicting future frames. Also, they prove that our approach outperforms the previous models on all the predicted frames. The memory flow algorithm based on 3rd order ConvLSTM with α = 0.9 increases the average PSNR over the 10 predicted frames by 3.15% by comparing it with PredRNN (Wang et al., 2021). However, it outperforms the same approach by 0.22% in terms of SSIM. This means that this metric could not be a good evaluation metric in this case. Moreover, Our approach performs favorably against the traditional ConvLSTM approach in terms of PSNR and SSIM. It brings 14.59% PSNR improvement and 26.95% SSIM improvement over ConvLSTM based frames prediction approach. These numerical results are confirmed by Figure 11 that shows the quality of the 10 predicted frames generated by the different approaches. Robust-ST-ConvLSTM outputs clearer frames. However, the state-of-theart algorithms produce blurry images. This means that Robust-ST-ConvLSTM is more precise and sure about the future variations which proves its robustness against the other longterm prediction algorithms mentioned above. we can notice also that the recursive least squares filters based approach in designing f has approximately similar results as the first approach and that for different values of α. Different from KTH dataset, the value of the parameter α does not affect the quality performances of the outputs but it affects the computational cost of our algorithm since a number of multiplications are added to the calculation process. This means that, for MNIST dataset, only the first approach of designing f , which is based on returning the mean value of the previous hidden states, is taken into consideration. we propose a new recurrent neural network (Robust-ST-ConvLSTM) for video prediction. It is based on new robust spatiotemporal unit inspired from the well-known ConvL-STM structure. This spatiotemporal unit rely on two different approaches in order to strengthen its prediction abilities: a memory flow to handle the spatiotemporal information and a higher order ConvLSTM approach that enable the cell states to decide their values from previous hidden states. Our approach outperforms the state-of-the-art research works on different datasets, including KTH dataset for human motion and Moving MNIST.

In conclusion, video prediction is a promising research direction and can be used in different applications such as video surveillance, video compression and intelligent decisionmaking systems. While great work has been done in video prediction, there is still a place for improvement especially with the continuous advancement in deep learning techniques.

Fig. 1 .

 1 Fig. 1. The structure of a standard LSTM module

Fig. 2 .

 2 Fig. 2. The structure of convolutional LSTM

Fig. 3 .

 3 Fig. 3. The main structure of Robust Spatiotempral LSTM

Fig. 8 .

 8 Fig. 8. Prediction examples on the KTH data set, where we predict 20 frames into the future based on the past 10 frames

Fig. 9 .

 9 Fig. 9. Frame-wise PSNR comparisons of different models on MNIST dataset after 100 000 iterations

Fig. 10 .

 10 Fig. 10. Frame-wise SSIM comparisons of different models on MNIST dataset after 100 000 iterations

TABLE I QUANTITATIVE

 I EVALUATION OF DIFFERENT ALGORITHMS ON THE KTH DATASET. THE METRICS ARE AVERAGED OVER THE 20 PREDICTED FRAMES. HIGHER SCORES INDICATE BETTER PREDICTION RESULTS

	Model	PSNR(dB)	SSIM
	ConvLSTM (Shi et al., 2015)	23.009	0.704
	PredRNN (Wang et al., 2017)	27.624	0.839
	PredRNN (Wang et al., 2021)	28.502	0.831
	Robust-ST-ConvLSTM	28.828	0.848
	Robust-ST-ConvLSTM with α = 0.9	28.992	0.854

TABLE II QUANTITATIVE

 II EVALUATION OF DIFFERENT ALGORITHMS ON THE MNIST DATASET. THE METRICS ARE AVERAGED OVER THE 10 PREDICTED FRAMES. HIGHER SCORES INDICATE BETTER PREDICTION RESULTS

	Model	PSNR(dB)	SSIM
	ConvLSTM (Shi et al., 2015)	28.380	0.705
	PredRNN (Wang et al., 2017)	30.569	0.869
	PredRNN (Wang et al., 2021)	31.525	0.893
	Robust-ST-ConvLSTM	32.490	0.894
	Robust-ST-ConvLSTM with α = 0.9	32.520	0.895

ACKNOWLEDGMENTS

This work was supported in part by the sensors generation project of Nouvelle Aquitaine region (2018-1R50214).


Alexandre Fenneteau

David Helbert
email: david.helbert@univ-poitiers.fr

Pascal Bourdon

Imane M'rabet

Christine Fernandez-Maloigne

Rémy Guillevin

A size-adaptative segmentation method for better detection of multiple sclerosis lesions

Keywords: Artificial intelligence, multiple sclerosis, CNN, small lesions, unbalanced segmentation, detection

Purpose: The automatic segmentation of multiple sclerosis lesions on magnetic resonance images is an open research task aiming to bring more reproducibility in the radiological visual assessment of the disease while reducing the burden of this time-consuming task. The development of artificial intelligence has led to significant improvements in computer-aided diagnosis tools for radiology. It exists several approaches for the segmentation of multiple sclerosis lesions using convolutional neural networks. However, the small lesions are frequently neglected by those algorithms despite their importance. We propose here an adaptable method to improve the detection of small lesions. Approach: The problem of small lesions detection mainly comes from the under-representation of those lesions at a voxel level and the segmentation loss function. The presented method consists in weighting the lesion importance during the training of a convolutional neural network depending on lesion size to correct the impact of voxel lesion imbalances.

Results: With our method, the lesion segmentation computed with the Dice score is only slightly improved but the detection sensitivity is significantly improved at the cost of a limited augmentation of lesion false positive rate. The F1 score has been substantially improved with the correct set of parameters. The improved prediction quality of segmentation maps has been confirmed visually with the help of a radiologist. Conclusions: The described method improves the lesion detection by giving more importance to small lesions during the multiple sclerosis lesions segmentation learning, bringing a better help for radiologists towards a better impact for the patient care.

Introduction

Multiple Sclerosis (MS) is a chronic and autoimmune disease affecting the central nervous system that causes disability and cognitive impairment. The disease has been estimated to touch more than 2,000,000 people worldwide in 2016 according to Wallin et al. [START_REF] Wallin | Global, regional, and national burden of multiple sclerosis 1990-2016: a systematic analysis for the global burden of disease study 2016[END_REF] and is considered as the commonest non-traumatic disabling disease to affect young adults. [START_REF] Dobson | Multiple sclerosis a review[END_REF] The disease causes inflammatory demyelination lesions in the brain and spinal cord that lead to axonal degeneration and eventually neuron death as described by Baecher et al. [START_REF] Baecher-Allan | Multiple sclerosis: mechanisms and immunotherapy[END_REF] 1 The diagnosis of the disease is based on the revised McDonald criteria [START_REF] Thompson | Diagnosis of multiple sclerosis: 2017 revisions of the mcdonald criteria[END_REF] in which the Magnetic Resonance (MR) exam plays an important role by allowing the radiologist to assess the dissemination of the lesions in space and time. Furthermore, the MR exam is particularly used for the follow-up of the disease to spot possible new and enlarging lesions. Combined with the clinical state evaluated by the EDSS score, the MRI follow-up assess response to treatment, allowing therapeutic adjustment if necessary. [START_REF] Brisset | New ofsep recommendations for mri assessment of multiple sclerosis patients: special consideration for gadolinium deposition and frequent acquisitions[END_REF] The MR T2 sequence is set as the reference sequence to detect MS lesions as hyper-intensities, [START_REF] Thompson | Diagnosis of multiple sclerosis: 2017 revisions of the mcdonald criteria[END_REF] but T2-fluid-attenuated inversion recovery (FLAIR) images are more used in practice since it suppresses the hyper-intensities caused by free water while keeping lesion-related hypersignal.

MS lesions screening is performed by radiologists by comparing visually successive MR images and is considered as a daily repetitive and time-consuming task. Moreover, this implies an obvious part of subjectivity, essentially depending on the level of expertise of the radiologist. Detection has been associated with high inter-and intraobserver variability. [START_REF] Cabezas | Improved automatic detection of new t2 lesions in multiple sclerosis using deformation fields[END_REF] A performant automatic segmentation tool for MS lesions would reduce valuable analysis time, bring more reproducibility and give more metrics previously inaccessible.

A large number of methods for automatic segmentation of MS lesions have been proposed in the last two decades. These methods are based on statistical models, atlases, machine learning models and more recently by making use of deep learning as described by Danelakis et al. [START_REF] Danelakis | Survey of automated multiple sclerosis lesion segmentation techniques on magnetic resonance imaging[END_REF] Deep learning has been extensively used for this segmentation problem for its performances with more than a hundred articles referenced by Shoeibi et al. [START_REF] Shoeibi | Applications of deep learning techniques for automated multiple sclerosis detection using magnetic resonance imaging: A review[END_REF] Other methods, considered less efficient, are gradually being abandoned.

It has been pointed out by Kaur et al. [START_REF] Kaur | State-of-the-art segmentation techniques and future directions for multiple sclerosis brain lesions[END_REF] that one of the challenges for MS lesion segmentation is the variability of lesion size. In particular small lesions are generally less detected or less segmented as shown by Vang et al., [START_REF] Vang | Synergynet: A fusion framework for multiple sclerosis brain mri segmentation with local refinement[END_REF] Coronado et al. [START_REF] Coronado | Deep learning segmentation of gadoliniumenhancing lesions in multiple sclerosis[END_REF] and Valverde et al. [START_REF] Valverde | One-shot domain adaptation in multiple sclerosis lesion segmentation using convolutional neural networks[END_REF] This means that small lesions are less segmented and detected than bigger ones because of lesion heterogeneity in representation, size and appearance. However, the problem of small lesion detection is seldom mentioned or addressed despite its impact on the lesion detection and the clinical value of those lesions. Among the top-performing methods, on the ISBI MS segmentation challenge [START_REF] Carass | Longitudinal multiple sclerosis lesion segmentation: resource and challenge[END_REF] leaderboard [START_REF] Wallin | Global, regional, and national burden of multiple sclerosis 1990-2016: a systematic analysis for the global burden of disease study 2016[END_REF] the lesion true positive rate (LTPR) or detection sensitivity is rarely above 55% for best methods with a low lesion false positive rate (LFPR) often lower than 20%. It means that, in general, models learn to favor detecting correct lesions than detecting a lot of them.

However, for the radiologist, it is more important to spot a maximum of MS lesions for computer-aided detection (CAD) tools and especially the smallest ones. As a matter of fact, big lesions are easy to see compared to small lesions that can be missed, especially when the lesion load is high. In addition, for the detection of new lesions on consecutive MR exam, new lesions are often small and difficult to spot despite their value for the diagnostic and for the treatment adjustment. It is also simpler to discard a falsely detected lesion than to spot a not detected one when using CAD.

A method for a better lesion detection sensitivity and particularly for small lesions is of interest for clinical use. To the best of our knowledge, the lesion size variability problem is taking into account during the training of a deep learning model for MS lesion segmentation only by Zhang et al. [START_REF] Zhang | All-net: Anatomical information lesion-wise loss function integrated into neural network for multiple sclerosis lesion segmentation[END_REF] In their method, they include a lesion-wise module in which all lesions are considered as spheres with fixed volume centered on the lesion to segment to allow the model to learn lesion detection without lesion size influence.

First we proposed a simple method for the better detection of small lesions in MS lesion seg-mentation using a Convolutional Neural Network (CNN). The method consists in weighting the loss function for learning MS lesion segmentation in order to prioritize the detection of those lesions. Then, we evaluate the weighting method on an in-house dataset and show how it can improve the lesion detection for a better radiological CAD. Finally, we evaluated our method on the ISBI-2015 dataset as an external validation of the method.

Materials and methods

The segmentation task consists in assigning for every pixel or voxel in 3D a class: lesion or not lesion in our case. Whereas the detection suggests that the object of interest is spatially detected without a voxel-wise constraint, and is most of the time associated with bounding boxes. We can see segmentation as a particular case of detection. Therefore, the segmentation quality has been used for white matter and MS lesion to reflect the detection though it is a wrong oversimplification of the problem as pointed out by Carass et al. [START_REF] Carass | Evaluating white matter lesion segmentations with refined sørensen-dice analysis[END_REF] A good segmentation score can be associated with a low detection although the detection is, in fact, the most radiologically valued. As a matter of fact, it is more important to spot any potential lesions than to be absolutely correct on their real volumes and edges, though both properties are of interest. It is then important to define radiologically relevant metrics to evaluate, design and train automatic algorithms for MS lesion segmentation and detection.

Metrics

For white matter and MS lesion segmentation, the Dice score metric has become a golden-standard for validating the segmentation quality. [START_REF] Carass | Evaluating white matter lesion segmentations with refined sørensen-dice analysis[END_REF] However, it is a segmentation measurement at a voxel point of view that does not necessarily reflect the lesion detection. The lesion detection is generally evaluated with the F1 score that relies on the same principles but is computed at the lesion scale instead of the voxel scale.

Both metrics come from the Sørensen formula:

SD = 2|X ∩ Y | |X| + |Y | , (1)
where X and Y are two sets and | • | is the cardinality operator. In the binary context, where only two classes exist, the Sørensen-Dice can be written:

SD = 2T P 2T P + F P + F N = 2 P P V × T P R P P V + T P R , (2)
where T P is the number of true positives, F P the number of false positives, F N the number of false negatives, T P R = T P T P +F N is the true positive rate or sensitivity and P P V = T P T P +F P is the positive predictive value or precision.

In this paper we adopt the definition of Dice score and F1 score as Carass et al. [START_REF] Carass | Longitudinal multiple sclerosis lesion segmentation: resource and challenge[END_REF] The segmentation Dice score is computed at voxel level whereas the F1 score is computed at a lesion level.

By considering voxel-wise previously described metrics with V prefix and L prefix the lesion-wise metrics, the Dice and F1 score can be written:

Dice = 2V T P 2V T P + V F P + V F N = 2 V P P V × V T P R V P P V + V T P R , (3)
F 1 = 2LT P 2LT P + LF P + LF N = 2 LP P V × LT P R LP P V + LT P R , (4)
where LT P counts lesion overlap between lesion masks, LF P counts predicted lesions that are not in the ground truth mask and LF N counts lesions in the ground truth mask that are not in the predicted mask as formalized by Carass et al. [START_REF] Carass | Longitudinal multiple sclerosis lesion segmentation: resource and challenge[END_REF] Both metrics are the harmonic mean of precision and sensitivity but computed at a voxel scale for the Dice score and at a lesion scale for the F1 score. They measure then the similarity of two sampled as a compromise between detecting a maximum of elements in the ground truth (sensitivity) and predicting a minimum of elements that are not in the ground truth (precision). A common metric used in lesion detection is the lesion false positive rate LF P R = LF P LF P +LT P = 1 -LP P V as described by Zhang et al. [START_REF] Zhang | All-net: Anatomical information lesion-wise loss function integrated into neural network for multiple sclerosis lesion segmentation[END_REF]

Lesion size weighting

When applied to MS lesion segmentation, both Dice and F1 score are expected to be maximized to have a good lesion and voxel lesion detection. However, the implementation of a differentiable loss function based on F1 score is in practice really complex and that explains why the Dice loss is generally used. [START_REF] Kamraoui | MSSEG-2 challenge proceedings: Multiple sclerosis new lesions segmentation challenge using a data management and processing infrastructure[END_REF] The problem of Dice loss is that it is computed at a voxel level, meaning that the lesion scale is not taken into account at all in the learning. When applied on MS lesion segmentation with different lesion sizes, this leads the model to learn with more voxel from big lesions than from voxel from small lesions since their volume is negligible compared to big ones. Since big lesions are easier to spot and the most prevalent, it conducts the trained model to avoid segmenting small lesions that do not account much in the loss for a big batch though false positive are still penalized.

It is then observed in most methods that the Dice score is high, the LF P R is close to 0 and the LT P R close to 0.5, meaning that the model segments almost every time big lesions but does not take the risk to detect more difficult one (the smallest).

The objective of the presented method is to formalize a differentiable loss function for segmen-tation that gets closer to the lesion detection F1 score and especially increases the lesion sensitivity or LT P R for small lesions. To meet those expectation, the method aims at increasing the importance of small lesions, which are neglected by most voxelwise segmentation trained models due to the over-representation of large lesion voxels and to the loss functions.

To give more importance to small lesions, we designed a lesion voxel-weighting method depending on lesion size. The visible lesions on classical MR images vary in volume in a large range from less than 1 mm 3 to a few cm 3 . A simple linear weighting is not relevant since the volume range is too large for such a mapping and would either not give enough importance to small lesions or to big ones. The proposed weighting method is designed:

• to be decreasing or stagnating with lesion size,

• to have a configurable decay to adjust the small and big lesions weight contrast,

• to have limited sensitivity to small size variations in the beginning nor to separate big from very big lesions from a certain threshold.

The designed weighting method is based on logistic functions for their interesting properties and "S" shape. Weighting function ω : R + → R + is defined as:

ω(v les ) = w max - (w max -w min ) 1 + α × e -kv les range , (5)
where v les is the lesion volume, w max is maximum weighting value, w min is the minimum weighting value, α sets the x-axis translation of the curve, range ∈ R + sets the v range between the minimum and maximum asymptotes of the curve and k defines the steepness of curve. Weighting function is illustrated in Figure 1.

w max -ω(0) = ω(range) -w min =⇒ α = √ e k .A max (6)
In practice, w max , k and α has been considered as constants and the two parameters of interest are:

• w min which sets the weighting contrast between small and big lesions,

• range that sets the decreasing speed and the threshold between small and big lesions.

In order to improve the lesion detection, and small lesion detection in particular, parameters w min and range of the weighting function are investigated in this study.

Loss function

The loss function used is basically a weighted Dice loss function calculated at a voxel level. To weight the lesion voxel depending on the lesion volume, the mapping of weighting function Ω is firstly generated assigning for each lesion voxel position the value of ω(v les ) where v les is the volume of the lesion in which the voxel belongs. Weighted Dice loss function is written:

W D(P, P, Ω) = -2 B b=1 H h=1 W w=1 D d=1 P b,h,w,d Pb,h,w,d Ω b,h,w,d B b=1 H h=1 W w=1 D d=1 P b,h,w,d + Pb,h,w,d , (7)
where B ∈ N + is the batch size, H ∈ N + the height, Ω ∈ N + the width, D ∈ N + the depth, P ∈ {0, 1} BHW D the ground truth segmentation and P ∈ [0, 1] BHW D the model prediction. The loss function, by definition, gives more importance to matched voxel with a high weight. It, then, leads the model to increase its sensitivity for high-weighted voxels that are voxels belonging to small lesions in this case.

CNN architecture

The U-net architecture, originally used for segmentation in electronic and photonic microscopy [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] has been extensively used for segmentation tasks and particularly in medical images. It has become a reference architecture for pixel-and voxelwise segmentation and has been successfully adapted for MS lesion segmentation by multiple research teams. The top-performing methods of MS lesion segmentation (from Zhang et al., [START_REF] Zhang | All-net: Anatomical information lesion-wise loss function integrated into neural network for multiple sclerosis lesion segmentation[END_REF] Isensee et al., [START_REF] Isensee | nnu-net: a self-configuring method for deep learning-based biomedical image segmentation[END_REF] Kang et al. [START_REF] Kang | Acu-net: A 3d attention context u-net for multiple sclerosis lesion segmentation[END_REF] ) on the ISBI MS segmentation challenge, [START_REF] Carass | Longitudinal multiple sclerosis lesion segmentation: resource and challenge[END_REF] use U-net like architectures.

We used the MPU-net++2CBND architecture, [START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF] for its good segmentation performance on small datasets [START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF] and its fast convergence because of its lightness.

Dataset

In-house dataset

We used an in-house dataset from the university hospital of Poitiers gathered by the I3M laboratory as dataset for the lesion size study. This dataset is used for the study since it contains high resolution FLAIR images and to be able to compare the segmentation prediction with the radiologist who annotated the images. In this way, the radiological feedback of the method comes from the annotator itself that can confront the prediction and its segmentation. The dataset contains FLAIR cerebral MR images from 35 different MS diagnosed patients of the university hospital of Poitiers. Only FLAIR images were gathered since it is the most clinically valued sequence for MS lesion detection and also the most important for learning. [START_REF] Feng | A self-adaptive network for multiple sclerosis lesion segmentation from multi-contrast mri with various imaging protocols[END_REF] The FLAIR images have been acquired

sagittally with a resolution of 1 × 0.5 × 0.5 mm 3 resolution by a Verio and a Skyra Siemens MR scanners with a 3 Tesla magnetic field. The segmentation masks were segmented slice by slice for each exam by the same radiologist on not preprocessed FLAIR images.

The training set is composed of 20 FLAIR images and the testing set by 15 FLAIR images. The size of the training set is enough with the light architecture used since 10 images are sufficient. [START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF] The MR images have been carefully distributed in each set equitably regarding the lesion profile (ie the lesion sizes, locations, density) to make training and testing set as much representative and comparable as possible.

For the training and testing images were preprocessed by 1 mm 3 isometric resampling to work on cubic voxel with a good compromise between memory space and lesion size. Images are skull-stripped with the BrainSuite extractor tool, histogram matched [START_REF] Nyúl | New variants of a method of mri scale standardization[END_REF] and the intensities are centered and reduced to having comparable and in a small range around zero.

Training and testing

Since the best set of parameters {w min , w max , range, k, α} is unknown, multiple configurations have been evaluated. In practice, parameter w max has been set to w max = 10, to ensure the loss to be big enough when w min is close to 0, given the rare occurrence of small lesion voxels, to avoid gradient issues due to the float precision when the model begins to be trained. The k parameter has been fixed in order to be close enough to each asymptote when v les = 0 and v les = range. We set k = 10 to be close to the asymptote at those points with a distance of ≈ 0.06. Then α is set to √ e 10 = e 5 following Equation 6.

In total, 32 combination of weighting parameters w min and range have been tested, the corresponding profiles of the corresponding weighting function are provided in Figure 7 of the appendix.

Each learning is initiated with the exact same model with the MPU-next architecture initialized uniformly [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] with the same values in order to evaluate the trained model performances without the initialization factor. Following previous work, [START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF] the model is trained with 32 × 32 × 32 patch with a 2, 048 batch size. The models are trained with 6, 144 patches randomly extracted in the brain volume of each training exams each epoch and are trained 40 epochs. The Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] is used with a learning rate of 0.0004.

The predictions are generated per patch regularly spaced by 8 voxels in each spatial direction and spatially averaged to obtain the predicted image. The prediction map is resampled and registered to the original resolution, it is then binarized with a threshold of 0.5 before being evaluated.

Every test segmentation prediction is evaluated with the Anima segmentation performance analyzer script [START_REF] Baecher-Allan | Multiple sclerosis: mechanisms and immunotherapy[END_REF] .

Results and analysis

The problem of detecting small lesions is not obvious in most evaluations because most segmentation are evaluated regardless of lesion size and metrics are mostly aggregated into statistical descriptors calculated by patient with truly distinct lesion profiles. We propose here to deepen the problem of lesion detection and segmentation through analysis at the patient and lesion level. In this section, we show how it is important to consider the lesion size and discuss the improvements brought by our method.

The lesion size

The size-adaptative method is based on the observation that the described model trained with the common Dice loss gives an acceptable Dice score but a bad F1 score when averaged by patients.

So, even if the voxel-wise segmentation is acceptable, the quality of lesion detection is low. The performances vary a lot from an exam to another and are particularly low for patients with small and few lesions. In Figure 3(a) 2% of lesions are less than 3mm 3 , 10% of lesions have a volume in [3, 5], 17% of lesions have a volume in [5, 10] and 55% of lesions have a volume inferior to 26 mm 3 in the test set. This indicates that most of the lesions are really small lesions with a volume inferior to 26 mm 3 . But, compared to the voxel representation of those lesion in Figure 3(b), the small lesions with a volume inferior to 26 mm 3 represents less than 4% of voxels whereas very large lesions with a volume superior to 1000 mm 3 represent 63% of lesion voxels whereas they represent only 3% of lesions. This phenomenon is a problem since segmentation losses is computed at a voxel level and results in a over-representation of big lesion voxels.

[ 0 , 3 ] [ 3 , 5 ] [ 5 , 1 0 ] [ 1 0 , 2 5 ]

Size-adaptative method

The size adaptative method is always compared to the performances of the model learned without weighting, with the Dice loss.

Patient-wise

Dice F1 LTPR LFPR 0.5990 0.4883 0.4490 0.4236 Table 1 Average performances of the model trained with the simple Dice loss for an exam of the in-house test set. LTPR is the lesion detection sensitivity and LFPR is the lesion false positive rate.

In Figure 4 and Table 1, we can see the voxel-wise segmentation performances in Dice score and detection performances in terms of F1 score, LTPR and LFPR averaged per patient.

For all metrics, we observe different behaviors when w min ≤ 0.01. For each metric, when w min ≤ 0.01 and range ≤ 100, the worst performances are obtained, lower than the scores obtained without weighting in Table 1. This indicates that the rapid and abrupt separation of large and small lesion weights (see Figure 7 in appendix) leads to an overall decrease in lesion detection and segmentation quality. When the range increases with the same values of w min , we observe an improvement which shows that by keeping the same values of weights, but spreading them more over the size range, the model learns to detect more lesions and segment them better. For each metric, the color blue indicates that the value is inferior to the corresponding metric of the model trained without weighted Dice loss (Table 1) and the red color indicates that it is superior.

In Figure 4, we observe that by decreasing w min to 2, and 1, and decreasing range to 100 and 50 there is an increase in Dice score, F1 and LTPR while there is no real trend on LFPR. The LFPR is even higher than without weighting for almost every weighting function evaluated. This observation remains consistent with the weighting which is designed to improve the detection of small lesions that are generally more difficult to spot and more easily confused with healthy brain structures.

The weighting functions tested in Figure 4 only slightly improve the Dice score up to 0.6103 for w min = 2 and range = 50 compared to 0.5990 without weighting in Table 1. Voxel segmentation is thus slightly improved, however, lesion detection is drastically improved with an F1 score up to 0.5943 for the same values of w min and range compared to an F1 score of 0.4883 without weighting. The gain in F1 score is explained by the drastic increase in LTPR brought by the method up to 0.7289 on average for w min = 1 and range = 50. On the other hand, the LFPR is slightly higher between 0.45 and 0.5 for the best F1 scores against 0.4236 without weighting.

This behavior is consistent with the objective of method to detect more lesions and particularly the smallest ones. In fact, without weighting, the learned model detects less than half of the lesions with a LFPR close to 42% whereas with the right set of weights we can detect more than 68% of the lesions or even 73% in average per patient with a higher LFPR of 3% to 8%. In these cases, weighting increases drastically the lesion detection with an F1 score around 0.59 against 0.49 without the proposed method.

On average (per patient), the best segmentations and lesion detections are achieved for w min = 2 with range = 50 with better (lower) LFPR also and for w min = 1 with range = 50 for better TPR. Matching these results with the weighting model, we have the best detection performance by weighting small lesions 5 to 10 times more than large lesions with a progressive separation of lesions by volume up to 50 mm 3 . ter LTPR. It means that more small lesions are detected with the method. However, the LFPR, ie the rate of false lesions detected among all predicted lesions, is slightly higher in general, but not as far as the sensitivity is improved. It means that, with the weighting method, the model learns a new lesion detection compromise between sensitivity and precision and detects far more lesions at the cost of detecting some not existing ones.

Lesion-wise

The F1 score, more commonly used to assess the detection performance, is the harmonic mean of LTPR and PPV (1 -LFPR). The F1 score of best methods depending on the lesion volume is illustrated in Figure 5(c). In this figure, the benefit of the proposed weighting method in lesion detection is mostly for small lesions with volume in [3, 25] mm 3 . However, compared to the distribution of lesion size in Figure 3(a), the lesions with volume in [3, 25] mm 3 represents more than 52% of lesions. So, the proposed method improves the detection of small lesions that represents more than half of lesions. The smallest lesions with a volume less than 3 mm 3 do not have a good F1 score since the LFPR is very high for those lesions despite the good LTPR up to 75% depending on weighting parameters range and w min .

To summarize, the proposed method, with the right set of parameters, improves the detection of small lesions, and particularly the detection sensitivity. Small lesions are the most represented lesions and also the more difficult to spot for a radiologist. The proposed method, then, helps to improve the value of automatic MS lesion segmentation tools for radiologists where it is the most valued.

Radiological confirmation of the method

The radiological assessment of the method was performed on multiple test patients with different weighting parameters sets. For clarity, we have illustrated the observations with the example in Figure 6.

In this figure, we observe, for this example, at first that the unweighted segmentation (the reference) and the radiologist segmentation are very close, except for a small lesion not segmented by the model near the gray matter of the right brain cortex. Some parts are over-segmented and others under-segmented compared to the reference segmentation. When w min = 0.1, we observe that the largest lesions are not segmented at all, indicating that the smallest lesions are overly favored during training. However, with range augmenting from 50 to 1,000, the model tends to segment bigger lesions which is consistent with weighting model.

The segmentation for w min = 6 is very close to the unweighted segmentation which is consistent with the results in Figure 4 and is explained by only a small difference in weighting between different lesion sizes.

For w min = 2, new small lesions are segmented that are not on the radiological segmentation. However, after verification, these are edges of lesions that were not segmented, very small lesions that were not seen on the radiological segmentation and one false positive lesion. With w min = 2, we observe that the increase of range from 50 to 1000 eliminates the smallest segmented lesions, which remains consistent with the weighting function, indeed, the passage from large to small lesions is more progressive and spread out, which decreases the importance of very small lesions.

From observation of the predicted segmentations, we see that with w min values approaching 0, there is a tendency to ignore large lesions and prefer only smaller ones. For w min values at 2 and 1 and range values at 50 and 100, we observed very good segmentations with some false negatives.

Several small lesions that are not detected without weighting are detected with this weighting.

More generally, when observing the different segmentations, we found that the best performing models tend to segment more widely at the edge of lesions than the reference radiologist segmentation, thus lowering the Dice score without influencing the radiological analysis of the segmentation itself. The noisiest images were also the images on which the models performed the least well, however, the models were not trained with data augmentation which could have alleviated this problem, at least in part. The learned models segmented fewer lesions at the base of the skull.

In fact, this area contains artifacts that creates local hyper-intensities that are not lesions and it contains few lesions in general. The models have therefore learned to recognize this area and to avoid to segment any hypersignals in this volume.

A few lesions have been confirmed to be missed by our radiologist but segmented by the model. This last observation points out an interesting fact: an automatic segmentation model, although imperfect, can help a practitioner by detecting lesions that are not very prominent or too small to be obvious, especially with our weighting method.

Evaluation on the ISBI-2015 challenge dataset

The method has been trained on the ISBI-2015 challenge dataset and evaluated by the organization on their test set. The models were trained with every available MR sequence provided (T1, T2, FLAIR and Proton Density) in order to stay consistent with the challenge and to allow the model to work not only with upsampled FLAIR image to grasp the brain anatomy in working resolution. Data augmentation has also been performed in order to help the model to generalize and to improve prediction on unseen data domain since the test dataset is not acquired with the same machine.

The data augmentation performed consists in random flip in 3 dimensions, random blur, random gaussian noise and random anisotropy as implemented in the TorchIO library. [START_REF] Pérez-García | TorchIO: a Python library for efficient loading, preprocessing, augmentation and patch-based sampling of medical images in deep learning[END_REF] The method is evaluated and compared on the ISBI-2015 challenge dataset on Table 2. Our method outperforms the other known published attempts in terms of lesion sensitivity with a LTPR of 0.623 for w min = 1 and range = 100. However, the LFPR reaches 0.457, F1 score is 0.521 and the Dice score is 0.576 which are outperformed by most of other recent approaches in Table 2.

Compared to the unweighted method, the weighted MPU-next with w min = 1 and range = 100 allows a substantial gain in Dice score (+0.060), F1 score (+0.106), LTPR (+0.301) for a loss of 0.215 in LFPR. Surprisingly, the other evaluated weighting parameters ({w min = 1, range = 50}, {w min = 2, range = 50}, {w min = 2, range = 100}) did not perform as well as in Figure 4 on the in-house dataset used for the rest of the study. This can be explained by the fact that the in-house dataset and the ISBI-dataset have a really different lesion size distribution and so the optimal weighting setting is different. As a matter of fact, the in-house dataset has been segmented on FLAIR images acquired with a fine resolution of 1 × 0.5 × 0.5 mm 3 against an acquisition resolution of 0.8 × 0.8 × {4.4, 2.2} mm 3 for the ISBI-2015 training FLAIR images which results in augmenting the minimal discernible lesion size.

Discussion

The automatic MS lesion segmentation can be included into CAD systems for helping radiologists in a daily, tiring and time-consuming task. Given, the context of the disease, it is radiologically more valuable to detect any lesion suspicion than to be sure on prediction. Indeed, it is easier to discard a false detected lesion (such as artifacts or non pathological thin hyperintensities in con-tact with the ventricles) than to spot undetected lesions. In clinical practice, it is essential for the follow-up not to miss any new lesion, even if it is very small, because that may change the medical care. Moreover, the detection is crucial when assessing the spatial and temporal dissemination of lesions that are required for the diagnostic. However, most performing methods for MS lesion segmentation, rely on voxel-wise classification losses that give more importance to big lesion voxels even if small ones have a better radiological value. It results, most of the time, in a competitive voxel-wise segmentation with a limited capacity to detect small lesions. In this work, we focused on giving more importance to small lesions, during model training, in order to improve small lesion detection for a better radiological value.

We showed that, with our weighting method, we can significantly improve the small lesion detection and the lesion detection in general since small lesions are the most frequent. The lesion detection is improved in terms of sensitivity with a LTPR improvement of 0.28 on our in-house dataset and of 0.30 on the ISBI-2015 challenge test set. This improvement is accompanied by an improved Dice score of about 0.01 on our in-house dataset and 0.06 on the ISBI-2015 test set and an increased F1 score from 0.49 to 0.59 on our in-house dataset and from 0.415 to to 0.521 on the ISBI test. The LFPR, however, increases up to 0.8%.

Small lesions are the more difficult to discover, especially when the lesion load is important.

Furthermore, the new lesions on successive MR exams are the most crucial to detect in spite of their often limited size. The presented method has an important radiological significance given the applied context, since less perceptible lesions are more spotted and the overall sensitivity is improved. The high LFPR on smallest lesions can be partially explained by the existing difficulty to detect those lesions by radiologists and consequently by the absence of some of these small lesions in the segmented ground truth. Cases of missing segmented lesions in our in-house dataset have been observed resulting in wrong false positive leading to a biased increase of LFPR. In the same way, in the ISBI-2015 dataset, some lesions, especially small ones, were spotted by only one radiologist segmentation over the two available pointing out the difficulty of spotting small lesions.

The weighting method has, however, to be adapted depending on the dataset and acquisition resolution. More the resolution is fine more the small lesions have to account in the segmentation loss of the learning model. The proposed weighting can be easily reused with other models and methods in order to improve the lesion detection and the radiological contribution for a better adoption of such tools for a better patient care.

Fig 1

 1 Fig 1 The lesion size weighting function ω(v les ) with w max , w min , range parameters illustrated

Fig 2

 2 Fig 2The architecture of the MPU-next model, also referred as MPU-net++2CBND in a previous article.[START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF]

2. 5

 5 .2 ISBI-2015 dataset The ISBI 2015 MS segmentation challenge 132 dataset has been used to compare and validate the method with other existing methods thanks to the online open submission. The ISBI training dataset is composed of 21 preprocessed MR exams from different time points (4 to 5) of five different patient. The exams consist in T1, T2, FLAIR and Proton Density (PD) images. Ground truth segmentations from two different radiologists are provided for each exam. The testing dataset is composed by 61 exams from 14 different patients acquired in another MR scanner. For those examples, the ground truth segmentation is not provided. All the images have been preprocessed by the challenge organization including registration steps and 1 × 1 × 1 mm 3 resampling. Note that the FLAIR image acquisition resolution is thicker than the in-house FLAIR images with voxel size of 0.8 × 0.8 × {4.4, 2.2} mm 3 .

Fig 3

 3 Fig 3The frequency of lesions depending on the lesion volume in the in-house dataset in (a) and the frequency of voxels depending on the lesion volume they belong to (b). The volume range is not linear to avoid flatten representation due to a large volume range. The frequencies are calculated on the test set to be comparable with other results.

Fig 4

 4 Fig 4 Average performances of the model trained with weighted Dice loss with the set of tested parameters w min and range for an exam of the in-house test set. Voxel-wise segmentation is measured in Dice score (a), the lesion detection is measured in F1 score (b), in LTPR (c) and in LFPR (d).For each metric, the color blue indicates that the value is inferior to the corresponding metric of the model trained without weighted Dice loss (Table1) and the red color indicates that it is superior.

Figure 5 Figure 5 Fig 5

 555 Figure5shows the detection capacity of the model trained with the best sets of weighting pa-

Fig 6

 6 Fig 6 Visualization of ground truth segmentation and predicted segmentations on an example of the in-house dataset without weighting and with weighting functions with range ∈ {50, 1000} and w min ∈ {0.1, 6}.

Fig 7

 7 Fig 7 Profiles of the tested weighting functions ω. w max is set to 10 and only w min and range are varying. range sets the distance between the maximal and minimal asymptotes and w min gives the minimum weighting value associated to the largest lesions.

Fig 8

 8 Fig 8 The lesion detection sensitivity (LTPR) depending on the lesion size with all tested weighting functions ω with varying w min and range in blue and without weighting in red. Computed on the in-house test set.

Fig 9

 9 Fig 9 The lesion detection false positive rate (LFPR) depending on the lesion size with all tested weighting functions ω with varying w min and range in blue and without weighting in red. Computed on the in-house test set.

Fig 10

 10 Fig 10 The lesion detection F1 score depending on the lesion size with all tested weighting functions ω with varrying w min and range in blue and without weighting in red. Computed on the in-house test set.

  computed

	Method	Dice	F1	LTPR LFPR
	Weighted MPU-next w min = 1, range = 100 (ours) 0.576 0.521 0.623 0.457
	Multi-view CNN 26	0.627 0.533 0.568 0.498
	All-Net 14	0.639 0.663 0.533 0.122
	nn-Unet 18	0.679 0.645 0.523 0.159
	Weighted MPU-next w min = 2, range = 100 (ours) 0.561 0.505 0.495 0.356
	Ensembling Models 27	0.622 0.645 0.491 0.151
	Multi-Dimensional GRU 28	0.629 0.605 0.487 0.201
	Geo-loss 29	0.643 0.618 0.480 0.132
	Weighted MPU-next w min = 1, range = 50 (ours) 0.463 0.524 0.476 0.293
	Location-aware CNN 30	0.501 0.426 0.429 0.577
	Imagine focal 31	0.584 0.569 0.414 0.087
	Multi-branch CNN 32	0.611 0.556 0.410 0.139
	Cascaded 3D CNN 33	0.630 0.512 0.367 0.153
	Weighted MPU-next w min = 2, range = 50 (ours) 0.449 0.451 0.365 0.255
	MPU-net 34	0.632 0.429 0.347 0.347
	MPU-next without weighting	0.516 0.415 0.322 0.242

Table 2

 2 Comparison of average performances on the ISBI-2015 challenge dataset with known published attempts and our method with weighting functions with range ∈ {50, 100} and w min ∈ {1, 2} and without weighting ordered by LTPR.

available at https://smart-stats-tools.org/lesion-challenge

leaderboard available at https://smart-stats-tools.org/lesion-challenge

https://anima.irisa.fr/

Acknowledgments

This work will not have emerged without the implication of Christophe Habas who passed away in May 2022 and we thank him for his contribution on the whole work built together.

We thank Siemens Healthineers that supported this work.

Appendix A: Additional figures and results

The appendix contains the profile of all tested weighted function ω in Figure 7. It also gather the performances in terms of LTPR (Figure 8), LFPR (Figure 9) and F1 score (10) when training with those weighted functions depending on the lesion size on the in-house dataset.

Disclosures

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The architecture of the MPU-next model, also referred as MPU-net++2CBND in a previous article. [START_REF] Fenneteau | Cnn for multiple sclerosis lesion segmentation: How many patients for a fully supervised method?[END_REF] 3

List of Figures

The frequency of lesions depending on the lesion volume in the in-house dataset in weighted Dice loss (Table 1) and the red color indicates that it is superior.