The proposed prototype system consists of three main steps that are executed consecutively (Figure 1). The first step is to localize the ROI (see details in section 2.2) from the input which is a 4CH CINE series. ROI can be chosen for any anatomical structures, that are of interest displayed in 4CH images, such as the RCA. The localization can be performed by neural networks trained for landmark detection tasks. In case of landmark detection, the output can be pixel coordinates that are used for cropping the image to the localized anatomy (see details in section 2.3). The cropped series containing the ROI are used for motion quantification (see details in section 2.4) by performing the elastic image registration (Chefd’Hotel et al., 2002). By taking the median of the magnitude of the deformation fields from consecutive frames, the motion values over time points are calculated. Finally, the RPs are classified by an absolute threshold, defined based on the correlations of the predicted RPs by varying the thresholds with the expert annotations (see details in section 2.5). The frames which have lower motion values than the threshold value, are selected as RPs. The detected RPs are then used to plan the subsequent static image acquisition.

In this work, the landmark detection neural network is built to automatically detect the location of RCA in a 4CH time-resolved series. The densely connected neural network (3-D DenseNet) architecture is trained to regress the x- and y- coordinates of the RCA over time and it is described in detail here. As a preprocessing step, the 4CH CINE series are interpolated to a fixed spatial and temporal size of $224\times 224\times 32$ to be independent from resolution. Further, the min-max pixel intensity normalization was applied to rescale the different intensity range in [0,1]. The 3-D DenseNet proposed in (Huang et al., 2017) was trained under supervised learning. The weights of the network were updated by using the Adam optimizer (Kingma and Ba, 2014) with $\lambda={10}^{-3}$ and the mean-squared-error loss function (MSE) as follows:

where $\hat{y}_{t,n}$ is the predicted pixel coordinates at the time $t$ from $n$ dataset, the $\hat{y}_{t,n}$ ground truth, $T$ the number of frames, and $N$ the number of the datasets. The ground truth is generated in a semi-supervised manner, where the RCA pixel coordinates in the first frame are manually annotated and propagated to the next frames using the deformation fields describing the displacement between $\hat{y}_{t,n}$ and $\hat{y}_{t+1,n}$. The deformation field are generated by using the elastic image registration (Chefd’Hotel et al., 2002). Each frame was then corrected manually if the propagated coordinates were not accurate. The total number of convolutional layers (Conv) is 122, and before each Conv a 3-D batch normalization (BN) (Ioffe and Szegedy, 2015) and rectified linear unit (ReLU) activation functions (Nair and Hinton, 2010) are applied. After the initial 3-D Conv and max pooling operator, the feature maps were forwarded through 4 concatenated dense blocks (DB) and transition blocks (TB). The number of 4 concatenated DBs are set to 6, 12, 24, 16, in which after each DB, a TB is applied. In each DB, 2 Convs, each followed by 3-D BN and ReLU operators with the increase of the feature maps with 12 are applied. Each layer obtains additional inputs from all preceding layers and forwards the feature maps to all subsequent layers. In each TB, a Conv with BN and ReLU followed by 3-D average pooling operator is applied for spatial and temporal down-sampling. As the last step, a global average pooling and $1\times 1\times 1$ convolutional operator are used to regress the coordinates from the extracted features maps. The detailed architectural details can be found in the Table 1.

From the output of the localization task, a ROI can be simply selected from the pixel coordinates. The detected pixel coordinates are transformed back to the original coordinate system, and then the cropping is performed. Given the predicted pixel coordinates of RCA by 3-D DenseNet, the bounding box is defined by taking the minimum and maximum x- and y-coordinates of the points in the coordinate plane from a time-resolved series and calculating the average of these x and y-coordinates. The size of the bounding box is selected based on prior knowledge about the size of the anatomy, chosen as $50\times 50\ mm^{2}$.

The motion values are quantitatively determined using elastic image registration (Chefd’Hotel et al., 2002). Consecutive frames of the CINE series, $s_{t}(x)$ and $s_{t+1}(x)$ for all timepoints $t$, are registered to obtain deformation fields $d_{t}(x)$ such that $s_{t+1}(d_{t}(x))$ minimizes the dissimilarity measure related to $s_{t}(x)$. The motion curve $m(t)$ describing the amount of RCA motion is then computed as the median of the weighted magnitudes of the deformation fields $\left\lVert d_{t}(x)\right\rVert$ as follows:

where $G_{t}(x)$ is a Gaussian weighting function centered at the midpoint of the detected location of the RCA between $\hat{y}_{t}$ and ${\hat{y}}_{t+1}$ at the time point t:

while $p_{t}$ denotes the midpoint of the detected RCA position. This Gaussian weighting ensures that the motion curve represents mainly the motion of the RCA, while still being robust to slight imprecisions of the localization results. The standard deviation was empirically chosen as $\sigma=12$. Figure 2 shows the Gaussian weighting functions overlaid on the anatomical images at each time point, as well as the weighted deformation fields corresponding to subsequent image pairs. The quantification of motion can be considered in different ways, and the detailed analysis of obtaining the RCA motion values can be found in section 3.2.

After the motion curve is quantified, a classification is required to obtain the RP. This classification is performed within a time interval $T_{valid}=[\alpha,\omega]$, i.e. the first $\alpha$ and the last $\omega$ ms of a cardiac cycle are excluded. $\alpha$ accounts for the time needed for preparation pulses before the data acquisition window can start. $\omega$ is a safety margin at the end of the cardiac cycle to account for heart rate variability, i.e. the detected resting phase should not be too late in the cardiac cycle in case subsequent cardiac cycles during the 3-D acquisition are shorter, which can lead to unstable measurements (Kim et al., 2001; Seifarth et al., 2007). Based on the ground truth annotation of RPs, the $\alpha$, and $\omega$ are empirically chosen to be 80 ms. The RPs can be determined with an absolute threshold from the motion values. The frames with motion values lower than the absolute threshold value are assigned as RPs which can be described as follows:

where $\tau$ is the absolute threshold value, and $m(t)$ is the obtained motion value at trigger time $t$. The threshold value is obtained based on sensitivity (true positive rate) and specificity (true negative rate) analysis from manual annotations. The optimal absolute threshold is chosen as the one achieved with best balanced accuracy.

The mean and standard deviation error between the prediction and ground truth of the fully convolutional 3-D DenseNet with 122 layers was 4.6 ± 1.8 mm. The box plot in Figure 5, shows the Distance Error in mm between the $\hat{p}$ and the $p$ in each frame. Robustness results for unseen datasets are presented in two ways, first in Figure 3 which shows the performance on oblique and diverse oriented cases (n=12) and second in Figure 4, which visualized the worst cases. The quantitative localization results of a part of the study dataset (n=9) acquired with breath-hold and free-breathing CINE sequences are shown in Figure 6 above. The first frame of each CINE series and the corresponding RCA cropped series are shown. Each case was visualized with the first frame of the corresponding CINE series marked with a red box showing the position of ROI defined based on the network prediction and beside it with the cropped series enclosing the area with the outer edge of the cross-section of the RCA overlayed with the generated heatmap.

The motion values obtained by calculating the distance between the predicted pixel coordinates in each adjacent time points achieved 61.1% accuracy for 1.5T, and 52.8% for 3T. As the motion quantification analysis in the Table 3 shows, the 50 percentile/median of Gaussian weighting achieved 90.1% accuracy, whereas the accuracy was 87.2% without weighting the deformation field. The median performed 91.0% accuracy for 1.5T and 88.9% for 3T. The best accuracy achieved by taking the mean metric was 89.3%, without Gaussian weighting. The motion values quantified based on the median approach with Gaussian weighting in 9 clinically acquired dataset are shown in Figure 6 bottom.

For each percentile and mean analysis, the threshold $\tau$ selected based on binary classification task is listed in the right column in the Table 3. From the motion values obtained by taking the median, the selected $\tau$ was 0.2. The resulting ROC curve is plotted in Figure 7, and the accuracy over each threshold step is plotted in Figure 7 in the top row, in which the threshold $\tau$ is marked by an orange vertical line. On the below row in Figure 7 shows the confusion matrices of each annotated datasets. On the above one, the performance of the threshold on the testing dataset is shown, while on the below one the performance of the threshold on the study dataset is displayed.

The detailed results about the performance of the predicted systolic and diastolic RP on the study datasets are listed in Table 4, Table 5. The Bland-Altmann plots showing the performance of start and end detected time point for each systolic and diastolic RP is shown in Figure 8. $MAE_{start,end\ windo{w/frame}_{sys}}$ and $MAE_{start,end\ windo{w/frame}_{dia}}$ for 1.5T datasets was 11.8 ± 16.5 ms (0.38 ± 0.53 frame) and 14.2 ± 18.9 ms (0.48 ± 0.63 frame) for 3T datasets. By using the selected $\tau$, the proposed system resulted in 93.4% accuracy, sensitivity at 90.1% and specificity 96.8% for 1.5T and 92.6%, 90.7% and 94.5% for 3T datasets. $MAE_{start,end\ windo{w/frame}_{sys}}$ and $MAE_{start,end\ windo{w/frame}_{dia}}$ was 13.6 ± 18.6 ms (0.46 ± 0.62 frame) when using the datasets independent from field strength for analysis. The accuracy, sensitivity and specificity achieved by the defined absolute $\tau$, <0.2, was 92.7%, 90.5% and 95.0%. The automatically classified RPs resulted in a mean max error of  30 ms, meaning that it deviates by roughly one frame.

The datasets with RPs with less than 30 ms were discarded (n=9). Further, there was no systolic RP annotated by the expert in 14 cases, and no diastolic RP in 12 cases. These phases were excluded from the analysis. The ${{\widehat{RP}}_{start\ window_{sys}}}$ matched with the annotation, or was off one frame in 93.6%, the ${{\widehat{RP}}_{end\ window_{sys}}}$ in 97.5%, the ${{\widehat{RP}}_{start\ window_{dia}}}$ in 93.8% and ${{\widehat{RP}}_{end\ window_{dia}}}$ in 96.3%, respectively. ${{\widehat{RP}}_{start\ window_{sys}}}$ was detected earlier/later than the expert’s annotation in 27.8%/5% and ${{\widehat{RP}}_{start\ window_{dia}}}$ in 20.9%/16.0%. ${{\widehat{RP}}_{end\ window_{sys}}}$ was detected earlier/later than the expert’s annotation in 11.4%/29.1% and ${{\widehat{RP}}_{end\ window_{end}}}$ in 15.2%/14.8%. In average, the ${\widehat{RP}}_{start\ window}$ was selected earlier than the ground truth in 24.3% of the cases and later in 10.6%. In 10 cases, outliers were present, off by 2 or more frames. In 10 cases, there were outliers present for the ${\widehat{RP}}_{start\ window}$, off by 2 or more frames. For the ${\widehat{RP}}_{end\ frame}$, it was off by 2 or more frames in 5 cases.
The range of annotated systolic RP was 61.1 ± 24.1 ms and the predicted range of systolic RP was 75.5 ± 32.9 ms. Further the range of annotated diastolic RP was 156.0 ± 102.1 ms and the predicted range was 158.2 ± 104.3 ms.

In Figure 9, the robustness of the proposed system is shown in which the system was tested in different sequences including a rescan from a single volunteer. The visualized different CINE outputs are acquired in the order from top to bottom and with  5min between the first and last acquisition. The predicted start and end systolic phases were matched with annotation in most cases, except in one CINE sequence, where the systolic RP was detected by the system but not by the expert, and in a repeat scan, the end time point was off one frame. The start diastolic RP was detected one frame earlier in 2 cases, and end diastolic RP was detected two frames off in real-time sequences. In an example case, the automatically detected RPs were used for the later 3-D static cardiac acquisition targeted to the RCA. The 3-D RCA visualization with the automatically classified RPs showed no residual motion artifacts (Figure 10). The computation time of the proposed system was averaged 1.5 seconds.