Convolutional LSTM model for cine image prediction of abdominal motion

The data used in this project was collected from a Philips Ingenia MR-simulator, with the sequence equivalent to Elekta Unity (Elekta Unity MR-Linac, Elekta AB 2023) to emulate MR-Linac imaging functionalities. The study and the data acquisition were conducted in accordance with all federal, state and institutional regulations regarding use of human subjects in research, under and institutional review board -approved protocol. Cine-MR images were acquired in multi-2D mode using a Philips balanced turbo field echo (bTFE) sequence. Unlike the more common interleaved acquisition, the complete data for each slice was acquired before the next slice acquisition started. The bTFE sequence is a special type of steady-state free precession MR technique that enables rapid and high signal-to-noise ratio (SNR) imaging. The contrast in the bTFE sequence depends on the T2/T1 ratio. Solid tissues typically have a small T2/T1 ratio, such as muscle (0.05) and liver (0.08), while fat (0.30) and fluids (0.70) tend to have a higher ratio. The fast acquisition speed and the contrast characteristics make bTFE an ideal sequence for real-time abdominal motion tracking.

The cine-MR images were acquired in a multi-planar setting, capturing images from different orientations (coronal, axial, and sagittal) sequentially at 200 ms acquisition time for each plane. In this study, only one orientation was used in the model, either coronal or axial, considering the location and shape differences of different organs.

The complete datasets were collected in three separate sessions. The first session included seven volunteer datasets, the second session included eight volunteer datasets, and the third session included data from two patients. No breathing abnormalities were found among all the datasets. Table 1 provides detailed information about the datasets. The first column displays a specific ID assigned to each dataset, where 'Dataset 1–7' indicates the seventh dataset of the first session, for example. The second column represents the matrix size of the 2D cine image, and the 'Time Interval' indicates the time gap between two consecutive cine images of the same plane. The fourth, fifth, and sixth columns show the physical spacing of each pixel, the thickness of each slice, and the total number of images in each time series, respectively. The total time length of one time series is calculated by multiplying the number of images by the time interval.

The LSTM neural network is renowned for its proficiency in handling time series data. Essentially, LSTM utilizes a memory cell to convert 'short-term memory' to 'long-term memory' through self-iteration over time steps. This is achieved by controlling multiple gates inside the LSTM cell. The basic structure of an LSTM cell is illustrated in figure 1. The LSTM cell consists of four components: a memory cell state c_t , an input gate i_t , a forget gate f_t , and an output gate o_t , where t denotes time frames. The memory cell state c_t is updated at each time frame. The update process of the LSTM cell is detailed as follows:

Zoom In Zoom Out Reset image size

The input vector of cell at current time stamp t is denoted by x_t , and the current hidden layer vector is denoted by h_t . W, U, and b represent recurrent weights, input weights, and bias respectively. ⨂ denotes Hadamard product (element-wise multiplication). The sigmoid function and hyperbolic tangent (tanh) function are used as activation functions for cell gates and cell input/output.

The input gate i_t determines the information that can flow into the current cell. It is computed using equation (1):

$$\begin{eqnarray}&&{i}_{t}=\sigma ({x}_{t}{U}_{i}+{h}_{t-1}{W}_{i}+{b}_{i}).\end{eqnarray} \tag{ 1 }$$

The forget gate f_t decides which information will be 'forgotten' or discarded by the current cell. It is determined by equation (2):

$$\begin{eqnarray}&&{f}_{t}=\sigma ({x}_{t}{U}_{f}+{h}_{t-1}{W}_{f}+{b}_{f}).\end{eqnarray} \tag{ 2 }$$

The output gate controls the output of information from the current cell. It is updated by equation (3):

$$\begin{eqnarray}&&{o}_{t}=\sigma ({x}_{t}{U}_{o}+{h}_{t-1}{W}_{o}+{b}_{o}).\end{eqnarray} \tag{ 3 }$$

The cell updates its state from c_t−1 to c_t using equation (4) and equation (5):

$$\begin{eqnarray}&&\tilde{{c}_{t}}=\tan h({x}_{t}{U}_{c}+{h}_{t-1}{W}_{c}+{b}_{c})\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{c}_{t}={f}_{t}\otimes {c}_{t-1}+{i}_{t}\otimes \tilde{{c}_{t}}.\end{eqnarray} \tag{ 5 }$$

Finally, the hidden layer vector h_t is updated using equation (6)

$$\begin{eqnarray}&&{h}_{t}={o}_{t}\otimes \tan h({c}_{t}).\end{eqnarray} \tag{ 6 }$$

ConvLSTM is a modification of LSTM that incorporates convolutional structures within its units. The basic unit of the convLSTM network is the convLSTM cell, which follows almost identical equations to the LSTM cell, except that the internal matrix multiplication is replaced by a convolution operation. Unlike the 1D array structure of LSTM, the convLSTM cell can handle data with higher dimensions, such as 2D or 3D data.

The convLSTM network implemented in this study comprises multiple layers of convLSTM cells followed by regular 2D convolutional layers. Based on a previous study (Lin et al 2019b), the number of convLSTM layers was set to 3. The network structure is illustrated in figure 2. The network inputs consist of consecutive cine images, and the output is a single image predicted several frames ahead of the latest input. The specific frame interval used for prediction as well as the number of input images are part of the hyperparameter selection process.

Zoom In Zoom Out Reset image size

Prior to inputting the data into the network, several pre-processing steps were performed on the cine images. First, all the cine images were resized to a dimension of 352 by 352, and specific regions of interest (ROIs) were cropped, with a separate ROI for each organ of interest. The selected ROIs included the stomach, liver, kidney, and pancreas. Coronal images were chosen for the stomach and liver, while axial images were selected for the kidney and pancreas. Detailed information regarding the ROIs is provided in table 2. It is important to note that the visibility and identification of certain organs varied across datasets, resulting in different numbers of datasets available for each ROI, as indicated in the last column of the table. For example, the pancreas is only available in 5 out of the 17 datasets.

After cropping the images, additional pre-processing steps were performed to enhance image quality, including contrast adjustment, median filtering, and other de-noising techniques. Subsequently, the dataset underwent data partitioning ,as depicted in figure 3, where the images were divided into training, validation, and test datasets using a fixed ratio of 8:1:2, respectively. This consistent partition ratio was employed across different datasets to ensure comparability in subsequent analyzes and evaluations. The images are further partitioned into input images, prediction frames, and output image in training data, validation data, and test data. Input images are used to generate the predicted output image. The predicted output image is then compared with the ground truth output image, as well as the last frame of input images, which is referred as image without prediction in the following text.

Zoom In Zoom Out Reset image size

The convLSTM model was implemented, trained, and evaluated using the Tensorflow 2.6.0 framework and the Adam optimizer in a Python 3.9.7 environment. The training process utilized an NVidia A100 GPU with 80 GB memory and an AMD EPYC 7742 64-Core CPU for computational tasks.

Prior to the individual training phase, a hyperparameter selection process was conducted on the Stomach dataset 2–3 for coronal images and Kidney dataset 2–3 for axial images. This dataset was chosen because the stomach in this dataset exhibits a more complex motion behavior compared to other datasets, and the kidney boundaries in the dataset are distinctly delineated.

Multiple hyperparameters were fine-tuned to determine the optimal configuration of the model, including the number of convolution filters, dropout rate, number of epochs, choice of loss function, kernel size, and input size. During the tuning process, each hyperparameter was adjusted while maintaining the other hyperparameters constant to ensure controlled experimentation. The initial hyperparameter setup consisted of 64 convolutional filters with a filter kernel size of 3 by 3 for each of the three convLSTM layers, the SSIM loss function, 30 epochs, 4 frames as prediction interval, 6 input frames, and 0 dropout rate. Once the optimal hyperparameters were determined, the resulting model was referred to as the Optimal Model and was subsequently applied to other datasets through individual training and transfer learning approaches.

In individual training, the training, validation, and test data were sourced from the same dataset. In transfer learning, 10 datasets were used for training and validation, while the remaining 6 datasets were employed for fine-tuning and testing. Transfer learning techniques were employed to transfer the globally trained model to the specific test sets. Therefore, the test set was completely unseen by the model if fine-tuning was not used. However, the model can be fine-tuned when being applied to the new dataset. Due to limitations in the available datasets, only the kidney and liver went through transfer learning phase. Each case was trained five times, and the average performance was calculated.

It's crucial to clarify that the term 'Optimal Model' does not mean the model is universally optimal across all datasets or organs. Separate Optimal Models were chosen for coronal and axial images, recognizing that a model excelling with coronal images might not perform as well for axial ones. The choice of optimal hyperparameter set is influenced by numerous factors, with the primary goal being consistent and reasonable outputs across diverse datasets. This paper does not delve into optimal hyperparameters for each individual dataset.

The details of hyperparameters that we have examined are listed below:

2.4.1. Prediction interval

The number of prediction frames indicates the interval of time that the model attempts to predict. For example, in Dataset 2–3, where the time interval between two frames is 0.58 s, a prediction frame length of 4 frames corresponds to a total prediction time interval of 2.32 s, which is approximately half of the usual human's breathing cycle. The largest motion-error due to respiration may occur at around this time offset as shown in figure 6. In this study, we have tested intervals ranging from 1 to 50 frames. Interval length of 4 frames was used as default when other hyper-parameters were being tested because it represents the most challenging time interval.

Zoom In Zoom Out Reset image size

2.4.2. Loss functions

Various loss functions were employed for comparison, including SSIM, mean square error (MSE), and mean absolute error (MAE). SSIM was used as default.

SSIM is an image comparison metric commonly used to assess the similarity between two images. The SSIM formula is defined as follows:

$$\begin{eqnarray}&&\mathrm{SSIM}(x,y)=\displaystyle \frac{(2{\mu }_{x}{\mu }_{y}+{C}_{1})(2{\sigma }_{{xy}}+{C}_{2})}{({\mu }_{x}^{2}+{\mu }_{y}^{2}+{C}_{1})({\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{C}_{2})}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{C}_{1}={\left({K}_{1}L\right)}^{2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{C}_{2}={\left({K}_{2}L\right)}^{2},\end{eqnarray} \tag{ 9 }$$

where μ_x and μ_y are the mean intensity, σ_x and σ_y are the standard deviation estimation of the images, σ_xy is the covariance between image x and image y , K₁ and K₂ are 0.01 and 0.03, respectively, based on the standard implementation (Wang et al 2004). L is the dynamic range of the pixel values. In actual usages, the SSIM of images are calculated based on the mean SSIM value of various windows of an image, with the default window size 11 by 11.

MSE is defined as

$$\begin{eqnarray}&&\mathrm{MSE}(x,y)=\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{\left({x}_{i}-{y}_{i}\right)}^{2}\end{eqnarray} \tag{ 10 }$$

where x_i and y_i are the pixel grayscale intensity values of corresponding pixels with index i. MAE, similar to MSE, is defined as

$$\begin{eqnarray}&&\mathrm{MAE}(x,y)=\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}| {x}_{i}-{y}_{i}| \end{eqnarray} \tag{ 11 }$$

2.4.3. Number of input frames

2.4.4. Number of epochs

2.4.5. Dropout

2.4.6. Convolution filters

2.4.7. Kernel size

After obtaining the optimal model, we applied it to each dataset for comparison. The ROIs included stomach, kidney, pancreas, and liver. The data partition method was mentioned in 2.3. The output images were saved for comparison with ground truth.

Transfer learning is based on the idea that if a model is trained on a large and comprehensive dataset, then it can effectively learn general features that can be applicable to other tasks. Instead of starting from scratch, transfer learning typically uses a a model pre-trained on a large dataset and then apply it to a new and smaller dataset, with or without fine-tuning. Transfer learning allows for leveraging pre-existing models, reducing both computational cost and training time. Also, it often leads to improved performance, especially in tasks where the available dataset is relatively small.

Transfer learning was implemented for the liver and kidney datasets. We used 10 datasets for initial training and validation, while the remaining 6 datasets were dedicated to testing and fine-tuning. For each training dataset, approximately $\tfrac{8}{9}$ of the data were used as the training set for transfer learning, and the remaining portion served as the validation set. This resulted in a total of 8284 training images and 1036 validation images. After the training process was finished, the trained model was then applied to test set via transfer learning. Regarding the test set, approximately $\tfrac{8}{11}$ of the data were utilized for fine-tuning of transfer learning, $\tfrac{1}{11}$ for validation, and the last $\tfrac{2}{11}$ for final testing. The final test set was used as the sole test set regardless of whether fine-tuning was involved or not.

To assess the performance of image prediction, we utilized several image comparison metrics. The metrics employed in our experiments include structural similarity index (SSIM), normalized root mean square error (NRMSE), NCC, and PSNR.

NRMSE is a normalized variant of MSE ranged from 0 to 1 for better visualization. $\bar{y}$ is the average intensity of the reference image. NRMSE and related normalized mean square error (RMSE) are defined as:

$$\begin{eqnarray}&&\mathrm{RMSE}=\sqrt{\mathrm{MSE}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\mathrm{NRMSE}=\displaystyle \frac{\mathrm{RMSE}}{\bar{y}}.\end{eqnarray} \tag{ 13 }$$

NCC is a normalized index measuring correlation between two variables. The 2D NCC is defined as

$$\begin{eqnarray}&&\mathrm{NCC}(x,y)=\displaystyle \frac{\displaystyle {\sum }_{m}\displaystyle {\sum }_{n}({x}_{{mn}}-\bar{x})({y}_{{mn}}-\bar{y})}{\sqrt{\left(\displaystyle {\sum }_{m}\displaystyle {\sum }_{n}{\left({x}_{{mn}}-\bar{x}\right)}^{2}\right)\left(\displaystyle {\sum }_{m}\displaystyle {\sum }_{n}{\left({y}_{{mn}}-\bar{y}\right)}^{2}\right)}}\end{eqnarray} \tag{ 14 }$$

where NCC is normalized correlation coefficient between two 2D variables x and y, m and n are the row and the column indices, respectively, $\bar{x}$ and $\bar{y}$ are the average value of x and y.

PSNR is defined as

$$\begin{eqnarray}&&\mathrm{PSNR}=10\cdot {\mathrm{log}}_{10}\left(\displaystyle \frac{{\mathrm{MAX}}_{I}^{2}}{\mathrm{MSE}}\right)\end{eqnarray} \tag{ 15 }$$

where MAX_I is the maximum possible pixel value of the image, typically 65 535 for 16 bit images.

Boxplot figures 4 and 5 show the prediction performances with various hyperparameters in order to find the optimal hyperparameter combinations for both coronal and axial image datasets. The boxes in boxplot display results based on the five-number summary: the maximum value, the 75% quartile, the sample median, the 25% quartile and the minimum value. The comparison metric used for selection is SSIM score. If one hyperparameter variation showed a higher score than the default setup and also showed a statistically significant difference, then the hyperparameter was selected instead of the default parameter. The statistical significance was calculated based on p-value with an alpha of 0.05 signifying a significant difference which was obtained by comparing the SSIM score of the test model with the default model. For coronal images, all hyperparameters in the default setup were selected for the optimal model. For axial images, the optimal model was selected with loss function of MAE, input number of 6, and the number of convolution filters of 96, 96, and 96.

Zoom In Zoom Out Reset image size

Figure 6 displays prediction results based on varying time intervals while maintaining other hyperparameters consistent with the optimal model. As time intervals increase, the SSIM score of prediction starts to steadily and slowly drop. In the default model, four frames which roughly equate to 2.4 s in most datasets was used as time interval. It is about half of the human breath cycle time, which represents the largest motion offset due to respiration. In order to prove the model's proficiency in predicting a challenging period of time, an interval of four frames was utilized as time interval in subsequent experiments. Figure 7 offers a visual comparison of the ground truth image (first row), predicted image (second row), and the image without any prediction (third row), as well as their pixel intensity disparities with the ground truth for stomach, liver, kidney, and pancreas. This predicted representation indicates motion projected 2.4 s ahead from the most recent input. The pixel-wise difference between the ground truth and the prediction is noticeably smaller than the difference between the ground truth and the latest input for all ROIs.

Zoom In Zoom Out Reset image size

The model performance was also displayed quantitatively using different metrics besides visual comparison. Tables 3–6 show the model performance scores in stomach, liver, kidney, and pancreas. Despite variations across datasets and organs, our predictions generally enhance the image similarity across most metrics, with only a few exceptions. The average SSIM scores for stomach, liver, kidney, and pancreas are 0.538, 0.641, 0.772, and 0.661, respectively.

Table 7 and table 8 present the prediction results for the kidney and liver datasets using transfer learning. The tables are divided into three sections: the first section shows the prediction results with fine-tuning, the second section shows the performance without fine-tuning, and the last section displays the image similarity without prediction.

The total computation time of model inference is much smaller than current latency time in MR-Linac The computation time to generate one predicted image for smaller ROIs like stomach is less than 10ms, and the computation time for larger ROIs (liver) is around 100ms, given single GPU.

4.1.1. Input image frames

Figure 4 and figure 5 illustrate that for a prediction time interval of 4 frames, the input number of 6 results in the best average SSIM performance for coronal images and the input frames of 10 works best for axial images. This finding suggests that for different organs with different orientations, the number of inputs that is needed to provide sufficient information for accurate motion prediction may be different. The optimal input image frames may be related to the breathing cycle which is about 8–10 frames in this case. For coronal images, in order to predict one image 4 frames ahead, an input number of 6 is more likely to capture the corresponding breathing state of ROI of the object image from the previous breathing cycle, compared to an input number of 5, which sometimes may fail to cover a full breathing cycle. For axial images, the number of inputs for optimal outcome is 10, which is around a full breathing cycle. The difference of the optimal input frames may due to the fact that the motion of ROIs is largely superior–inferior.

Zoom In Zoom Out Reset image size

4.1.2. Convolutional filters

4.1.3. The choice of loss function

4.1.4. Dropout

4.1.5. Other hyperparameters

Figure 6 demonstrates the effect of increasing prediction time interval on performance. It is evident that the model achieves the highest SSIM score when the interval is 1 frame. However, even at larger intervals, beyond approximately 10 frames, the prediction SSIM scores stabilize. Although the scores for these larger intervals are not as high as those for smaller intervals, they still exhibit a higher similarity with the ground truth compared to the corresponding image without prediction. Notably, a discernible pattern is observed in the SSIM score without prediction, indicating the presence of periodic motion in the organs. This pattern is also evident in the liver, pancreas, and kidney datasets (figures 6(b)–(d)), which has a period about 10 frames or 6 s, likely mainly contributed by respiration cycle. The similarity score is higher when the prediction interval time aligns with respiration, while it decreases when the prediction interval time corresponds to the opposite phase of respiration.

Different organs display varying prediction performances, with the kidney generally exhibiting better performance than the liver. Several factors, such as image pixel size, image orientation, and organ motion patterns, could contribute to this disparity in prediction accuracy. Nevertheless, the SSIM score without prediction curve demonstrates a similar pattern across all four organs (figure 6). This suggests that all four organs are at least partially affected by a common type of motion, likely respiratory motion, as indicated by the peak values around frames 9 and 19.

The SSIM score in our study was relatively lower compared with some studies using PCA and other methods. For example, recently, results from a PCA modelling paper using 4DCT achieved SSIM of 0.88 for chest region (Asano et al 2023). For comparison, our SSIM values were 0.5–0.6. However, it is important to note that the PCA result was based on whole uncropped 4DCT chest images, which included mostly stationary anatomy. Our SSIM is based on a cropped ROI which includes mostly moving anatomy, so our lower SSIM values are to be expected. Most importantly, the LSTM model results are for predicted images for future target position rather than current position. As well, the PCA model deals with simpler target motion pattern in lung compared to abdominal motion which is more complex.

To assess the robustness of the model, we evaluated the SSIM scores for all images in the test set (figure 8). The robustness of a test set is defined as the ratio of effective predictions (prediction SSIM score higher than the score without prediction) to the total number of images in the test set. The robustness scores for all test sets are presented in table 9. Notably, the stomach test set achieved a robustness of 79.75%, indicating that the model was effective for around 80% of the images in the stomach dataset. Moreover, one liver dataset used in transfer learning achieved a robustness score of one, indicating that the prediction model successfully improved the image similarity for all images in this dataset. The overall average test set robustness was 82.1%.

Zoom In Zoom Out Reset image size

Table 7 and table 8 compare the prediction results for kidney and liver datasets using transfer learning. The tables are divided into three sections: performance with fine-tuning, performance without fine-tuning, and image similarity without prediction. It can be observed that transfer learning with fine-tuning has a beneficial impact on prediction accuracy. The SSIM scores are consistently higher when fine-tuning is applied compared to the case without prediction, especially in table 7. While the transfer learning model performs relatively poorly in some datasets without fine-tuning, the application of fine-tuning largely compensates for the performance difference.

A primary issue is that our current model necessitates cropping the region of interest (ROI) from cine-MR images. Using full-size images without this cropping results in diminished performance in organ motion tracking. Another limitation is that our model only processes one plane from a three-planar dataset, thereby ignoring potential insights from the remaining planes. Furthermore, the 2D model's intrinsic constraint is its inability to monitor organ motion beyond its current plane. To address these drawbacks, we're in the process of designing a 3D cine-MR motion prediction algorithm, which would allow motion tracking in all three dimensions. Although 3D cine-MR is still not a mature technique and is faced with multiple issues, some promising results have been published (Wu et al 2023). This 3D model development is our imminent research focus.