# Convolutional module for heart localization and segmentation in MRI

D. Lima (1, 2),  C. Graves (2),  M. Gutierrez (2),  B. Brandoli (3),  J.
Rodrigues-Jr (1),
(1) ICMC-USP – Institute of Mathematical and Computer Sciences, University of
Sao Paulo
Avenida Trabalhador Sao-carlense 400, 13566-590, Sao Carlos, SP, BR
(2) Biomedical Informatics Lab, InCor/HC.FMUSP – Heart Institute, University
of Sao Paulo
Avenida Doutor Eneas de Carvalho Aguiar 44, 05403-900, Sao Paulo, SP, BR
(3) Dalhousie University, 1459 Oxford Street, B3H 4R2, Halifax, NS, CA
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Magnetic resonance imaging (MRI) is a widely known medical imaging technique
used to assess the heart function. Deep learning (DL) models perform several
tasks in cardiac MRI (CMR) images with good efficacy, such as segmentation
Bernard2018 , estimation xue2018full , and detection of diseases
khened2017densely . Many DL models based on convolutional neural networks
(CNN) were improved by detecting regions-of-interest (ROI) either
automatically or by hand. In this paper we describe _Visual-Motion-Focus_
(VMF), a module that detects the heart motion in the 4D MRI sequence, and
highlights ROIs by focusing a Radial Basis Function (RBF) on the estimated
motion field. We experimented and evaluated VMF on three CMR datasets,
observing that the proposed ROIs cover 99.7% of data labels (Recall score),
improved the CNN segmentation (mean Dice score) by 1.7 ($p<.001$) after the
ROI extraction, and improved the overall training speed by 2.5 times (+150%).

_Keywords_ cardiac mri $\cdot$ deep learning $\cdot$ motion $\cdot$
localization $\cdot$ segmentation

##  1 Introduction

Magnetic resonance imaging (MRI) is a medical imaging technique used to
capture volumetric image sequences of internal soft tissues, such as cardiac
muscles. In comparison to X-Ray imaging (XR) and Computer Tomography (CT), MRI
provides images with improved structural details via finer spatial
resolutions. Cardiac MRI (CMR) focuses on the heart, allowing trained
cardiologists to measure heart parameters, for example the mass of the cardiac
muscle (myocardium mass), the volumes of blood cavities (atrial and
ventricular volumes) and the amount of blood pumped per heartbeat (ejection
fraction) peng2016 . Those parameters are used to assess how healthy is the
heart, by recognizing early conditions and signs before the onset of infarcts
and other complications.

Due to the size and complexity of CMR sequences, complex techniques are
required to produce detailed analyses; one of these techniques is deep
learning (DL). Many of the tasks and goals related to the cardiac functional
analysis – for example segmentation of structures Bernard2018 , estimation of
heart parameters xue2018full , and detection of diseases khened2017densely  –
have benefited from DL methods. For even better results, research in DL has
pointed out that models based on convolutional neural networks (CNN) have had
a higher efficacy when provided with regions-of-interest (ROI) either
explicitly or implicitly xue2018full . The detection of ROIs, usually named
ROI proposal, is a preprocessing step whose goal is to identify the most
prominent regions of an image (frame) for discovering clinically-relevant
artifacts.

The explicit ROI proposal approaches usually follow a combination of methods,
for example: (a) pipelining a segmentation and a regression network; (b)
preprocessing the input with a region proposal algorithm he2015spatial  or
with a CNN wu2020lv ; or (c) by using manual cropping xue2017full . The
implicit ROI detectors are added to the DL network abstracted as additional
operators and variables; e.g. multi-scale “Inception” szegedy2015going  and
attention vaswani2017attention  modules, which benefit from the ROIs to down-
weight less-informative neurons and inputs inside the network. Inception
modules weight convolutions of different sizes, while attention modules assign
a weight to each feature channel. This additional neural-network information
processing guides which input features or channels shall have more weight so
to improve the accuracy of the output layer.

In this paper we develop a module that highlights regions in the image
sequence by analyzing the motion field using a Radial Basis Function (RBF). In
our experiments we analyze our method by using the RBF for cropping the input
before having it processed by a pretrained segmentation CNN. Our methodology
is an innovation in the task of region proposal for CMR analysis; as we
demonstrate, we achieved results that justify the use and further
investigation of the employed principles. We named it after its working
mechanism as VMF – Visual-Motion-Focus. Section 2 presents the theory and
related works regarding ROI proposal and motion detection. Section 3 details
VMF, presenting its architecture and describing how each step works. Section 4
shows the experiments over VMF, where we demonstrated excellent performance
for cropping the ROI, according to metric _Recall_. Additionally, we also
report our results regarding both _Dice_ and _Speed-up_ metrics. In
conclusion, this work demonstrates the application of technique VMF in
objective tasks and how to achieve such a setting. Section 7 discusses
conclusions and future works.

##  2 Theory and related work

###  2.1 Cardiac MRI

MRI is the most precise medical imaging technique for examination of the heart
structures, it allows the recording of heart images along a complete heartbeat
cycle earls2002cardiac ; seraphim2020quantitative . In practice the
magnetization signal is triggered by a reference pulse, then captured several
times for noise reduction and, finally, reconstructed by inverse FFT. The
resulting image is usually visualized in slices along a positional axis: long-
axis has a frontal or lateral view of the heart, and short-axis aligns to a
cross-sectional plane.

Figure 1: Picture of heart, and examples of short-axis CMR slices in the Base,
Middle and Apex regions. The medial slice shows the right and left ventricles
(RV and LV). LV boundaries are highlighted in the three slices – epicardium is
the outer boundary and endocardium is the inner boundary. Pictures adapted
from the datasets used in iaizzo2013importance ; saber2020multi .

The short-axis view is split in three regions: the base or basal region near
the top where blood vessels connect to the heart (slice B); the middle or
medial in between (slice M); and the apex or apical region at the bottom tip
of the heart (slice A) – refer to Figure 1. The normal human heart has four
chambers: right atrium (RA), right ventricle (RV), left atrium (LA) and left
ventricle (LV). The atria receive blood and the ventricles pump it out of the
heart. Even though all chambers are important, the LV is of special interest
because it is the cardiac muscle that does the “heavy lifting” of pumping
oxygenated blood from the lungs to the whole body. In the short-axis view, the
LV appears as a ring shape, whose thickness and internal volume measurements
are essential to estimate the myocardium mass and ejection fraction,
respectively.

###  2.2 Region proposal and Image segmentation

Region proposal is at the core of image analysis; it is responsible for
determining the direction to where, in the image, the analytical algorithms
should focus at. Region proposal methods range from classical algorithms
(statistical features, signal processing, binary morphology, etc.) to very
complex models such as modern Convolutional Neural Networks (CNNs). The goal
of such algorithms varies from small outputs such as center points or bounding
boxes (e.g. quad-tree segmentation, R-CNN and YOLO spann1985quad ;
liu1994multiresolution ; girshick2014rich ; redmon2016you ) to fine contours
(e.g. active contour models, region-growing segmentation, U-net and Mask-RCNN
ranganath1995contour ; muhlenbruch2006global ; he2017mask ).

In many methodologies, the detected ROIs are the starting point of supervised
learning with labeled data. However, region-proposal is also organized
according to the region under analysis, e.g. image-level, polygon-level, and
pixel-level. In this paper we use the nomenclature from lin2014microsoft : (a)
_Image classification_ assigns labels to the whole image (all of its objects);
(b) _Object identification_ assigns labels to the bounding box region of the
detected objects; (c) _Semantic segmentation_ assigns labels in polygonal or
pixel-level regions, producing finer contours; (d) _Instance segmentation_ not
only detects the objects’ regions but also each individual instance of a
detected label. Thus, in order to clear any confusion, we consider that
_parameter estimation_ is a regression task, and that _image segmentation_ is
a polygon or pixel-level classification task as mentioned above.

###  2.3 Image features and Motion analysis

Some diseases may decrease the extension and intensity of heart movements; in
such cases, the motion features will have little to no information. A way to
address this issue is by combining the static visual features of the image and
the dynamic features. The static visual features are computed by classical
feature extractors, then are adjusted to detect the heart that, in short-axis
CMRs, appears as a gray bubble with dark borders. The dynamic features are
computed by motion estimation.

Many categories of image feature extractors are present in the literature
gonzalez2002digital , such as statistical moments (mean, standard deviation,
skewness, kurtosis), robust statistics (median, quartiles, interquartile
range, percentiles, rank, median absolute error), signal processing filters
(Gaussian, Sobel, Laplacian, Gabor) and transforms (Hough, FFT, DCT, DWT),
among other measures of information distribution (Shannon entropy, Minkowski-
Bouligand box-counting dimension, Hausdorf fractal dimension).

Motion estimation refers to a set of techniques to explore the temporal
redundancy in sequences of images, with applications in video compression,
analysis and processing Oliveira2006 . One classic formulation of motion
analysis is optical flow horn1981determining , which estimates the
displacement field of an object by smooth derivatives of brightness, assuming
that the brightness of moving objects does not change between frames. Other
approaches include statistical and block-based algorithms; change detection by
individual or accumulated frame difference; and parametric methods
Oliveira2006 . Motion estimation is mainly used for general 2D video coding,
camera motion compensation thanou2015graph ; zhang2018fast ; bao2019memc , and
for medical imaging stemkens2016image ; mohsin2017accelerated ; yan2018left ;
yu2020foal .

###  2.4 Radial Basis Functions

A Radial Basis Function (RBF) is a real-valued function of the radius (or
distance) from an origin to a given point majdisova2017radial . They have
several applications in function approximation and neural networks
broomhead1988radial ; buhmann2003radial . RBFs can be understood as a
flashlight beam on a wall, which is intense at the center and fades sideways.
RBF kernels are built by computing the distance of each pixel position from
the center, then transforming this distance by a decaying function (e.g. exp,
log, Gaussian). Many center functions can be used as the RBF center, such as
statistical measures of central tendency (e.g. simple and weighted mean,
median, center-of-mass). Examples of distance functions include the standard
deviation (also known as Malahanobis distance, or z-score difference),
Minkowski (e.g. $L_{1}$, $L_{2}$, $L_{\infty}$), Jaccard, and Cosine
distances.

###  2.5 Computer methods for CMR analysis

Computer methods for functional analysis of CMR as reviewed by Peng et al.
peng2016  were organized in three ways: image-driven, model-driven, and by
direct estimation. Further subdivisions of the LV segmentation spans five
groups: (1) image processing methods such as thresholding, morphology
operators, and region-growing; (2) pixel/voxel-based classification by
gaussian mixture models, neural networks, k-means, k-nearest neighbors, or
support vector machines (SVM); (3) active contours (snakes), deformable
models, level sets, and motion tracking; (4) principal or independent
component analyses (PCA and ICA) with strong priors from statistical models of
the heart anatomy; and (5) direct estimation by, e.g., latent discriminant
analysis (LDA) combined with SVM. This work refers to category 2 as it employs
a pixel-based classification mediated by a neural network.

###  2.6 Deep learning region proposal and LV segmentation

Recent approaches for LV segmentation use CNN models (such as U-net)
experimented over diverse datasets and methodological combinations. U-net
ronneberger2015u  is a general segmentation model which combines a tower of
downscaled-then-upscaled deep representations by concatenation with the input
features, instead of addition like FPN lin2017feature  and LinkNet
chaurasia2017linknet , with success in LV segmentation. The approach
introduced by Smistad et al. smistad2017  departs from the pretraining of an
U-net CNN using a Kalman-filter segmentation; it achieves 86% Dice
(epicardium) in ultrasound images – Dice (DSC) will be detailed in Section
4.1. For CMR, U-net displayed even better results (89% Dice) when trained from
min-cut priors guo2018cardiac . An U-net architecture with residual blocks and
optical flow information yan2018left  achieved 89%, 95% and 85% Dice in the
base, middle and apex regions of the heart respectively. In the work of Wu et
al. wu2020lv , the authors combine a custom CNN for region proposal with an
U-net segmentation to achieve 95% Dice. Overall, the combination of region
proposal to U-net had good results in particular datasets, but still has room
for improvement when the evaluation generalizes across multiple datasets.
Different to former works, our methodology, VMF, uses RBFs to propose ROIs
that will aid a CNN processing in the task of LV segmentation. As we
demonstrate, we achieve superior results over an ample set of experiments.

##  3 Material and methods

VMF starts with a 4D image input $\mathbf{x}=I(t,x,y,z)$, that is, a sequence
of volumes (frames) each one with a time $t$ coordinate. Initially, we
normalize $\mathbf{x}$ to produce sequence $\mathbf{x}^{*}$ with frames in a
format more adequate for Neural Network processing – see Section 3.1. From
$\mathbf{x}^{*}$, we extract visual features to produce $\mathbf{x}_{s}$, and
motion estimation to produce $\mathbf{x}_{t}$, detailed in Section 3.2. Next,
we apply two sets of weights: $w_{s}$ the weights related to visual features;
and $w_{t}$, the weights related to the motion, or time; combining both
features in tensor $\mathbf{v}$. Then, as presented in Section 3.3, we compute
the center voxel $\mu_{v}$ defined by the energy-weighted sum of all the
voxels’ coordinates; then we produce a segmentation map $\mathbf{y_{S}}$ by
applying a threshold to $\mathbf{v}$ and extracting the bounds of non-zero
voxels; afterwards we compute the scale $\sigma_{v}$ given by the standard
deviation of all the voxels’ distances from center $\mu_{v}$. At this point,
refer to Section 3.4, we can apply a Radial Basis Function at center $\mu_{v}$
with radius $\sigma_{v}$, computing $\mathbf{y_{L}}$, then we scale the region
defined by $\mathbf{y_{L}}$ to the CNN input shape. In Section 3.6, we explain
the Neural Network processing, its parameters and training.

###  3.1 Intensity normalization

We normalized the image intensities between 0 and 1, by subtracting the
minimum value and dividing by the range of values, where $\epsilon$ is a small
constant to avoid division-by-zero:

$\mathbf{x}^{*}=\frac{\mathbf{x}-\mathrm{min}(\mathbf{x})+\epsilon}{\mathrm{max}(\mathbf{x})-\mathrm{min}(\mathbf{x})+\epsilon}$
(1)

###  3.2 Visual features and motion estimation

Figure 2: Sequence of images in the medial slice, showing six frames (upper
row) and the respective absolute derivatives in the time dimension (lower
row). The example in the figure only exhibits the medial slice because it is
easier to visualize and interpret, but notice that VMF is defined for volumes.

For visual features extraction, we consider the statistical mean and standard
deviation obtained with 3 $\times$ 3 $\times$ 3 kernels. The mean image
$I_{\mu}$ is computed using the convolution operation with mean kernel $M$,
defined as:

$M=\frac{1}{1\cdot 3\cdot 3\cdot 3}\cdot J_{1,3,3,3}\\\ $ (2)

where $J_{1,3,3,3}$ is a 1 $\times$ 3 $\times$ 3 $\times$ 3 matrix-of-ones.
That is, $M$ is just the arithmetic average of a 3 $\times$ 3 $\times$ 3
matrix arranged for convolution. With kernel M, we compute $I_{\mu}$ as:

$I_{\mu}=I\ast M$ (3)

In turn, the standard deviation image $I_{\sigma}$ is computed by taking the
differences between the image $I$ and mean image $I_{\mu}$, then squaring the
differences element-wise (Hadamard power), after that by convolving with the
mean kernel $M$, and then taking the Hadamard square-root, as follows:

$I_{\sigma}=[(I-I_{\mu})^{(2)}\ast M]^{(1/2)}$ (4)

Accordingly, the mean image $I_{\mu}$ and the standard deviation image
$I_{\sigma}$ define $\mathbf{x}_{s}=\\{I_{\mu},I_{\sigma}\\}$. Notice that we
express the Hadamard powers using the definitions and notations defined in the
work of Fallat and Johnson fallat2007hadamard .

We used the motion estimate function $E(I)$ – refer to Equation 5, given by
the root-mean-squared differences of intensity along the time coordinate,
where $T$ is the time interval (or number of frames) in image $I$, and $S_{t}$
is the Sobel kernel w.r.t. time, instead of the default $S_{x}$ and $S_{y}$
spatial Sobel kernels. This function is related to the magnitude of the
optical flow horn1981determining  vector field in each voxel, as follows:

$\mathbf{x}_{t}=E(I)=\sqrt{\frac{1}{T}\int_{t=0}^{T}\left(\frac{\partial
I}{\partial t}\right)^{2}dt}\approx\sqrt{\frac{1}{T}\sum_{t=0}^{T-1}[I(t)\ast
S_{t}]^{(2)}}$ (5)

a)  $E(I)$ b)

Figure 3: a) Motion estimate $E(I)$ along the whole cardiac cycle. b)
Cumulative histogram of the computed values. The histogram shoulder lies near
the percentile $p=0.9$, which is used for thresholding.

Figure 2 illustrates a sequence of six CMR frames spaced by $1/6$ of the
cardiac cycle (upper row in the figure), and the respective absolute
derivatives in each point along the time dimension (lower row) as computed
with Equation 5. Figure 3a shows the total motion estimate in the sequence,
while Figure 3b presents the cumulative energy histogram. By combining both
static and motion features $\mathbf{x}_{s}$ and $\mathbf{x}_{t}$, we obtain:

$\mathbf{v}=w_{s}\mathbf{x}_{s}+w_{t}\mathbf{x}_{t}$ (6)

The weights $w_{s}$ and $w_{t}$ are initialized to $0.1$ and $0.9$
respectively for $\mathbf{x}_{s}$ and $\mathbf{x}_{t}$, i.e., although we
emphasize the motion features, we also include the static visual features,
which will address the problematic cases when the heart has limited motility,
which might be the case for heart complications.

###  3.3 Center and scale computation

The center of energy $\mu_{v}$, formalized by Equation 7, is the voxel defined
by the energy-weighted sum of all the voxel coordinates $\mathbf{r}$ in
dimensions $x,y,z$; that is, every voxel $i$ has $\mathbf{r}_{i}=\langle
x_{i},y_{i},z_{i}\rangle\ |\ 0\leq x_{i}<\mathrm{width},0\leq
y_{i}<\mathrm{height},0\leq z_{i}<\mathrm{slices},1\leq i\leq N$, for $N$ is
the number of voxels in the image, as follows:

$\mu_{v}=\langle\bar{x},\bar{y},\bar{z}\rangle=\frac{\sum_{i=1}^{N}\mathbf{v}_{i}\cdot\mathbf{r}_{i}}{\sum\mathbf{v}}$
(7)

The segmentation $\mathbf{y_{S}}$ is given by the thresholded image from the
previous step, see Figure 4. The scale estimate $\sigma_{v}$ is the cube root
of the volume (in voxels) of the thresholded image considering the values
above the 90% percentile:

$\displaystyle(\mathbf{y_{S}})_{i}$
$\displaystyle=\mathbbm{1}[\mathbf{v}_{i}>Q(0.9,\mathbf{v})]$ (8)
$\displaystyle\sigma_{v}$
$\displaystyle=\frac{3}{\|\mathbf{r}_{\max}\|}\sqrt[3]{\sum_{i=1}^{N}\nolimits(\mathbf{y_{S}})_{i}}$
(9)

where the indicator function $\mathbbm{1}[c]$ returns $1$ if the condition $c$
is true, or $0$ otherwise; and $Q(p,v)$ is the quantile function, returning
the maximum of the lowest $p$ (%) values in $v$.

###  3.4 Segmentation and localization

The segmentation $\mathbf{y_{S}}$ is derived from thresholded $\mathbf{v}$
(Equation 8). The localization $\mathbf{y_{L}}$ is found by fitting an RBF to
$\mathbf{y_{S}}$. The radius $d_{i}$ is the distance from the center $\mu_{v}$
to each voxel $i$. The Euclidean distance ($L_{2}$-norm $\|\cdot\|$) was used:

$d_{i}=\left\|\frac{\mathbf{r}_{i}-\mu_{v}}{\mathbf{r}_{\max}}\right\|=\sqrt{\left(\frac{x_{i}-\bar{x}}{x_{\max}}\right)^{2}+\left(\frac{y_{i}-\bar{y}}{y_{\max}}\right)^{2}+\left(\frac{z_{i}-\bar{z}}{z_{\max}}\right)^{2}}$
(10)

The chosen RBF is a Gaussian $\phi$ of the voxel distances:

$(\mathbf{y_{L}})_{i}=\phi(i)=\exp\left[-(d_{i}/\sigma_{v})^{2}\right]$ (11)

a)  b)

Figure 4: a) $E(I)$ above 90% of Figure 3a, convolved by Gaussian kernel with
$\sigma=5$. b) The RBF is fitted to this region with center $\mu_{v}$ and
radius $\propto\sigma_{v}$.

Both outputs $\mathbf{y_{S}}$ and $\mathbf{y_{L}}$ are illustrated in Fig. 4.
According to the framework depicted so far, we can design a focus for many
different objects in the images by changing the functions for image feature
extraction, motion estimation, center, scale and RBF. It is also possible to
detect multiple objects or objects of complex shapes by a mixture of RBF
models.

###  3.5 Crop, scaling and fine CNN segmentation

This step performs the final preparation of the image so that it will fit the
CNN input shape; this is necessary because the MRI images may have diverse
dimensions. After all the RBF estimation steps executed in the original
resolution, the proposed region is cropped then adjusted, followed by an
intensity normalization. For CNNs with fixed input dimensions, we rescale
images using bicubic spline interpolation. For CNNs with variable input shape,
we only adjust the image proportions as requested by the model, e.g. U-net has
five max-pooling layers with down-scaling factor 2, which means the input
dimensions should be multiples of 25 in order to have at least 1 pixel in the
encoder output.

###  3.6 Neural Network ROI proposal

VMF was developed in combination with a 2D U-net CNN ronneberger2015u  for
simplicity, but it is suited to any other convolutional network. The U-net is
a general encoder-decoder segmentation architecture with success in LV
segmentation, as explained in Section 2.6. The CNN was trained with outputs
obtained with VMF during 30 epochs using an adaptive momentum optimizer,
initial learning rate $\eta$ = 0.001, Nesterov $\beta_{1}$ = 0.9,
$L_{\infty}$-decay $\beta_{2}$ = 0.999 and the loss below, which is binary
cross-entropy plus DSC score:

$\mathcal{L}(y,p)=-[y\cdot\log(p)+(1-y)\cdot\log(1-p)]+DSC(y,p)$ (12)

where DSC is defined by Equation 14 in Section 4.1.

The VMF pipeline executes in the whole cardiac volume and, after cropping, the
cropped ROI is passed to the segmentation CNN, which in this case is the 2D
U-net.

###  3.7 Datasets

To evaluate the VMF framework we used three CMR datasets specified across the
short-axis orientation. The metadata in the datasets include: a) binary masks
for the LV and RV; b) physiological parameters such as myocardium mass,
ventricle area, volume, ejection fraction, thickness, and dimensions of
structures; c) image acquisition parameters such as spatial resolution (mm),
temporal resolution (frames per cardiac cycle) and slice gaps (mm). Not all of
the datasets encompass the same information; following we provide more
details, with a summary presented in Table 1.

  * •

LVSC \- Cardiac Atlas Project (CAP) 2011 LV Segmentation Challenge
SUINESIAPUTRA2014 . We used the 100 patients in the training set, which have
LV masks on all the frames;

  * •

ACDC \- MICCAI 2017 Automated Cardiac Diagnosis Challenge Bernard2018 ,
training dataset. It was created from clinical data, including sequences of
100 patients from the University Hospital of Dijon (France) over 6 years, it
includes RV and LV masks in two frames;

  * •

M&Ms \- MICCAI 2020 Multi-Centre, Multi-Vendor & Multi-Disease Cardiac Image
Segmentation Challenge saber2020multi . This database was collected from six
hospitals in Spain, Canada, and Germany using several MRI scanners (Siemens,
GE, Philips and Canon), it includes RV and LV masks in two frames. We used 320
patients for which the labels are publicly available.

Dataset | LVSC | ACDC | M&Ms
---|---|---|---
n | 100 | 100 | 320
width | 138-512 | 154-428 | 196-548
height | 138-512 | 154-512 | 192-512
slices | 8-24 | 6-18 | 6-20
frames | 18-35 | 12-35 | 18-36
Sxy | 0.68-2.14 | 0.70-1.92 | 0.68-1.82
Sz | 6-10 | 5-10 | 5-10
Table 1: Overview of the three CMR datasets: n = number of patients;
(width,height,slices) = spatial image dimensions in (x,y,z) pixels; frames =
number of images in the time coordinate; Sxy = spatial resolution (pixel
spacing) in the axial plane (mm/pixel); Sz = slice resolution (mm/pixel).

##  4 Experiments

In order to evaluate VMF, we analyzed the selected datasets by applying VMF,
then compared the results of two CNNs: one base CNN on the raw images without
VMF, and another CNN combined with the images processed by VMF. After the
execution of VMF and cropping of ROIs, only the slices and frames with label
masks of at least 25 pixels are passed to the U-net (i.e. stray markings are
discarded). The labels are obtained by subtracting epicardium and endocardium
masks, i.e. they contain the myocardium region. The datasets were used
individually for training and validation of models, then combined in an all-
versus-all scheme.

###  4.1 Metrics

We compared three indices: (1) Recall, the proportion of the labels that was
preserved in the VMF region proposals; it is also known as Sensitivity or True
Positive Rate – TPR as defined by Equation 13, which aims to verify if the
bounding boxes cover the labels entirely; (2) the Sørensen-Dice Coefficient –
DSC as defined by Equation 14, which is equivalent to the F1-score (average of
Precision and Recall) that refers to the segmentation output when not-using vs
using method VMF; (3) speedup, the ratio of the times taken by the
segmentation CNN when not-using vs using method VMF – it is defined as
($t_{\mathrm{base}}/t_{\mathrm{ours}}$).

$TPR(y,\hat{y})=\frac{|y\cap\hat{y}|}{|y|}=\frac{TP}{TP+FN}$ (13)
$DSC(y,\hat{y})=\frac{2|y\cap\hat{y}|}{|y|+|\hat{y}|}=\frac{2TP}{2TP+FP+FN}$
(14)

Recall ranges from 0.0, when none of the marked voxels are detected; to 1.0,
when 100% of the marked voxels are detected. DSC ranges from 0.0, for no
intersection; to 1.0, for a perfect match between $y$ and $\hat{y}$. Speedup
is a positive ratio with value 1.0 when both methods take the same time;
$>$1.0 when VMF is faster; and $<$1.0 when VMF is slower.

###  4.2 Test Setup

The experiments ran in a rack-mounted machine with Intel i7-7700k CPU and
NVIDIA Titan-Xp GPU. Our scripts were written in Python, scipy,
keras/tensorflow, and matplotlib. The steps were: load an image; compute
features; fit RBF; crop and resize both the input image and the label mask;
Base CNN prediction on the original image and VMF-CNN prediction on the
cropped image; then, evaluate the predicted masks against the labels using the
aforementioned metrics.

##  5 Results

Table 2 presents the results of VMF region proposals (Recall) and VMF-CNN
combination (Dice score). After applying method VMF, the U-net CNN trained on
dataset ACDC had the worst Dice scores, with a mean performance decrease of
8.3% when considering all the datasets. On the other hand, the model trained
on dataset M&Ms was significantly improved when using the VMF ROI proposal,
with a mean performance increase of +7.2 (percent points) considering all the
datasets. With dataset LVSC, we observed a lower performance when validating
over dataset LVSC (-3.5), but an improvement when validating with datasets
ACDC and M&Ms (+4.0 and +2.5, respectively).

Train Set | Test | Base | VMF | Dice
---|---|---|---|---
(Recall) | Set | Dice | Dice | Diff.
| LVSC | 50.1% | 44.8% | -5.7
ACDC | ACDC | 74.6% | 69.9% | -4.7
Recall of 98.84% | M&Ms | 67.6% | 47.4% | -20.2
| all | 56.6% | 48.3% | -8.3
| LVSC | 62.4% | 70.9% | +8.5
M&Ms | ACDC | 60.8% | 67.0% | +6.2
Recall of 99.75% | M&Ms | 81.3% | 82.2% | +0.9
| all | 65.4% | 72.6% | +7.2
| LVSC | 74.3% | 70.8% | -3.5
LVSC | ACDC | 75.7% | 79.7% | +4.0
Recall of 99.75% | M&Ms | 67.0% | 69.5% | +2.5
| all | 74.5% | 72.7% | -1.8
| LVSC | 76.5% | 77.5% | +1.0
all | ACDC | 85.9% | 86.8% | +0.9
Recall of 99.69% | M&Ms | 77.2% | 83.0% | +5.8
| all | 77.7% | 79.4% | +1.7
Table 2: Performance results. In the first (left-most) column, the name of the
dataset and the Recall metric (Section 4.1. In the remaining columns, the Dice
scores (metric DSC in Section 4.1) for CNN segmentation without VMF (Base) and
combined with VMF.

We applied the McNemar’s chi-squared test with continuity correction from
package statsmodels 0.12.2, whose results asserted that the Base predicitons
and predictions after VMF were significantly different, with $p<.001$ (see
Table 3).

|  VMF (T) |  VMF (F)
---|---|---
Base (T) | 6,041,403 | 2,543,222
Base (F) | 818,157 | 709,767,778
| $\chi^{2}=885305.1$ | $(p<10^{-302})$
Table 3: Contingency matrix considering all voxels in the configuration
train:all versus test:all, refer to Table 2; and the corrected McNemar’s chi-
squared test comparing CNN without VMF (Base) versus with VMF asserted
significantly different predictions, with $p<.001$.

From a practical perspective we also compared the training speed, which is
paramount to larger experiments such as hyper-parameter search. Table 4 shows
the training time of the networks without and with method VMF; in all the
cases our method improved the training speed by 150%.

Train Set | Base | VMF | Speed-up
---|---|---|---
ACDC | 21.0s | 11.6s | 1.81x
M&Ms | 93.2s | 37.6s | 2.48x
LVSC | 254.7s | 95.1s | 2.68x
all | 427.2s | 170.3s | 2.51x
Table 4: Mean training time of 30 epochs for CNN segmentation without VMF
(Base) and after VMF. In all the cases, method VMF accelerates the training
speed; this is because the CNN processes only the data in the region of
interest.

##  6 Discussion

Based on our results, the main achievement of method VMF is the ability to
identify and crop the ROIs with a very small error (Recall = 99.69%).
Furthermore, by cropping the ROIs before the CNN training, we observed that
the training process was remarkably accelerated, increasing the training speed
by 150%. When considering metric DSC, we observed significant improvements
when considering the entire bundle of experiments, as presented in Table 2
but, for dataset ACDC when used for training, the results were not promising.

This inefficiency with ACDC is possibly related to imaging artifacts and
standard operating procedures used to label the ACDC images (refer to Figure
5). In general, VMF is able to automatically focus on the correct region and
guide the segmentation CNN (Figure 5, P27). However, in some cases the CNN
lost track of the LV when fed with a slightly displaced ROI (P32). In the
rarest of cases the CNN is unable to detect anything, even when fed with a
centered ROI (P91), and this is a known problem when segmenting the cardiac
apex. These findings should be investigated in future work.

Figure 5: Examples of slices with DSC$<$0.2 found in ACDC patients number 27,
32 and 91 (rows). Left column shows the original label for each patient. The
center column shows the prediction by the Base CNN. The right column shows
both the output of VMF (bright boxes) which is used for cropping, and the
prediction of the CNN after VMF.

Most of these issues were subsequently reviewed in M&Ms saber2020multi  by: a)
complete coverage, including papillary muscles, of RV and LV; b) no
interpolation of masks at the base; c) guarantees of larger RV surface in end-
diastole compared to end-systole, excluding the pulmonary artery; d) sampling
data from several hospitals, different patients and different equipment
vendors. The outcome of these corrections is clear when comparing the results
from training in ACDC and testing M&Ms (-20.2) vs training on M&Ms and testing
on ACDC (+6.2). M&Ms is very likely to contain images following the same
patterns of ACDC and others, but not the opposite.

Overall, method VMF demonstrated significant improvements in the task of ROI
detection; from the results, it became evident that the performance depends on
the training dataset and on the network model – image quality, format and
intensities should be similar across datasets, and labelling standard
operating procedures should be compatible.

##  7 Conclusions

In this paper, we proposed method VMF, a novel approach based on convolution
operations and on the use of a Radial Basis Function to detect the ROIs in
cardiac magnetic resonance images. We validated VMF with a U-net CNN comparing
our results to those of the canonical U-net and of the VMF-CNN in three public
reference datasets. According to our results, VMF was able to crop 99.69%
(Recall metric) of the ROI voxels in all the datasets, being suitable to
preprocess the data for CNN segmentation. VMF accelerated the training process
by 150%, and also increased metric Sørensen-Dice Coefficient in the majority
of our test cases ($p<.001$). Further improvement is possible by extended
preprocessing of the training datasets, and change the U-net CNN after the ROI
proposal step to a more advanced architecture with fine-tuning to the dataset
specifics, specially in cases such as ACDC.

As future work, we intend to expand this pipeline considering possibilities:
(1) embed this method as the first layer of a CNN, so that the parameters are
adjusted automatically; (2) use other networks in the segmentation step, such
as 3D Unet or Feature Pyramid Networks; (3) experiment with more datasets,
preprocessing, and augmentation methods; (4) test different feature extractors
for the static visual features, such as Gray-Level Co-occurrence Matrices and
Gabor filters, other motion estimation methods, and Radial Basis Functions.

## Acknowledgment

We thank Dr. Diego Cardenas and Dr. Thiago Costa for reviewing equations,
Nvidia Corporation for donating the GPU used in this work, and USP, FAPESP,
CAPES, CNPq, Zerbini Foundation and FoxConn for supporting this research.

## Declarations

Conflict of interest: The authors declare that they have no conflict of
interest.

Research data: The datasets used in this research are available upon
registration on https://www.ub.edu/mnms/, https://acdc.creatis.insa-
lyon.fr/description/databases.html and
https://www.cardiacatlas.org/challenges/lv-segmentation-challenge/

## References

  * [1] O. Bernard, A. Lalande, C. Zotti, F. Cervenansky, X. Yang, P. A. Heng, I. Cetin, K. Lekadir, O. Camara, M. A. Gonzalez Ballester, G. Sanroma, S. Napel, S. Petersen, G. Tziritas, E. Grinias, M. Khened, V. A. Kollerathu, G. Krishnamurthi, M. M. Rohé, X. Pennec, M. Sermesant, F. Isensee, P. Jäger, K. H. Maier-Hein, P. M. Full, I. Wolf, S. Engelhardt, C. F. Baumgartner, L. M. Koch, J. M. Wolterink, I. Išgum, Y. Jang, Y. Hong, J. Patravali, S. Jain, O. Humbert, P. M. Jodoin, Deep learning techniques for automatic mri cardiac multi-structures segmentation and diagnosis: Is the problem solved?, IEEE Transactions on Medical Imaging 37 (11) (2018) 2514–2525.  doi:10.1109/TMI.2018.2837502.
  * [2] W. Xue, G. Brahm, S. Pandey, S. Leung, S. Li, Full left ventricle quantification via deep multitask relationships learning, Medical image analysis 43 (2018) 54–65.
  * [3] M. Khened, V. Alex, G. Krishnamurthi, Densely connected fully convolutional network for short-axis cardiac cine mr image segmentation and heart diagnosis using random forest, in: International Workshop on Statistical Atlases and Computational Models of the Heart, Springer, 2017, pp. 140–151.
  * [4] P. Peng, K. Lekadir, A. Gooya, L. Shao, S. E. Petersen, A. F. Frangi, A review of heart chamber segmentation for structural and functional analysis using cardiac magnetic resonance imaging, Magn Reson Matter Phy 19 (2016) 155–195.
  * [5] K. He, X. Zhang, S. Ren, J. Sun, Spatial pyramid pooling in deep convolutional networks for visual recognition, IEEE transactions on pattern analysis and machine intelligence 37 (9) (2015) 1904–1916.
  * [6] B. Wu, Y. Fang, X. Lai, Left ventricle automatic segmentation in cardiac mri using a combined cnn and u-net approach, Computerized Medical Imaging and Graphics 82 (2020) 101719.  doi:https://doi.org/10.1016/j.compmedimag.2020.101719.
URL https://www.sciencedirect.com/science/article/pii/S0895611120300227

  * [7] W. Xue, A. Lum, A. Mercado, M. Landis, J. Warrington, S. Li, Full quantification of left ventricle via deep multitask learning network respecting intra- and inter-task relatedness, in: M. Descoteaux, L. Maier-Hein, A. Franz, P. Jannin, D. L. Collins, S. Duchesne (Eds.), Medical Image Computing and Computer Assisted Intervention - MICCAI 2017, Springer International Publishing, Cham, 2017, pp. 276–284.
  * [8] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, A. Rabinovich, Going deeper with convolutions, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 1–9.
  * [9] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, in: Advances in neural information processing systems, 2017, pp. 5998–6008.
  * [10] J. P. Earls, V. B. Ho, T. K. Foo, E. Castillo, S. D. Flamm, Cardiac mri: recent progress and continued challenges, Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine 16 (2) (2002) 111–127.
  * [11] A. Seraphim, K. D. Knott, J. Augusto, A. N. Bhuva, C. Manisty, J. C. Moon, Quantitative cardiac mri, Journal of Magnetic Resonance Imaging 51 (3) (2020) 693–711.
  * [12] P. A. Iaizzo, R. H. Anderson, A. J. Hill, The importance of human cardiac anatomy for translational research, Journal of cardiovascular translational research 6 (2) (2013) 105–106.
  * [13] M. Saber, D. Abdelrauof, M. Elattar, Multi-center, multi-vendor, and multi-disease cardiac image segmentation using scale-independent multi-gate unet, in: International Workshop on Statistical Atlases and Computational Models of the Heart, Springer, 2020, pp. 259–268.
  * [14] M. Spann, R. Wilson, A quad-tree approach to image segmentation which combines statistical and spatial information, Pattern Recognition 18 (3-4) (1985) 257–269.
  * [15] J. Liu, Y.-H. Yang, Multiresolution color image segmentation, IEEE Transactions on Pattern Analysis and Machine Intelligence 16 (7) (1994) 689–700.
  * [16] R. Girshick, J. Donahue, T. Darrell, J. Malik, Rich feature hierarchies for accurate object detection and semantic segmentation, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2014, pp. 580–587.
  * [17] J. Redmon, S. Divvala, R. Girshick, A. Farhadi, You only look once: Unified, real-time object detection, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 779–788.
  * [18] S. Ranganath, Contour extraction from cardiac mri studies using snakes, IEEE transactions on medical imaging 14 (2) (1995) 328–338.
  * [19] G. Mühlenbruch, M. Das, C. Hohl, J. E. Wildberger, D. Rinck, T. G. Flohr, R. Koos, C. Knackstedt, R. W. Günther, A. H. Mahnken, Global left ventricular function in cardiac ct. evaluation of an automated 3d region-growing segmentation algorithm, European radiology 16 (5) (2006) 1117–1123.
  * [20] K. He, G. Gkioxari, P. Dollár, R. Girshick, Mask r-cnn, in: Proceedings of the IEEE international conference on computer vision, 2017, pp. 2961–2969.
  * [21] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, C. L. Zitnick, Microsoft coco: Common objects in context, in: European conference on computer vision, Springer, 2014, pp. 740–755.
  * [22] R. C. Gonzalez, R. E. Woods, et al., Digital image processing (2002).
  * [23] J. C. de Oliveira, Motion Estimation, Springer US, Boston, MA, 2006, pp. 426–434.
URL https://doi.org/10.1007/0-387-30038-4_131

  * [24] B. K. Horn, B. G. Schunck, Determining optical flow, Artificial intelligence 17 (1-3) (1981) 185–203.
  * [25] D. Thanou, P. A. Chou, P. Frossard, Graph-based motion estimation and compensation for dynamic 3d point cloud compression, in: 2015 IEEE International Conference on Image Processing (ICIP), IEEE, 2015, pp. 3235–3239.
  * [26] Y. Zhang, C. Zhang, R. Fan, Fast motion estimation in hevc inter coding: An overview of recent advances, in: 2018 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), 2018, pp. 49–56.  doi:10.23919/APSIPA.2018.8659650.
  * [27] W. Bao, W.-S. Lai, X. Zhang, Z. Gao, M.-H. Yang, Memc-net: Motion estimation and motion compensation driven neural network for video interpolation and enhancement, IEEE transactions on pattern analysis and machine intelligence (2019).
  * [28] B. Stemkens, R. H. Tijssen, B. D. De Senneville, J. J. Lagendijk, C. A. Van Den Berg, Image-driven, model-based 3d abdominal motion estimation for mr-guided radiotherapy, Physics in Medicine & Biology 61 (14) (2016) 5335.
  * [29] Y. Q. Mohsin, S. G. Lingala, E. DiBella, M. Jacob, Accelerated dynamic mri using patch regularization for implicit motion compensation, Magnetic resonance in medicine 77 (3) (2017) 1238–1248.
  * [30] W. Yan, Y. Wang, Z. Li, R. J. Van Der Geest, Q. Tao, Left ventricle segmentation via optical-flow-net from short-axis cine mri: preserving the temporal coherence of cardiac motion, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2018, pp. 613–621.
  * [31] H. Yu, S. Sun, H. Yu, X. Chen, H. Shi, T. S. Huang, T. Chen, Foal: Fast online adaptive learning for cardiac motion estimation, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 4313–4323.
  * [32] Z. Majdisova, V. Skala, Radial basis function approximations: comparison and applications, Applied Mathematical Modelling 51 (2017) 728–743.
  * [33] D. S. Broomhead, D. Lowe, Radial basis functions, multi-variable functional interpolation and adaptive networks, Tech. rep., Royal Signals and Radar Establishment Malvern (United Kingdom) (1988).
  * [34] M. D. Buhmann, Radial basis functions: theory and implementations, Vol. 12, Cambridge university press, 2003.
  * [35] O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, in: International Conference on Medical image computing and computer-assisted intervention, Springer, 2015, pp. 234–241.
  * [36] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, S. Belongie, Feature pyramid networks for object detection, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2017, pp. 2117–2125.
  * [37] A. Chaurasia, E. Culurciello, Linknet: Exploiting encoder representations for efficient semantic segmentation, in: 2017 IEEE Visual Communications and Image Processing (VCIP), IEEE, 2017, pp. 1–4.
  * [38] E. Smistad, A. Østvik, B. O. Haugen, L. Løvstakken, 2d left ventricle segmentation using deep learning, in: 2017 IEEE International Ultrasonics Symposium (IUS), 2017, pp. 1–4.  doi:10.1109/ULTSYM.2017.8092573.
  * [39] F. Guo, M. Ng, G. Wright, Cardiac mri left ventricle segmentation and quantification: a framework combining u-net and continuous max-flow, in: International Workshop on Statistical Atlases and Computational Models of the Heart, Springer, 2018, pp. 450–458.
  * [40] S. M. Fallat, C. R. Johnson, Hadamard powers and totally positive matrices, Linear Algebra and its Applications 423 (2-3) (2007) 420–427.
  * [41] A. Suinesiaputra, B. R. Cowan, A. O. Al-Agamy, M. A. Elattar, N. Ayache, A. S. Fahmy, A. M. Khalifa, P. Medrano-Gracia, M.-P. Jolly, A. H. Kadish, D. C. Lee, J. Margeta, S. K. Warfield, A. A. Young, A collaborative resource to build consensus for automated left ventricular segmentation of cardiac mr images, Medical Image Analysis 18 (1) (2014) 50–62.  doi:https://doi.org/10.1016/j.media.2013.09.001.
URL https://www.sciencedirect.com/science/article/pii/S1361841513001217


11institutetext: Department of Radiology, Lausanne University Hospital (CHUV)
and University of Lausanne (UNIL), Lausanne, Switzerland 22institutetext: CIBM
Center for Biomedical Imaging, Switzerland 33institutetext: BCNatal Fetal
Medicine Research Center (Hospital Clínic and Hospital Sant Joan de Déu),
Universitat de Barcelona, Spain 44institutetext: IDIBAPS and CIBERER,
Barcelona, Spain 55institutetext: Department Woman-Mother-Child, CHUV,
Lausanne, Switzerland

# FetMRQC: Automated Quality Control
for fetal brain MRI††thanks: TS is supported by the Era-net Neuron MULTIFACT
project (SNSF 31NE30_203977), OE is supported by the Swiss National Science
Foundation (SNSF #185872), the NIMH (RF1MH12186), and the CZI
(EOSS5/‘NiPreps’). YG acknowledges support from the SICPA foundation and EE is
supported by the Instituto de Salud Carlos III (ISCIII) (AC21_2/00016). We
acknowledge access to the facilities and expertise of the CIBM Center for
Biomedical Imaging, a Swiss research center of excellence founded and
supported by CHUV, UNIL, EPFL, UNIGE and HUG.

Thomas Sanchez 1122 Oscar Esteban 11 Yvan Gomez 3355 Elisenda Eixarch 3344
Meritxell Bach Cuadra 2211

###### Abstract

Quality control (QC) has long been considered essential to guarantee the
reliability of neuroimaging studies. It is particularly important for fetal
brain MRI, where large and unpredictable fetal motion can lead to substantial
artifacts in the acquired images. Existing methods for fetal brain quality
assessment operate at the slice level, and fail to get a comprehensive picture
of the quality of an image, that can only be achieved by looking at the entire
brain volume. In this work, we propose FetMRQC, a machine learning framework
for automated image quality assessment tailored to fetal brain MRI, which
extracts an ensemble of quality metrics that are then used to predict experts’
ratings. Based on the manual ratings of more than 1000 low-resolution stacks
acquired across two different institutions, we show that, compared with
existing quality metrics, FetMRQC is able to generalize out-of-domain, while
being interpretable and data efficient. We also release a novel manual quality
rating tool designed to facilitate and optimize quality rating of fetal brain
images.

Our tool, along with all the code to generate, train and evaluate the model
will be released upon acceptance of the paper.

###### Keywords:

Image quality assessment Fetal brain MRI

##  1 Introduction

Magnetic Resonance Imaging (MRI) of the fetal brain is increasingly
complementing ultrasound imaging to diagnose abnormalities, thanks to its
unmatched soft tissue contrast [1, 2]. However, inherent noise sources and
imaging artifacts, such as fetal motion, can degrade the quality of the
acquired images and jeopardize subsequent analyses. Indeed, insufficient MRI
data quality has been shown to bias neuroradiological assessment and
statistical analyses [3, 4, 5]. Therefore, establishing objective image
quality assessment and control (QA/QC) protocols for neuroimaging studies has
long been considered critical to ensure their reliability, generalization, and
replicability [6, 7]. More recently, some QA/QC tools have been popularized
for the adult brain and attempt to automate and aid exclusion decisions [8, 9,
10, 11, 12]. Unfortunately, these techniques have found generalization across
imaging devices extremely challenging [8]. Moreover, they are inapplicable to
fetal MRI, as they rely on invalid prior knowledge, e.g., assuming that the
head is surrounded by air or the relative position of the brain with respect
to the axis of the bore of the scanner.

In fetal brain MRI, QA/QC has been approached implicitly within super-
resolution reconstruction (SRR) methods [13, 14, 15, 16, 17]. SRR builds a
high-resolution, isotropic 3D volume from several differently-oriented,
consecutive stacks of 2D slices with substantially lower resolution on the
through-plane axis (i.e., anisotropic resolution). Typically, SRR methods rely
on a manual selection of stacks that are deemed adequate for reconstruction.
These methods may incorporate an automated QC stage for outlier rejection that
excludes sub-standard slices or pixels from the input low-resolution stacks
[13, 14, 18, 17]. Outside of fetal brain MRI, Uus et al. [19] have explored
methods to automate the rejection of poorly registered stacks. However, bad
quality series can remain detrimental to the final quality of the
reconstruction, even when SRR pipelines include outlier rejection schemes.

Figure 1: Examples of clinical image quality: left, excellent quality, and
right, low quality. On the right, the stack would be excluded from further
analysis due to significant intensity changes between many slices and strong
signal drop; in the through-plane view, strong inter-slice motion makes it
difficult to discern the structure of the brain.

Recent work has applied deep learning on QA/QC of fetal brain MRI [20, 21,
22]. These solutions automatically identify problematic slices for exclusion
(QC), and, if streamlined with the acquisition, their short inference time
permits re-acquiring corrupted slices [23] (QA). However, not all artifacts
can be seen by analyzing single slices. For instance, inter-slice motion
(visible on the right of Figure 1), a strong bias field in the through-plane
direction, or an incomplete field of view can often clearly be seen only when
considering the entire volume.

We propose FetMRQC, a framework for QA/QC of clinical fetal brain MRI,
operating on stacks of LR, T2-weighted (T2w) images. Inspired by MRIQC [8],
FetMRQC generates a visual report corresponding to each input stack, for
efficient screening and manual QA. The tool also automatically extracts an
ensemble of image quality metrics (IQMs) that reflect some quality feature.
Finally, we propose a learning framework to automatically predict image
quality from the IQMs corresponding to new images. We manually assessed more
than 1000 LR stacks acquired across two different institutions and several MRI
scanners within each. Using these manual quality annotations, we
quantitatively assess FetMRQC in two QA/QC tasks, namely, regression and
binary classification (inclusion/exclusion). Our results demonstrate the
feasibility of automated QA/QC of T2w images of the fetal brain and the
substantial improvement of SRR after QC of subpar stacks.

##  2 Methodology

FetMRQC comprehends two major elements to implement QA/QC protocols of
unprocessed (LR stacks) fetal brain MRI data. First, the tool builds upon
MRIQC’s paradigm and generates an individual QA report for each LR stack to
assist and optimize screening and annotation by experts. Second, FetMRQC
proposes a learning framework (Figure 3) based on image quality metrics (IQMs)
extracted from the data to automate the assessment.

###  2.1 Data

We retrospectively collected LR T2w stacks from 146 subjects retrieved from
two existing databases at two different institutions, SITE1 and SITE2. The
data were acquired on different Siemens scanners and a common MR scheme (Half
fourier Single-shot Turbo spin-Echo; HASTE), at 1.5T and 3T, with both normal
and pathological cases. The full list of parameters is detailed in Table 4 in
the Supplementary material. The corresponding local ethics committees
independently approved the studies under which data were collected, and all
participants gave written informed consent. SITE1 provided 61 subjects, with
an average of 7.9$\pm$3.0 LR stacks per subject. The SITE2 provided 85
subjects, $5.8\pm 3.4$ stacks per subject. The aggregate sample size is
$N$=1010 LR series.

###  2.2 Manual QA of fetal MRI stacks

Akin to MRIQC [8], FetMRQC generates an HTML-based report adapted to the QA of
fetal brains for each input LR stack (Figure 2). The input dataset is required
to comply with the Brain Imaging Data Structure (BIDS)[24], a format widely
adopted in the neuroimaging community. The reports are generated using an
image with a corresponding brain mask. This mask can be extracted
automatically, and in this work, we used MONAIfbs [25]. Each individual-stack
report has a QA utility (the so-called rating widget, with which raters can
fill in an overall quality score, the in-plane orientation, and the presence
and grading of artifacts showcased by the stack. We chose also to use an
interval (as opposed to categorical) rating scale with four main quality
ranges: $[0,1]$: exclude – $(1,2]:$ poor – $(2,3]$: acceptable – $(3,4]$:
excellent. Interval ratings simplify statistical modeling, set lower bounds to
annotation noise, and enable the inference task where a continuous quality
score is assigned to input images, rather than broad categories.

We generated the individual report corresponding to each of the 1010 LR
stacks. Two raters independently annotated each 555 of the dataset using the
proposed tool, which yielded 100 stacks annotated by both raters to test
inter-rater variability. These 100 stacks were randomly selected. Rater1 is a
maternal-fetal physician with 5 years of experience in neuroimaging, and
Rater2 is an engineer with 20 years of experience in neuroimaging.

Figure 2: Visual quality assessment tool proposed in this work. The HTML
report gives some general information on the anonymized scan, displays all
slices with brain content, as well as both orthogonal through-plane views, as
seen on Figure 1.

###  2.3 IQMs extraction and learning

Figure 3: Proposed automated IQA pipeline. Given an input image and the
corresponding brain mask, an ensemble of IQMs are extracted, and then
processed in a prediction model. IQMs are normalized and irrelevant entries
are removed before being fed into a regression or classification model, that
produces the final outcome.

####  2.3.1 IQMs in fetal brain MRI.

Leveraging the same workflow that generates the individual reports, a number
of IQMs are extracted from the full stack and from the information within the
brain mask. Since only a few IQMs defined by MRIQC can be applied to fetal
brains, we implemented a set of IQMs specific to the application at hand. In
[14], Kainz et al. proposed a rank_error based on the estimated
compressibility of an LR stack. Later on, brain mask-derived measures were
used such as the volume of the brain mask, mask_volume, in [18] or its
centroid in [26]. Recently, deep learning-based IQM slice-wise and stack-wise
have been proposed, dl_slice_iqa[21] and dl_stack_iqa [27], which respectively
do slice-wise and stack wise quality assessment111We use the pre-trained
models throughout these experiments, as we want to test the off-the-shelf
value of these IQMs222The method of Liao et al. [22] was not included because
their code is not publicly available, and we could not get in contact with the
authors..

We propose additional IQMs for quality prediction that have not previously
been used in the context of fetal brain MRI. They can roughly be categorized
into two groups, and full details are provided in Table 4 in the Supplementary
material. In a nutshell, intensity-based IQMs directly rely on the voxel
values of the image. These include summary statistics [8] like mean, median,
and percentiles. We also re-use metrics traditionally used for outlier
rejection [13, 14, 18] to quantify the intensity difference between slices in
a volume. We also compute entropy [8], estimate the level of bias using N4
bias field correction [28] and estimate the sharpness of the image with
Laplace and Sobel filters. The second type of metrics are shape-based and
operate directly on the automatically extracted brain mask. We propose to use
a morphological closing to detect through-plane motion, as well as edge
detection to estimate the variation of the surface of the brain mask.

Variation on the metrics. All the IQMs operate on LR images or masks and can
be modified by various transformations of the data. For instance, Kainz et al.
[14] evaluated their metrics only on the third of the slices closest to the
center of a given LR volume. We include various variations on the IQMs,
including considering the whole ROI instead of the centermost slices (_full).
Other variations include keeping or masking the maternal background,
aggregating point estimates using mean, median, or other estimators, computing
information theoretic metrics on the union or intersection of masks, etc.
Finally, slice_loss metrics can be either computed as a pairwise comparison
between all slices (by default) or only on a window of neighboring slices
(_window). With all the different variations, we arrive at a total of $75$
different features333See Figure 8 in the Supplementary material for a cross-
correlation matrix between all metrics on our entire dataset..

Feature-wise normalization. We explored both global and subject-wise
normalization. In our experiments, we did not notice substantial changes of
using one over another, and for the sake of simplicity, we adopted subject-
wise standard scaling in the rest of this work.

Feature-selection and dimensionality reduction. We explore feature variance
and correlations. We remove irrelevant features ($\text{variance}=0$ for a
given subset of data) and drop IQMs that are highly correlated with each other
(with threshold $0.8$ and $0.9$). We also remove features which do not
contribute more than noise using the Winnow algorithm [29] with extremely
randomized trees [8]. Finally, we also explore using principal component
analysis to construct orthogonal features.

Model selection and evaluation through nested cross-validation. Nested cross-
validation (CV) is a fully automated mechanism to perform model selection and
evaluation without introducing optimistic biases [30]. In this framework, an
outer CV loop benchmarks the expected performance of the family of models
represented by the inner loop. An inner CV loop is executed for each fold of
the outer loop to select the best-performing model on average. We set up the
nested CV framework with five folds in both the inner and outer loops ($80\%$
train, $20\%$ test at each level), for both the regression and classification
tasks We ensured to group together all the LR stacks of each subject to avoid
data leaking between the training and testing fold. For the regression task,
we evaluated linear regression, gradient boosting, and random forests. For
classification, we considered logistic regression, random forests, gradient
boosting, and AdaBoost. We primarily optimized feature selection strategies
and model fitting parameters. A detailed description of the combinations of
models and parameters optimized is available in Supplementary Table 5. The
experiments were implemented with Python 3.9.15 and Scikit-learn 1.1.3 [31].

Baseline models and evaluation. We evaluated three variants of our model.
First of all, we assessed the predictive power of a subset of individual
features that have previously been proposed and used for fetal brain QA/QC
(namely, rank_error, rank_error_full mask_volume, centroid,
centroid_full,dl_slice_iqa and dl_stack_iqa). The individual features were
scaled and then fitted with a linear regression or logistic regression model
for prediction. We then reported the performance of the best-performing
feature, which consistently was the dl_stack_iqa [27]. Secondly, we considered
a base variant of our IQA pipeline, using the same subset of features above.
Finally, we considered a model using all available features, referred to as
FetMRQC. Our regression results are evaluated using mean absolute error (MAE)
and Spearman rank correlation. Our classification results use F1-score and the
area under the receiver operating characteristic curve (AUC ROC).

##  3 Results and discussion

####  3.0.1 FetMRQC’s reports optimize stack screening and bolster inter-rater
reliability.

Figure 6 summarizes the high agreement between the two raters, with a Pearson
correlation coefficient of $0.79$ and $0.83$ for SITE1 and SITE2. On SITE1,
120 series are manually rated below the exclude threshold $\text{Quality}<1$),
and 414 above. On SITE2, 151 series are rated excluded, and 393 rated above.
The total include-to-exclude ratio is 2.98. The total rating time was around
5h 40min for Rater 1 (median of 37s per volume), and around 6h10 for Rater 2
(median of 40s per volume).

Figure 4: Inter-rater correlation: blue and yellow data are respectively from
SITE1 and SITE2.

Figure 5: MAE as a function of the size of the training set. The individual
metric, dl_stack_iqa_full is constant because it only performs inference.

Table 1: In-domain evaluation

| Regression | Classification
---|---|---
| MAE ($\downarrow$) | Spearman ($\uparrow$) | F1 ($\uparrow$) | AUC ($\uparrow$)
dl_stack_iqa | $0.72\pm 0.05$ | $0.37\pm 0.11$ | $0.85\pm 0.02$ | $0.53\pm 0.02$
Base | $0.58\pm 0.05$ | $0.60\pm 0.03$ | $0.88\pm 0.02$ | $0.71\pm 0.06$
FetMRQC | $\mathbf{0.53\pm 0.09}$ | $\mathbf{0.71\pm 0.05}$ | $\mathbf{0.90\pm 0.02}$ | $\mathbf{0.77\pm 0.07}$

Table 2: Out of domain evaluation

| Regression | Classification
---|---|---
| MAE ($\downarrow$) | Spearman ($\uparrow$) | F1 ($\uparrow$) | AUC ($\uparrow$)
dl_stack_iqa | $0.75\pm 0.03$ | $0.42\pm 0.10$ | $0.83\pm 0.04$ | $0.57\pm 0.07$
Base | $0.75\pm 0.01$ | $0.38\pm 0.03$ | $0.85\pm 0.02$ | $0.67\pm 0.06$
FetMRQC | $\mathbf{0.68\pm 0.06}$ | $\mathbf{0.50\pm 0.08}$ | $\mathbf{0.89\pm 0.01}$ | $\mathbf{0.77\pm 0.05}$

Figure 6: Feature importance averaged on regression (left) and classification
models (right), for the $20$ largest contributors. Gray features are not
correlated with each other ($|\text{corr}|<0.75$). Features with the same
color are correlated ($>0.75$) with each other. Hatched features are features
that were proposed in this work.

FetMRQC’s new IQMs reach the best performance when evaluated within-domain.
Table 6 shows that the additional features substantially increase the
performance of the models, both in the case of regression and classification.
We used subjects from both sites on the training fold and different subjects
from both sides on the validation fold.

Using more IQMs allows to generalize with fewer training data within domain.
On Figure 6, we see that not only using more features reduces the overall
error on the validation set, but allows to reach a good testing performance
using as little as $30\%$ ($\sim 40$ subjects, $300$ LR series) of the data
for training.

FetMRQC’s classifier is robust when trained on one site and evaluated on the
other. On Table 6, we see that classification is retained when training on one
of the sites and evaluating on the other. Although regression performance is
hurt by the domain gap between sizes, our proposes model still generalizes
best. We used nested CV to tune the hyperparameters on the data from one site
and evaluated the model on the other site. In the inner loop, data were still
split by subject.

Different features matter for different tasks. An advantage of feature-based
methods over deep learning-based methods is that their decisions are based on
a combination of simple features, which make them more interpretable. On
Figure 6, we see a comparison of the feature importance for regression and
classification tasks. These were obtained as by averaging the feature
importance for the model selected at each fold of the nested CV. We observe
that the importance of features varies largely with the task at hand. While
the closing of the brain mask (closing_mask_full and the nMAE are the most
predictive of the regression quality), the exclusion of slices relies more on
the compressibility of the centermost slices of the image (rank_error_center)
and the edges of the brain mask (laplace_mask_full).

Figure 7: Illustration of usefulness of QC on the SR reconstruction, using
NiftyMIC [18], and NeSVoR [17]. FetMRQC+NiftyMIC successfully removes 6 series
out of the 13 available, and FetMRQC+NeSVoR removes 2 out of 5. FetMRQC leads
to a substantially greater image quality.

FetMRQC’s classifier can drastically improve the outcome of SRR methods.
Figure 7 shows how bad quality stacks can drastically impact the final quality
of the SRR, using NiftyMIC [18] and NeSVoR [17]. The images reconstructed
without QC displayed significant artifacts. Conversely, the two SRR methods
showcased improved quality when subpar stacks were excluded.

##  4 Conclusion

We propose FetMRQC, which adapts the approach of a popular tool, MRIQC, to
implement reliable QA/QC of low-resolution stacks of T2w MRI of the fetal
brain. Using FetMRQC tools, two experts assessed the quality of 1010 stacks
with high inter-rater reliability; and these annotations were then applied to
automated regression (for QA) and classification (QC). Nested cross-validation
of a set of models and hyper-parameters showed how QA and QC can be automated.
Objective and reliable QA/QC procedures are critical to ensure the reliability
and repeatability of neuroimaging studies, and FetMRQC demonstrates how
existing approaches can readily be applied on fetal brain MRI.

## References

  * [1] A. Gholipour et al., “Fetal MRI: a technical update with educational aspirations,” _Concepts in Magnetic Resonance Part A_ , vol. 43, no. 6, pp. 237–266, 2014.
  * [2] S. N. Saleem, “Fetal MRI: An approach to practice: A review,” _Journal of advanced research_ , vol. 5, no. 5, pp. 507–523, 2014.
  * [3] J. D. Power et al., “Spurious but systematic correlations in functional connectivity mri networks arise from subject motion,” _Neuroimage_ , vol. 59, no. 3, pp. 2142–2154, 2012.
  * [4] M. Reuter et al., “Head motion during mri acquisition reduces gray matter volume and thickness estimates,” _Neuroimage_ , vol. 107, pp. 107–115, 2015\.
  * [5] A. Alexander-Bloch et al., “Subtle in-scanner motion biases automated measurement of brain anatomy from in vivo mri,” _Human brain mapping_ , vol. 37, no. 7, pp. 2385–2397, 2016.
  * [6] B. Mortamet et al., “Automatic quality assessment in structural brain magnetic resonance imaging,” _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , vol. 62, no. 2, pp. 365–372, 2009.
  * [7] G. Niso et and R. Botvinik-Nezer et al., “Open and reproducible neuroimaging: from study inception to publication,” _NeuroImage_ , p. 119623, 2022.
  * [8] O. Esteban et al., “MRIQC: Advancing the automatic prediction of image quality in MRI from unseen sites,” _PloS one_ , vol. 12, no. 9, p. e0184661, 2017.
  * [9] E. T. Klapwijk et al., “Qoala-t: A supervised-learning tool for quality control of freesurfer segmented mri data,” _Neuroimage_ , vol. 189, pp. 116–129, 2019.
  * [10] C. Vogelbacher et al., “Lab–qa2go: a free, easy-to-use toolbox for the quality assessment of magnetic resonance imaging data,” _Frontiers in neuroscience_ , vol. 13, p. 688, 2019.
  * [11] Z. R. Samani et al., “Qc-automator: Deep learning-based automated quality control for diffusion mr images,” _Frontiers in neuroscience_ , vol. 13, p. 1456, 2020.
  * [12] M. Garcia et al., “BrainQCNet: a deep learning attention-based model for multi-scale detection of artifacts in brain structural mri scans,” _bioRxiv_ , 2022.
  * [13] M. Kuklisova-Murgasova et al., “Reconstruction of fetal brain MRI with intensity matching and complete outlier removal,” _Medical Image Analysis_ , vol. 16, no. 8, 2012.
  * [14] B. Kainz et al., “Fast volume reconstruction from motion corrupted stacks of 2D slices,” _IEEE Transactions on Medical Imaging_ , vol. 34, no. 9, 2015\.
  * [15] S. Tourbier et al., “Medical-image-analysis-laboratory/mialsuperresolutiontoolkit: MIAL super-resolution toolkit v2.0.1,” _Zenodo_ , 2020. [Online]. Available: https://zenodo.org/record/4392788
  * [16] A. Uus et al., “Retrospective motion correction in foetal MRI for clinical applications: existing methods, applications and integration into clinical practice,” _The British Journal of Radiology_ , vol. 95, p. 20220071, 2022\.
  * [17] J. Xu et al., “NeSVoR: Implicit neural representation for slice-to-volume reconstruction in MRI,” _IEEE Transactions on Medical Imaging_ , 2023.
  * [18] M. Ebner et al., “An automated framework for localization, segmentation and super-resolution reconstruction of fetal brain MRI,” _NeuroImage_ , vol. 206, 2020.
  * [19] A. U. Uus et al., “Automated 3d reconstruction of the fetal thorax in the standard atlas space from motion-corrupted mri stacks for 21–36 weeks ga range,” _Medical image analysis_ , vol. 80, p. 102484, 2022.
  * [20] S. Lala et al., “A deep learning approach for image quality assessment of fetal brain mri,” in _ISMRM, Montréal, Québec, Canada_ , 2019, p. 839.
  * [21] J. Xu et al., “Semi-supervised learning for fetal brain mri quality assessment with roi consistency,” in _MICCAI_. Springer, 2020, pp. 386–395.
  * [22] L. Liao et al., “Joint image quality assessment and brain extraction of fetal mri using deep learning,” in _MICCAI_. Springer, 2020, pp. 415–424.
  * [23] B. Gagoski et al., “Automated detection and reacquisition of motion-degraded images in fetal haste imaging at 3 t,” _Magnetic Resonance in Medicine_ , vol. 87, no. 4, pp. 1914–1922, 2022.
  * [24] K. J. Gorgolewski et al., “The brain imaging data structure, a format for organizing and describing outputs of neuroimaging experiments,” _Scientific data_ , vol. 3, no. 1, pp. 1–9, 2016.
  * [25] M. Ranzini et al., “Monaifbs: Monai-based fetal brain mri deep learning segmentation,” _arXiv preprint arXiv:2103.13314_ , 2021.
  * [26] P. de Dumast et al., “Translating fetal brain magnetic resonance image super-resolution into the clinical environment,” _European Congress of Magnetic Resonance in Neuropediatrics_ , 2020.
  * [27] I. Legoretta, S. Samal, and al., “Github repository: Automatic fetal brain MRI quality assessment.” [Online]. Available: https://github.com/FNNDSC/pl-fetal-brain-assessment
  * [28] N. Tustison et al., “N4itk: Improved N3 bias correction,” _IEEE Transactions on Medical Imaging_ , vol. 29, no. 6, 2010.
  * [29] N. Littlestone, “Learning quickly when irrelevant attributes abound: A new linear-threshold algorithm,” _Machine learning_ , vol. 2, pp. 285–318, 1988\.
  * [30] G. Varoquaux et al., “Assessing and tuning brain decoders: cross-validation, caveats, and guidelines,” _NeuroImage_ , vol. 145, pp. 166–179, 2017.
  * [31] F. Pedregosa et al., “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.

## Supplementary material

Table 3: Detailed description of the data used in the. Field refers to the
magnetic field of the scanner, TR is the repetition time and TE is the echo
time. All scanners used a Half-Fourier Acquisition Single-shot Turbo spin Echo
imaging (HASTE) sequence.

SITE1
---
Model | Field [T] | $(n_{\text{subjects}},n_{\text{LR}})$ | TR [ms] | TE [ms] | Resolution [$mm^{3}$]
Siemens Aera | 1.5 | $(34,281)$ | 1200 | 90 | $1.12\times 1.12\times 3.3$
Siemens MAGNETOM Sola | 1.5 | $(17,138)$ | 1200 | 90 | $1.1\times 1.1\times 3.3$
Siemens MAGNETOM Vida | 3 | $(2,14)$ | 1100 | 101 | $0.55\times 0.55\times 3$
Siemens Skyra | 3 | $(8,77)$ | 1100 | 90 | $0.55\times 0.55\times 3$
Siemens Avanto | 1.5 | $(1,5)$ | 1000 | 82 | $1.2\times 1.2\times 4$
SITE2
Model | Field [T] | $(n_{\text{subjects}},n_{\text{LR}})$ | TR [ms] | TE [ms] | Resolution [$mm^{3}$]
Siemens Aera | 1.5 | $(16,158)$ |  |  |
|  | $(6,80)$ | 1500 | 82 | $0.55\times 0.55\times 2.5$
|  | $(4,34)$ | 1000 | 137 | $0.59\times 0.59\times 3.5$
|  | $(4,33)$ | 1000 | 81 | $0.55\times 0.55\times 3.15$
|  | $(2,11)$ | 1200 | 94 | $1.72\times 1.72\times 4.2$
Siemens MAGNETOM Vida | 3 | $(11,56)$ | 1540 | 77 | $1.04\times 1.04\times 3$
Siemens TrioTim | 3 | $(59,322)$ |  |  | // 4 outliers
|  | $(24,97)$ | 1100 | 127 | $0.51\times 0.51\times 3.5$
|  | $(15,108)$ | 990 | 137 | $0.68\times 0.68\times 3.5-6.0$
|  | $(14,71)$ | 2009 | 137 | $0.51\times 0.51\times 3.5$
|  | $(1,14)$ | 3640 | 137 | $0.51\times 0.51\times 3.5$

Table 4: Detailed description of the metrics proposed for FetMRQC.

Intensity-based metrics
---
slice_loss |  Use metrics commonly used for outlier rejection [13, 14, 18] to compute the difference between slices in the volume. We considered (normalized) mean averaged error, (normalized) mutual information, normalized cross correlation, (normalized) root mean squared error, peak signal-to-noise ration, structural similarity and joint entropy.
sstats [8] |  Compute the mean, median, standard deviation, percentiles $5\%$ and $95\%$, coefficient of variation and kurtosis on brain ROI.
entropy [8] |  Measure the overall entropy of the image.
bias |  Level of bias estimated using N4 bias field correction [28]
filter_image |  Estimate the sharpness by using Laplace and Sobel filters (commonly used for edge detection)
Shaped-based metrics
closing_mask |  Morphological closing of the brain mask in the through-plane direction, to detect inter-slice motion. Report the average difference with the original mask.
filter_mask |  Estimate the sharpness of the brain mask using Laplace and Sobel filtering. In an ideal case, the brain mask would be smoothly varying, especially in the through-plane direction.

Figure 8: Correlation matrix between all 75 features, evaluated on the entire dataset. Blue refers to negative correlations, and red to positive ones. Table 5: Parameters automatically optimized by the inner loop of the nested CV. All models were trained with the default parameter from scikit-learn 1.1.3. Model step |  Parameters
---|---
Remove corr- elated features |  $\text{Threshold}\in\\{0.8,0.9\\}$; Disabled
Winnow algorithm | Enabled, Disabled
PCA | Enabled, Disabled
Regression models |  Linear regression, Gradient boosting, Random Forest
Classification models |  Logistic regression, Random Forest, Gradient Boosting, AdaBoost