# Histology to 3D In Vivo MR Registration for Volumetric Evaluation of MRgFUS
Treatment Assessment Biomarkers

Blake E. Zimmerman  Department of Biomedical Engineering, University of Utah,
Salt Lake City, UT  Sara L. Johnson  Department of Biomedical Engineering,
University of Utah, Salt Lake City, UT  Utah Center for Advanced Imaging
Research, University of Utah, Salt Lake City, UT  Henrik A. Odéen  Utah Center
for Advanced Imaging Research, University of Utah, Salt Lake City, UT  Jill E.
Shea  Department of Surgery, University of Utah, Salt Lake City, UT  Rachel E.
Factor  Department of Pathology, University of Utah, Salt Lake City, UT
Sarang C. Joshi  Department of Biomedical Engineering, University of Utah,
Salt Lake City, UT  Scientific Computing and Imaging Institute, University of
Utah, Salt Lake City, UT  Allison H. Payne  Utah Center for Advanced Imaging
Research, University of Utah, Salt Lake City, UT

###### Abstract

Advances in imaging and early cancer detection have increased interest in
magnetic resonance (MR) guided focused ultrasound (MRgFUS) technologies for
cancer treatment. MRgFUS ablation treatments could reduce surgical risks,
preserve organ tissue/function, and improve patient quality of life. However,
surgical resection and histological analysis remain the gold standard to
assess cancer treatment response. For non-invasive ablation therapies such as
MRgFUS, the treatment response must be determined through MR imaging
biomarkers. However, current MR biomarkers are inconclusive and have not been
rigorously evaluated against histology via accurate registration. Existing
registration methods rely on anatomical features to directly register in vivo
MR and histology. For MRgFUS applications in anatomies such as liver, kidney,
or breast, anatomical features independent from treatment features are often
insufficient to perform direct registration. We present a novel MR to
histology registration workflow that utilizes intermediate imaging and does
not rely on these independent features. The presented workflow yields an
overall registration accuracy of 1.00 +/- 0.13 mm. The developed registration
pipeline is used to evaluate a common MRgFUS treatment assessment biomarker
against histology. Evaluating MR biomarkers against histology using this
registration pipeline will facilitate validating novel MRgFUS biomarkers to
improve treatment assessment without surgical intervention.

##  1 Introduction

Improved early cancer detection has driven the development of more
conservative, less invasive cancer treatment alternatives to surgical
intervention. These minimally and non-invasive tumor treatments have the
potential to reduce or eliminate surgical risks, preserve organ tissue and
function, and improve patient quality of life. Magnetic resonance (MR) guided
focused ultrasound (MRgFUS) ablation is a rapidly growing technology for non-
invasive tumor treatment with applications including liver [1], prostate [2],
and breast tumors [3, 4, 5]. In any ablative treatment, gold-standard
determination of treatment response requires tissue resection and histological
analysis [6]. However, one of the key benefits of minimally and non-invasive
ablation treatments is the elimination of traditional surgical intervention.
Therefore, in non-invasive ablation therapies such as MRgFUS, treatment
response must be determined through MR imaging biomarkers alone. Before MR
biomarkers can be used to assess MRgFUS treatment response, they must be
validated against histology to determine accuracy [7]. The most robust way to
validate MR biomarkers is through spatially accurate and precise comparison
between in vivo MR biomarkers and histology.

Directly comparing MR images of in vivo tissue (in vivo MR) to histology is
challenging. Excising treated tissue, slicing into smaller tissue blocks, and
preparing tissue samples for histology all destroy the spatial relationship
between in vivo MR and histology images [8]. Because of this difficulty,
previous studies have avoided restoring this spatial relationship and evaluate
MRgFUS MR biomarkers via qualitative or indirect metrics, such as non-viable
tumor fractions that can be derived relative to the MR and histology spaces
independently [9]. Qualitative analysis and indirect metric evaluation have
been used to evaluate several MR biomarkers, including the commonly used non-
perfused volume (NPV) measured on contrast-enhanced (CE) T1-weighted (T1w) MR
imaging and other biomarkers such as T2w imaging, thermal dose, and diffusion
imaging [9, 10, 11, 12, 13]. Although the current spatially non-specific
qualitative analysis and indirect metrics demonstrate a strong correlation
between some MR biomarkers and histology, the conclusive spatial accuracy
cannot be determined with these methods. For example, conflicting studies show
that the NPV biomarker assessed immediately ($\sim$ 1 hour) after treatment
both over- and under estimates histology lesion size [7, 14]; however, NPV has
increased correlation with histology lesion size several days post treatment
[7, 15]. The spatial relationship between in vivo MR and histology must be
restored to establish the predictive accuracy of in vivo MR biomarkers against
histology.

Several studies have developed MR to histology registration methods to restore
the spatial relationship between in vivo diagnostic or treatment evaluation MR
images and 2D histology imaging [16, 6, 17, 18, 19, 20, 21, 22]. With accurate
registration, in vivo MR biomarkers can be directly compared to histology
using a variety of spatial metrics. Sufficient registration accuracy depends
on the resolution of the biomarker being evaluated. Current in vivo MR
biomarkers are typically acquired with an in-plane resolution of 1-2 mm, for
which the ideal registration accuracy is less than one millimeter for accurate
evaluation [7, 23]. A registration technique with this level of accuracy
allows development of more accurate, novel biomarkers. However, previous
registration methods rely on features that correlate directly between in vivo
MR images and histology that are independent of the diagnostic or treatment
features being evaluated [21]. These independent features can be used to
facilitate registration without influencing the final comparison between MR
treatment features and the corresponding histology response. For example,
registration methods for prostate applications excise the entire organ and use
whole-mount histology to allow a correlation between the organ boundary on
both in vivo MR and histology [16, 17]. Many registration methods additionally
rely on correlating anatomical landmarks to support the registration [6].
However, for MRgFUS ablation applications in liver, breast, or kidney,
features independent from the ablation treatment features are typically not
available for input to registration algorithms. This lack of sufficient
correlating features between in vivo MR and histology independent of the
treatment features prohibits the use of previously developed registration
techniques. Potential treatment dependent features are visible on both in vivo
MR and histology, but the use of these features to drive registration biases
the final result of biomarker evaluation. The need is critical for a new
registration method that uses intermediate steps in lieu of feature
correlation to accurately register in vivo MR and histology to evaluate MRgFUS
treatment biomarkers.

In this paper, we present a novel, multi-step MR to histology registration
workflow that corrects all deformations without assuming any direct feature
correlation between in vivo MR and histology, but rather between three
intermediate stages. We estimate the deformation from every tissue processing
step and compose the deformations to generate a full 3D map between any
histology image and in vivo MR biomarkers. We densely sample the histology
space and map all histology annotations to in vivo MR to form a volumetric
histological annotation of tissue necrosis in the in vivo MR space. This novel
registration workflow facilitates an evaluation of MR biomarkers that was not
possible before, using spatial metrics such as precision, recall, DICE
coefficient, and Hausdorff distance [24]. We demonstrate the capabilities of
this registration pipeline by evaluating the NPV biomarker following MRgFUS
ablation of a VX2 rabbit tumor model at two time points: 1) acute NPV ($\sim$
1 hour after ablation) and 2) post-NPV (3-5 days after treatment). We compare
these spatially quantified results to the findings of qualitative and indirect
findings in the literature. The presented registration method will facilitate
accurately validating a wide range of in vivo MR biomarkers against histology.
Validating in vivo MR biomarkers will increase their clinical viability for
MRgFUS treatment assessment and facilitate minimally invasive treatments
becoming feasible alternatives to surgical intervention for cancer treatments.

Figure 1: Flow chart of the end-to-end workflow required to register in vivo
MR to histology. The steps for destructive histopathology pipeline are
indicated by D1, D2, and D3, and the steps for restoring registration pipeline
are indicated by R1, R2, and R3. Dashed lines indicate information extracted
from digital imaging during the destructive histopathology pipeline for use in
the restoring steps.

##  2 Overview of Approach

(a)(b)(c) Figure 2: Depiction of the block reconstruction process. The
sequential block reconstruction is shown in (a). Two blocks are registered
together using corresponding faces as in (b), where the light blue face is the
target face and does not move, and the light green face is the moving face.
This registration results in a transformation $\varphi$ for the moving block.
The reconstructed block is registered to the red ex vivo surface in (c) to
generate a transformation $\varphi$ between the reconstructed blocks and ex
vivo MR imaging.

The presented registration approach utilizes intermediate imaging to
comprehensively correct for deformations introduced throughout the tissue
processing workflow. Intermediate imaging provides digital representations of
the tissue both before and after specific deformations have been introduced
during processing, eliminating the need for correlating features between in
vivo MR and histology to drive the registration, as has been done in previous
methods. A depiction of the general workflow is shown in Figure 1. The
pipeline developed and described in this work combines previously developed
registration strategies and contributes novel registration methods [6]. The
tissue processing and registration pipelines were evaluated on four subjects
of a large animal tumor model described in section 5.1.

The (D)estructive histopathology pipeline seen in Figure 1 consists of three
main steps that introduce deformation into the tissue: tissue excision (D1),
gross slicing (D2), and microtome sectioning (D3). Surgical excision and
formalin fixation of the treated tissue introduces orientation uncertainty and
non-linear deformation into the tissue sample. The excised tissue is typically
too large to be processed as a whole, so the sample is grossly sliced into
smaller tissue blocks, usually $\sim$3-5mm. The gross slicing step introduces
non-linear deformations to each tissue block and destroys the spatial
relationship between each block and the original volume. In typical histology
processing, after paraffin wax embedding, each block face is trimmed using a
microtome until a full section of tissue is exposed (’facing the block’), and
then sections are collected for histology and the remaining tissue is set
aside [6]. However, the presented method requires microtome sectioning through
the entire tissue block, collecting sections for histology throughout the
block. Microtome sectioning and subsequent staining and mounting cause non-
linear deformations from tissue shearing, tearing, and stretching. The end
product of the destructive histopathology pipeline is a series of stained 2D
microscopic images that must be re-aligned with in vivo MR.

The (R)estoring registration pipeline independently estimates the deformation
between each of three main destructive histopathology steps via blockface
registration (R1), ex vivo registration (R2), and in vivo registration (R3).
Blockface registration (R1) estimates the 2D deformations from microtome
sectioning, staining, and mounting. After blockface registration, each 2D
histology image must be correlated with the 3D tissue volume before gross
slicing. In this step, most other registration methods depend on 2D slice
correlation between MR (in vivo or ex vivo) and histology to restore the
relationship between the original tissue volume and histology [6]. However,
our novel ex vivo registration (R2) step uses 3D models from ex vivo MR and
blockface imaging to embed each histology section into the original 3D tissue
volume without assuming slice-to-slice correspondence or correlation between
MR and histology features. A general overview of this reconstruction step is
shown in Figure 2. We emulate the 2D slice correspondence assumption in our ex
vivo registration step and compare the error of the presented method against
this assumption to show the benefit of the 3D model reconstruction. Finally,
in vivo registration (R3) estimates deformations from tissue excision by
registering 3D models of corresponding features from ex vivo and in vivo MR
images. The deformations between any two steps is represented by a
diffeomorphism, and these deformations can be composed to generate a single
diffeomorphism between multiple stages of imaging. A single diffeomorphism
from histology to in vivo MR inserts each histology section into the in vivo
MR space and enables the generation of fully registered, volumetric
histological necrosis annotations in the in vivo MR space.

The necessary accuracy of registration is dependent on the native resolution
of the evaluation space, which is in vivo MR imaging. Ideally, the
registration accuracy should be at or below the resolution of the acquired MR
biomarkers being evaluated. Common clinical biomarkers are acquired with an
in-plane resolution of 1-2 mm, making 1 mm of total error an ideal constraint
for the total registration error [7, 23]. The accuracy of the presented
registration pipeline was evaluated at every stage using the Euclidean
distance between corresponding landmarks following registration, known as the
target registration error (TRE). These landmarks are only between intermediate
stages, not directly between in vivo MR and histology, and are not used to
drive registration. For example, features to evaluate R3 between in vivo T2w
MR and ex vivo T2w MR can be identified because they are the same imaging
modality. Additionally, landmarks, such as the the edge of the tissue, can be
selected between histology sections and blockface images to evaluate R1.
However, the edge of the tissue is not identifiable on in vivo MR (because the
tissue was trimmed) and the appearance of the T2w MR features in histology may
not be known. Unlike other applications, selecting landmarks that correlate
between in vivo MR and histology is not possible in this application, so we
evaluate the individual stages of registration. The total accuracy of the
pipeline is estimated by summing the errors from each stage to obtain the
cumulative error from registration. To determine the potential impact of the
2D slice projection on the registration accuracy, the landmarks from ex vivo
registration (R2) are deformed using the proposed registration method and the
2D slice assumption method independently and compared. The acute NPV and post-
NPV biomarkers are spatially evaluated against the volumetric histological
annotation of necrosis resulting from registration using spatially specific
metrics, including the precision, recall, DICE coefficient, and Hausdorff
distance.

(a)(c)(b) (d) (e) (f)

Figure 3: MR biomarker to histology registration with accurate spatial
comparisons of the acute NPV (blue), post-NPV (green), and histology necrosis
(pink) labels. (a, d-f) show volumetric overlays of each of the labels in the
in vivo MR space for all four subjects. The contours on the axis walls are
projections for each respective volume from the slice with the largest acute
NPV area in each dimension. The space between each grid line is 5 mm. (b)
Label contours from the subject shown in (a) sampled into the native histology
space. (c) Follow-up contrast-enhanced MR image resampled onto the
corresponding histology slice in (b) with the each of the labels overlaid,
demonstrating that the estimated transformations are invertible.

##  3 Results

The goal of the registration pipeline was to provide a 3D volumetric label of
histological necrosis registered to the in vivo MR imaging space. Figure 3
includes a volumetric comparison between the acute NPV (blue) and post-NPV
(green) MR biomarkers and the histology necrosis (pink) for each of the four
subjects analyzed. The histology volume represents the label of tissue
necrosis against which the acute and post NPV MR biomarkers were independently
compared to determine their accuracy. The spatial metrics calculated from
these volumes are presented in Table 1. The TRE accuracy for each stage and
specific stages under different assumptions is displayed in Figure 4.

Table 1: The total volume and different spatial evaluation metrics for each label, the acute NPV and post-NPV, compared against the histology volume for each subject analyzed. # |  Acute NPV Vol. ($mm^{3}$) |  Post NPV Vol. ($mm^{3}$) |  Histology Vol. ($mm^{3}$) |  Acute NPV Precision |  Post NPV Precision |  Acute NPV Recall |  Post NPV Recall |  Acute NPV DICE |  Post NPV DICE |  Acute NPV Hausdorff ($mm$) |  Post NPV Hausdorff ($mm$)
---|---|---|---|---|---|---|---|---|---|---|---
$1$ | $1{,}317.35$ | $915.74$ | $675.59$ | $0.39$ | $0.62$ | $0.76$ | $0.85$ | $0.52$ | $0.72$ | $7.61$ | $3.17$
$2$ | $2{,}039.13$ | $1{,}688.73$ | $1{,}281.84$ | $0.49$ | $0.69$ | $0.78$ | $0.91$ | $0.60$ | $0.78$ | $6.09$ | $3.31$
$3$ | $3{,}859.13$ | $2{,}627.20$ | $2{,}391.17$ | $0.49$ | $0.86$ | $0.80$ | $0.94$ | $0.61$ | $0.90$ | $7.51$ | $3.05$
$4$ | $5{,}542.83$ | $3{,}969.91$ | $3{,}477.31$ | $0.58$ | $0.81$ | $0.92$ | $0.92$ | $0.71$ | $0.86$ | $12.14$ | $3.07$

###  3.1 Registration Pipeline Accuracy

1234$0$$0.5$$1$$1.5$Subject NumberTRE (mm)Proposed Pipeline Per Stage ErrorR1
Error R2 ErrorR3 ErrorCum. Error

1234$0$$0.5$$1$$1.5$Subject NumberTRE (mm)2D Slice Projection ComparisonR2
Error 2D Slice

Figure 4: Target registration error (TRE) computed on a subject-specific
basis. (left) Subject-specific TREs for the three stages of the registration
pipeline: blockface (R1), ex vivo (R2), and in vivo (R3) registration. The
cumulative error is the summation of the mean of three stages for each
subject. (right) Subject-specific TREs from ex vivo registration (R2) for the
presented method versus the 2D slice correlation method. The dark region of
the 2D slice bar shows the portion of the total error resulting only from the
Z dimension.

The accuracy of each stage of the registration pipeline was independently
evaluated using landmarks; the R1 landmarks were different from the R2
landmarks, which were both different from the R3 landmarks. TREs for the
different stages of the registration pipeline are plotted in Figure 4 for the
four subjects analyzed. The TRE for blockface registration (R1), ex vivo
registration (R2), and in vivo registration (R3) ranged from 0.096 to 0.139 mm
(mean $\pm$ s.d. = 0.117 $\pm$ 0.045 mm), 0.21 to 0.32 mm (0.24 $\pm$ 0.10
mm), and 0.56 to 0.69 mm (0.64 $\pm$ 0.19 mm), respectively. The sum of the
error for the different stages ranged from 0.91 to 1.13 mm (1.00 $\pm$ 0.13
mm) for the four subjects analyzed.

###  3.2 Comparison with 2D Slice Correlation

The ex vivo registration (R2) was run while restricting the Z dimension
(perpendicular to the histology images) motion to only translation, enforcing
the assumption that the blockface and 2D histology slices correlate with an ex
vivo MR slice. The TRE for the presented ex vivo registration method and the
2D slice correlation method can be seen in Figure 4. The presented method
ranged from 0.21 to 0.32 mm (0.24 $\pm$ 0.10 mm) and significantly (p = 0.025,
N=20) outperformed the 2D slice correlation method, which ranged from 0.65 to
1.19 mm (0.84 $\pm$ 0.22 mm). The 3D error from the presented method even
significantly (p=0.015, N=20) outperformed the Z dimension error of the 2D
slice correlation method alone, ranging from 0.44 to 0.54 mm (0.48 $\pm$ 0.05
mm).

###  3.3 NPV Biomarker Evaluation

The total volume of each label and the spatial metrics comparing the acute NPV
and the post NPV MR biomarkers against the histology necrosis label are shown
in Table 1. The post-NPV precision significantly outperformed the acute NPV
precision (0.0079, N=4) when compared to the histology necrosis labels,
ranging from 0.62 to 0.86 (0.74 $\pm$ 0.09) and 0.39 to 0.58 (0.49 $\pm$
0.07), respectively. The post-NPV recall ranged from 0.76 to 0.94 (0.91 $\pm$
0.03) but was not significantly different (p=0.079, N=4) from the acute NPV
that ranged from 0.76 to 0.92 (0.82 $\pm$ 0.06). The DICE similarity
coefficient was significantly better (p=0.011, N=4) between post-NPV and
histology, which ranged from 0.72 to 0.9 (0.82 $\pm$ 0.07), than between the
acute NPV and histology, which ranged from 0.52 to 0.71 (0.61 $\pm$ 0.07).
Finally, the Hausdorff distance for post-NPV ranged from 3.05 to 3.31 mm (3.15
$\pm$ 0.10 mm) and was significantly smaller (p = 0.007, N=4) than the
Hausdorff distance for the acute NPV, ranging from 6.09 to 12.14 mm (8.34
$\pm$ 2.28 mm).

##  4 Discussion

The rigorous spatial evaluation of the acute and post-NPV MR biomarkers
against the gold-standard histology label is only possible because of the
presented MR to histology registration method. This evaluation revealed that
NPV biomarker acquired 3-5 days after treatment was a more accurate predictor
of the MRgFUS induced tissue necrosis than the acute NPV. Although this result
is consistent with prior studies that have used qualitative or indirect
metrics to compare NPV and histology [14, 7], the presented results
quantitatively demonstrate the significant difference between the two metrics.
The acute NPV over predicted the histology label for every subject, as can be
seen in Figure 3, reinforced by the low precision score, high recall score,
and larger volume. Looking at the contour comparisons in Figure 3, the acute
NPVs for subjects 2, 3, and 4 have distinct features that are not present in
the histology label. These large, mislabeled regions are reflected in the
large Hausdorff distances between the acute NPV and histology, with the most
obvious mislabeled structure in subject 4 with a Hausdorff distance of 12.14
mm. These regions and general overestimation of the acute NPV could be a
result of transient effects previously described in the literature, such as
vascular occlusion [7]. These results demonstrate the acute NPV biomarker
measured on CE-T1w MR imaging is not an accurate predictor of histological
tissue viability following MRgFUS ablation for the presented model.

The post-NPV MR biomarker significantly outperformed the acute NPV in every
spatial metric except for recall. This result is expected because the
overestimation yields a high recall score, which explains the lack of a
significant difference between the acute NPV and the post-NPV recall. However,
the mean precision of the post-NPV biomarker was significantly higher,
indicating that the post-NPV biomarker had fewer false positives compared to
the histology label than the acute NPV. The qualitative alignment of the post-
NPV to histology volume in Figure 3 is quantitatively reflected in the
significantly lower and more consistent reported Hausdorff distances relative
to the acute NPV. It is important to note that these spatial metric results
cannot be used to evaluate the accuracy of the registration because the post-
NPV MR biomarker may not perfectly match the histology volume even with
perfect registration. These correlations between in vivo MR and histology are
what we aim to understand and therefore cannot be used to drive or evaluate
the registration without biasing the final result.

The average TRE of 1.00 $\pm$ 0.13 mm resulting from the presented
registration workflow is lower than the TRE reported by state-of-the-art MR to
histology registration methods used in similar applications [22]. The TRE of
the presented registration workflow is on the order or the resolution of
several clinical MR biomarkers, making the workflow suitable for evaluating
these biomarkers [7, 23]. Additionally, the reported registration accuracy is
achieved without relying on features that correlate directly between in vivo
MR images and histology that are independent of the diagnostic or treatment
features being evaluated. R3 registration between in vivo MR and ex vivo MR in
the presented method may include treatment features. Correlating features
between in vivo and ex vivo MR is possible because the image contrast is
consistent, and the deformation is minimal between the two stages. However,
the direct relationship between in vivo and the H&E stained tissue sections is
unknown, so features between in vivo MR and histology may or may not directly
correlate. Therefore, selecting features directly between in vivo MR and
histology for registration would bias the evaluation of MR biomarkers against
the histological response. Additionally, evaluating the registration from
landmarks directly between in vivo MR and histology might not accurately
evaluate the accuracy of registration because the selected features may not
correlate. To address this problem, we use different landmarks for each stage
for registration evaluation as the correlation between two stages is better
defined. Although our method may include treatment features, we minimize the
dependence to one stage and only between images of the same MR modality. Using
intermediate imaging makes the presented registration workflow applicable for
not only improving registration in prior applications, but also validating of
a wider range of applications, including MRgFUS ablations in the liver,
breast, or kidneys.

A direct comparison between registration methods is not possible due to large
differences in the overall workflows and differences in the included
anatomies. Prior registration methods often rely on a 2D correlation between
in vivo MR and histology [16, 6, 17, 18, 19, 20, 21, 22]. However, this study
emulated the 2D slice correspondence assumption in our ex vivo registration
step and compared the error of the presented method against this assumption to
demonstrate the advantage of the 3D block reconstruction method. Prior studies
have shown that this 2D slice correlation assumption limits the accuracy of
the registration to a minimum of $\sim$ 1mm [8]. We find similar limitations
to this assumption in our results. For our pipeline, the 2D slice assumption
introduces significant amounts of error into the registration pipeline, as
shown in Figure 4. Although the X and Y dimension error could be improved, the
amount of unrecoverable error introduced in the Z dimension using the 2D
assumption is still significantly (p = 0.015) more than the total error using
our novel block reconstruction method. The high registration accuracy, novel
block reconstruction method, and independence of correlating MR and histology
features clearly show the advantages of the presented multi-stage registration
workflow.

The presented registration workflow and metric evaluation do not come without
limitations. The time required for tissue processing and data reconstruction
is not trivial. During microtome sectioning (D3), sequentially sectioning
through each tissue block while acquiring digital blockface images every 50
$\mu m$ is currently time consuming. In this work, each tissue block must be
fully sectioned in order to obtain a 3D model for the ex vivo registration
step (R2). However, there are some methods to generate 3D models of the tissue
blocks that are more automated and less time intensive, such as micro-CT.
Micro-CT can be used to extract the 3D surface of each tissue block without
sectioning, leaving the amount of microtome sectioning and blockface imaging
up to the discretion of the user. The postprocessing relies on semi-automatic
segmentation of blockface imaging and histology sections, which is time
consuming at the high resolutions needed to generate high-fidelity models.
However, the neural networks used to provide initial segmentations in this
study will become more accurate with additional samples, leading to more
automatic processing. Although the presented registration workflow requires
more time compared to other workflows, the improved registration accuracy
results in a reduction of the necessary number of samples to effectively
evaluate MR biomarkers [8, 22, 6]. Although the small sample size limits the
statistical significance of the NPV evaluation, this work demonstrates how the
presented registration workflow can be used to rigorously evaluate any kind of
MR biomarkers. Future studies utilizing the presented registration method can
be added to these results to improve the statistical power of the biomarker
evaluation. Additionally, as more data becomes available, current and novel
biomarkers can be prospectively studied to fully understand their accuracy for
predicting histology necrosis.

Despite these limitations, the presented registration workflow facilitated
extensive comparisons between the commonly utilized MRgFUS NPV biomarker and
gold-standard histology for viability. Contextualize the meaning of the
precision and recall scores for MRgFUS ablation treatments of oncological
targets is important. If an MR biomarker has low precision with high recall,
such as the acute NPV biomarker, a large number target pixels may be labeled
as treated, although they are still viable. Relying solely on these MR
biomarkers may result in untreated cancer remaining after MRgFUS ablation.
Although the acute NPV biomarker is not optimal for immediate MRgFUS treatment
assessment, several promising acute MR biomarkers will need to be evaluated
against histology, including T2w imaging, MR temperature imaging, and
diffusion imaging [9, 10, 11, 12, 13]. The presented registration methods can
facilitate evaluating and validating these proposed biomarkers in future
studies.

We presented a novel, multi-step MR to histology registration workflow that
corrects all deformations without assuming any feature correlation between in
vivo MR and histology, but rather between intermediate imaging steps. The
accuracy of the presented registration workflow has a lower error than
previously reported state-of-the-art in vivo MR to histology registration
methods. The novel contribution of our presented workflow is the block
reconstruction step depicted in Figure 2. The presented registration workflow
facilitated novel spatially accurate evaluation of MRgFUS NPV biomarkers
against the gold-standard label of tissue viability and can be used to
evaluate and develop novel, acute in vivo MR metrics that can more accurately
predict the tissue viability measured with histology. Validating in vivo MR
biomarkers will increase their clinical viability for MRgFUS treatment
assessment and facilitate clinical translation of minimally invasive
treatments for oncological applications.

##  5 Methods

The presented MR to histopathology registration pipeline was evaluated in a
large animal tumor model. Supplemental Video 1 provides an animation of the
entire tissue processing and reconstruction workflow. To evaluate both the
acute NPV and the post-NPV as MR imaging biomarkers, the acute NPV was
registered to the post-NPV using methods described in Zimmerman et al. [25]
The tissue from the subject was excised immediately after follow-up imaging
according to the procedures outline in section 5.2-5.4. The histology was
reconstructed and registered to the follow-up in vivo MR images according to
the procedures outlined in section 5.5-5.7.

###  5.1 Subject Model

A VX2 cell suspension (1×106 cells in 50% media/Matrigel) was bilaterally
injected intra-muscularly into the quadriceps of four New Zealand white
rabbits and grown for 1-2 weeks. Anesthesia was induced with a
ketamine/xylazine injection (IM, 25/5 mg/kg), and the animal was then
intubated, allowing anesthesia to be maintained with inhaled isoflurane
(0.5-4.0%) for the duration of the MRgFUS treatment. Animal vitals, including
temperature and respiration, were monitored throughout treatment. Hair on the
treated quadriceps was removed via clippers and a depilatory cream (Nair™) to
enable acoustic coupling. Using a pre-clinical MRgFUS system (Image Guided
Therapy, Inc.), ablation was performed on one tumor and the surrounding
quadriceps muscle tissue with a 256-element phased-array transducer (Imasonic,
Voray-sur-l’Ognon, France; 10-cm focal length, 14.4 x 9.8 cm aperture, f=940
kHz) inside a 3T MRI scanner (PrismaFIT Siemens, Erlangen, Germany). Ablation
details for the four subjects are outlined in Table 2. A single loop MR
receiver coil was incorporated into the MRgFUS system table to improve the
image signal to noise ratio (SNR) around the targeted quadriceps. The ablation
procedure was monitored in real time with MR thermometry imaging (MRTI) using
3D-segmented echo planar imaging sequences. Following ablation therapy, the
animal was recovered and monitored for 3 days. After 3 days, the animal was
re-anesthetized and follow-up (post) imaging was performed. For more details
on the in vivo MR imaging, we refer the reader to Zimmerman et al.[25] The
animal was euthanized immediately following imaging. The University of Utah
Institutional Animal Care and Use Committee (IACUC) approved all procedures
(#17-08012, September 7, 2017). All methods and procedures were performed in
accordance with the IACUC guidelines and regulations.

1.1  Follow-up Number of Acoustic Power Total #. Duration Sonications Mean
$\pm$ 1 Std. (W) Energy (kJ) 1 5 days 11 57 $\pm$ 17 23.14 2 3 days 12 69
$\pm$ 18 26.00 3 5 days 14 44 $\pm$ 9 18.59 4 5 days 10 56 $\pm$ 9 18.55

Table 2: Ablation details for each subject including time between treatment MR
(where the acute NPV is measured) and follow-up imaging (where post-NPV is
measured), ablation points, acoustic power output, and total energy achieved
for each subject.

###  5.2 D1: Tissue Excision

The treated quadriceps was surgically excised immediately following
euthanization. Pathology inks were applied to the tissue during excision to
maintain the MR orientation, and the excised tissue was submerged in 10%
formalin solution for 14 days. The fixed ex vivo tissue was mounted in a
custom tissue-processing box and encapsulated in 3.5% agar solution for ex
vivo imaging and subsequent gross slicing. The custom box and pathology inks
facilitated orienting the tissue as close to the in vivo MR orientation as
possible. The tissue-processing box incorporated a single-loop MR coil around
the box to improve ex vivo MR image SNR. T1w and T2w MR images of the agar-
embedded tissue were acquired using a 3T MRI scanner (PrismaFIT Siemens,
Erlangen, Germany). The field of view and voxel size for in vivo and ex vivo
T1w MR images used in the registration pipeline was 256 x 56 x 192 mm with 0.5
x 1.0 x 0.5 mm spacing, and for T2w MR the field of view and spacing was 256 x
52 x 192 with $1.0$ mm isotropic spacing.

###  5.3 D2: Gross Slicing

Following ex vivo imaging, the agar embedded tissue block was grossly sliced
along the head-foot axis of the ex vivo MR imaging in $\sim$3 mm increments
using a deli slicer (Backyard Pro SL110E). The surrounding agar was removed
from each sliced tissue block without disturbing the tissue blocks. The
exposed tissue faces from gross slicing were re-inked with pathology inks to
indicate the head and foot surfaces and maintain the orientation of each block
relative to the original ex vivo sample. Each tissue block was placed into an
individual whole-mount tissue cassette for further formalin fixation (5 days).
After the additional fixation, each block was embedded in paraffin wax with
consistent orientation.

###  5.4 D3: Microtome Sectioning

Each paraffin wax tissue block was sectioned at 10 $\mu$m increments with a
microtome (Leica RM2255, Leica Microsystems, Wetzlar, Germany). Digital images
of the blockface were acquired every 50 $\mu$m starting from the very
beginning of the block using a digital single-lens reflex (DSLR) camera (Nikon
D7100; Macro 1:1 105mm Lens; 2.0x teleconverter). The image size and
approximate resolution for a blockface image are $6000x4000$ with $\sim 0.018$
mm isotropic spacing. Paraffin wax is slightly transparent, so tissue from
behind the exposed face would show through on the images. To address this, two
images were acquired at every 50 $\mu$m with two different lighting
conditions. The first image had the light approximately aligned with the
camera whereas the second image had the light perpendicular to the camera.
Taking the difference between two lighting conditions shows only the tissue
that is exposed on the surface. These blockface images were acquired
automatically with an Arduino board controlled with a Python GUI. The Python
GUI automatically triggered the camera and lights, transferred files to the
computer, recorded image information (name, section depth, camera settings,
etc), and backed up each image to cloud storage to ensure data retention.
Three 5 $\mu$m thick tissue sections were retained on glass slides for future
histopathology analysis every 250 $\mu$m. Sequential sampling was repeated
until there was no tissue remaining in the paraffin wax. After sectioning, the
first section of each group of three was stained with hematoxylin and eosin
(H&E). The remaining sections were reserved for additional future stains. All
H&E stained sections were imaged with a brightfield microscope (Axio Scan,
Zeiss, Oberkochen, Germany) at 2.5 magnification. The native resolution of the
microscopic images is $\sim 0.0076mm$, but the images were down-sampled to the
resolution of the blockface images for use in registration. Labels of necrotic
tissue on the H&E stained sections were semi-automatically generated using
Gaussian mixture modeling and expert manual segmentation.

###  5.5 R1: Histology to Block Registration

Each histopathology section was registered to the nearest incremental 50
$\mu$m blockface image. The digital histology images were down-sampled to the
resolution of the blockface image ($\sim$ 0.02$\times$0.02$\times$0.05 mm). An
affine transform between the histology and blockface was solved via automatic
intensity-based affine registration of the blockface segmentation and the
histology segmentation. Following affine registration, a multi-scale,
intensity-based registration was used to deformably register the
histopathology segmentation images to the corresponding blockface image
segmentations. Each image was registered by minimizing the following energy
$E$:

$E=\int_{\Omega}\left\lVert
I_{1}(\varphi^{-1}(\vec{x},t))-I_{0}(\vec{x})\right\rVert^{2},$ (1)

where $I_{1}$ is the histopathology segmentation, $I_{0}$ is the blockface
image segmentation, $\varphi^{-1}(\cdot)$ is a diffeomorphism, and $\Omega$ is
the image domain. Equation 1 was optimized using a gradient flow algorithm
with a Cauchy-Navier operator [26]. This registration step provided a
diffeomorphism between a histology section and its corresponding blockface
image that corrected deformations introduced during microtome sectioning.

###  5.6 R2: Tissue Blocks to Ex Vivo Registration

Surface registration techniques were used to drive the restoration of the
blocks to their original morphology and account for deformations from gross
slicing. A 3D surface, represented by a triangular mesh object, was
constructed for each tissue block and the ex vivo tissue from their
segmentations using a marching cubes algorithm. The exterior tissue surface
was segmented from ex vivo MR to generate a 3D surface of the ex vivo tissue
that was used as the target for block reconstruction. For each tissue block,
the 2D blockface images were stacked to create a 3D blockface volume.
Automatic intensity-based 2D affine registration was used to register
sequential blockface images to account for small shifts in camera position
during imaging. Finally, blockface volumes were semi-automatically segmented
using a custom-trained 3D V-Net neural network and manual segmentation to
generate the 3D surfaces of each tissue block [27].

Each tissue block was represented by a surface, $S_{b}$ where $b$ indicates
the block number ranging from one to the number of total blocks from the gross
slicing, $nb$. For a single tissue block, the mesh was semi-automatically
separated using surface normals and manual segmentation into a head surface
$S_{b}^{h}$, foot surface $S_{b}^{f}$, and exterior surface $S_{b}^{e}$, such
that $S_{b}=S_{b}^{h}\cup S_{b}^{f}\cup S_{b}^{e}$. Examples of the
corresponding faces used for registration are shown in Figure 2 (b). The
corresponding surfaces were registered together using surface-based
registration, which is outlined in section 5.8. The algorithm in Figure 5
provides an overview of the block reconstruction to ex vivo process.
$d(\cdot,\cdot)$ refers to the surface distance metric defined in section 5.8.
In general, a center block was chosen as the starting point for the
reconstruction, and deformations were solved for block by block, propagating
outwards from the center block first in the head direction, then in the foot
direction. Once the deformation was determined for one block, the deformed
block then became the target of registration for the next sequential block.
Each block had an associated affine transformation and diffeomorphism, $A_{b}$
and $\varphi^{-1}_{b}$, respectively. After all blocks were registered and
stacked back together, the deformed exterior surfaces from each block were
joined together and registered with the ex vivo tissue surface.

1:procedure Rebuild Blocks($\\{S_{1}\dots S_{nb}\\}$)

2: Center Block for Target $S_{t}$where $t=nb\;//\;2$

3: for i = 1; i $\leq$ nb // 2; i++ do

4: b $=$ (nb // 2) - i

5: $\hat{S}_{b}^{f}=\min_{A_{b}}d(\>S_{t}^{h},\>A_{b}\circ S_{b}^{f}\>)$

6:
$\bar{S}_{b}^{f}=\min_{\varphi^{-1}_{b}}d(\>S_{b}^{h},\>\varphi^{-1}_{b}\circ\hat{S}_{b}^{f}\>)$

7: $S_{t}^{h}=\varphi^{-1}_{b}\circ A_{b}\circ S_{b}^{h}$

8: end for

9: Reset $S_{t}$ to center block $t=nb\;//\;2$

10: for i = 1; i $\leq$ nb // 2; i++ do

11: b $=$ (nb // 2) + i

12: $\hat{S}_{b}^{h}=\min_{A_{b}}d(\>S_{t}^{f},\>A_{b}\circ S_{b}^{h}\>)$

13:
$\bar{S}_{b}^{h}=\min_{\varphi^{-1}_{b}}d(\>S_{b}^{f},\>\varphi^{-1}_{b}\circ\bar{S}_{b}^{h}\>)$

14: $S_{t}^{f}=\varphi^{-1}_{b}\circ A_{b}\circ S_{b}^{f}$

15: end for

16: Union deformed exterior surfaces $S^{e}=\bigcup\\{\bar{S}_{b}^{e};b\in
np\\}$

17: $\hat{S}^{e}=\min_{A}d(\>S_{\textit{ex vivo}},\>A\circ S^{e}\>)$

18: $\bar{S}^{e}=\min_{\varphi^{-1}}d(\>S_{\textit{ex
vivo}},\>\varphi^{-1}\circ\hat{S}^{e}\>)$

19:end procedure

Figure 5: Algorithm overview for reconstruction of the tissue blocks to the ex
vivo tissue surface.

###  5.7 R3: Ex Vivo to In Vivo Registration

Expert segmentations of corresponding features, such as blood vessels, tumor,
or treatment features, identified from the T2w images from the ex vivo and in
vivo MR data were used to generate 3D surfaces in both spaces. An affine
transform and diffeomorphism between ex vivo and in vivo MR was solved using
3D surface-based registration defined in section 5.8. The full transformation
from the in vivo to ex vivo environment was comprised of an affine matrix and
a diffeomorphism. All deformations from the three restoring registration
pipeline steps were composed to yield a final diffeomorphism that was a
spatial mapping between in vivo MR and any histopathology section.

###  5.8 Surface-Based Registration

Surface-based registration was used throughout the reconstruction pipeline to
solve for affine transforms and diffeomorphisms between two 3D surfaces. The
surfaces did not have corresponding landmarks that could be identified from
the images or surfaces. We implemented a previously developed unlabeled point-
set matching algorithm to solve for a diffeomorphism between the two surfaces
[28]. All surface objects were defined as triangular mesh objects. For a
single triangular mesh object $S$, the kernel norm (K-norm) is defined as

$\left\lVert\sum_{p=1}^{n\in
S}\eta(p)\delta_{c(p)}\right\rVert_{K}^{2}=\sum_{p=1}^{n}\sum_{p^{\prime}=1}^{n}\;\langle\eta(p),\eta(p^{\prime})\rangle\;K(c(p),c(p^{\prime}))$
(2)

where $p$ is a triangular element of $S$, $n\in S$ is the number of faces in
surface $S$, $\eta(p)$ the normal vector to $p$, $c(p)$ is the center point of
triangle $p$, $K(\cdot,\cdot)=k(\cdot,\cdot)I$ where $I$ is a $3\times 3$
identity matrix and $k(\cdot,\cdot)$ is a scalar valued Cauchy kernel. The
dissimilarity between two surfaces $S_{1}$, $S_{2}$ is driven by the
difference between the sum of the vector valued Dirac masses centered at the
triangle center with

$\begin{split}\left\lVert\sum_{p=1}^{n\in
S_{1}}\eta(p)\delta_{c(p)}-\sum_{q=1}^{m\in
S_{2}}\eta(q)\delta_{c(q)}\right\rVert_{K}^{2}=&\sum_{p=1}^{n}\sum_{p^{\prime}=1}^{n}\;\langle\eta(p),\eta(p^{\prime})\rangle\;K(c(p),c(p^{\prime}))-2&\sum_{p=1}^{n}\sum_{q=1}^{m}\;\langle\eta(p),\eta(q)\rangle\;K(c(p),c(q))\\\
+&\sum_{q=1}^{m}\sum_{q^{\prime}=1}^{m}\;\langle\eta(q),\eta(q^{\prime})\rangle\;K(c(q),c(q^{\prime})).\end{split}$
(3)

For a more detailed explanation on how affine transformations and
diffeomorphisms act on these surfaces objects, we refer the reader to Glaunes
et al [28].

###  5.9 NPV Biomarker Analysis

After performing registration, the acute NPV biomarker, post-NPV biomarker,
and the histology necrosis label were co-registered in the follow-up MR
imaging space (where post-NPV was measured). The spatial accuracy of the acute
NPV and the post-NPV biomarkers were evaluated against the histology necrosis
label using precision, recall, DICE coefficient, and Hausdorff distance. The
acute NPV vs histology and post-NPV vs histology spatial metrics were compared
using two-sample-independent t-tests.

###  5.10 Registration Accuracy

Landmarks were used to validate the accuracy of the registration via the
target registration error (TRE). Selecting landmarks directly between MR and
histology is difficult for two reasons: first, selecting landmarks between 2D
and 3D spaces is challenging, and the relationship between MR treatment
features and histology treatment features is being evaluated and cannot be
used to evaluate the registration. Consequently, we selected three sets of
N=20 landmarks (five per subject) for each of the four subjects to evaluate
each step of the registration pipeline individually. The TRE was calculated
via the Euclidean distance between the deformed landmarks and the target
landmarks for each stage of registration. The upper bound of the total
registration pipeline was estimated by summing the mean error from each
individual stage.

The novel contribution of our registration pipeline is the 3D block
reconstruction to insert each histology slice into the 3D ex vivo MR space.
Prior registration methods assumed a 2D relationship between 2D histology and
a 2D MR slice. We emulate this assumption during our 3D block reconstruction
by restricting the Z dimension (perpendicular to the histology sections) block
motion to only translation, effectively assuming that the 2D histology slices
align with a 2D ex vivo MR slice. Using the landmarks, we calculate the TRE
under this 2D projection assumption to demonstrate the improved accuracy of
the presented 3D block reconstruction. The TREs of the two different methods
were compared using two-sample-related t-tests on the Euclidean error between
the deformed landmarks and target landmarks.

## References

  * [1] Kennedy, J. _et al._ High-intensity focused ultrasound for the treatment of liver tumours.  _Ultrasonics_ 42, 931–935 (2004).
  * [2] Colombel, M. & Gelet, A.  Principles and results of high-intensity focused ultrasound for localized prostate cancer.  _Prostate cancer and prostatic diseases_ 7, 289–294 (2004).
  * [3] Merckel, L. G. _et al._ Eligibility of patients for minimally invasive breast cancer therapy based on mri analysis of tumor proximity to skin and pectoral muscle.  _The breast journal_ 24, 501–508 (2018).
  * [4] Sun, F. _et al._ _Therapies for Clinically Localized Prostate Cancer_ (Agency for Healthcare Research and Quality (US), 2014).
  * [5] Kaiser, W. A., Pfleiderer, S. O. & Baltzer, P. A.  Mri-guided interventions of the breast.  _Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine_ 27, 347–355 (2008).
  * [6] Pichat, J., Iglesias, J. E., Yousry, T., Ourselin, S. & Modat, M.  A survey of methods for 3d histology reconstruction.  _Medical image analysis_ 46, 73–105 (2018).
  * [7] Hectors, S. J., Jacobs, I., Moonen, C. T., Strijkers, G. J. & Nicolay, K.  Mri methods for the evaluation of high intensity focused ultrasound tumor treatment: Current status and future needs.  _Magnetic resonance in medicine_ 75, 302–317 (2016).
  * [8] Gibson, E. _et al._ 3d prostate histology image reconstruction: quantifying the impact of tissue deformation and histology section location.  _Journal of pathology informatics_ 4 (2013).
  * [9] Hectors, S. J., Jacobs, I., Strijkers, G. J. & Nicolay, K.  Multiparametric mri analysis for the identification of high intensity focused ultrasound-treated tumor tissue.  _PloS one_ 9, e99936 (2014).
  * [10] Plata, J. C. _et al._ A feasibility study on monitoring the evolution of apparent diffusion coefficient decrease during thermal ablation.  _Medical physics_ 42, 5130–5137 (2015).
  * [11] Mannelli, L. _et al._ Assessment of tumor necrosis of hepatocellular carcinoma after chemoembolization: diffusion-weighted and contrast-enhanced mri with histopathologic correlation of the explanted liver.  _American Journal of Roentgenology_ 193, 1044–1052 (2009).
  * [12] Haider, M. A. _et al._ Dynamic contrast-enhanced magnetic resonance imaging for localization of recurrent prostate cancer after external beam radiotherapy.  _International Journal of Radiation Oncology* Biology* Physics_ 70, 425–430 (2008).
  * [13] Wu, Q., Whitman, G. J., Fussell, D. S. & Markey, M. K.  Registration of dce mr images for computer-aided diagnosis of breast cancer.  In _2006 Fortieth Asilomar Conference on Signals, Systems and Computers_ , 826–830 (IEEE, 2006).
  * [14] Payne, A. _et al._ In vivo evaluation of a breast-specific magnetic resonance guided focused ultrasound system in a goat udder model.  _Medical physics_ 40 (2013).
  * [15] Wijlemans, J. W. _et al._ Evolution of the ablation region after magnetic resonance–guided high-intensity focused ultrasound ablation in a vx2 tumor model.  _Investigative radiology_ 48, 381–386 (2013).
  * [16] Li, L. _et al._ Co-registration of ex vivo surgical histopathology and in vivo t2 weighted mri of the prostate via multi-scale spectral embedding representation.  _Scientific reports_ 7, 1–12 (2017).
  * [17] Losnegrd, A. _et al._ Intensity-based volumetric registration of magnetic resonance images and whole-mount sections of the prostate.  _Computerized Medical Imaging and Graphics_ 63, 24–30 (2018).
  * [18] Iglesias, J. E. _et al._ Joint registration and synthesis using a probabilistic model for alignment of mri and histological sections.  _Medical image analysis_ 50, 127–144 (2018).
  * [19] Orczyk, C. _et al._ Preliminary experience with a novel method of three-dimensional co-registration of prostate cancer digital histology and in vivo multiparametric mri.  _Clinical radiology_ 68, e652–e658 (2013).
  * [20] Dickinson, L. _et al._ Image-directed, tissue-preserving focal therapy of prostate cancer: a feasibility study of a novel deformable magnetic resonance-ultrasound (mr-us) registration system.  _BJU international_ 112, 594–601 (2013).
  * [21] Goubran, M. _et al._ Registration of in-vivo to ex-vivo mri of surgically resected specimens: a pipeline for histology to in-vivo registration.  _Journal of neuroscience methods_ 241, 53–65 (2015).
  * [22] Rusu, M. _et al._ Registration of presurgical mri and histopathology images from radical prostatectomy via rapsodi.  _Medical Physics_ (2020).
  * [23] Schmitz, A. C. _et al._ Precise correlation between mri and histopathology–exploring treatment margins for mri-guided localized breast cancer therapy.  _Radiotherapy and Oncology_ 97, 225–232 (2010).
  * [24] Taha, A. A. & Hanbury, A.  Metrics for evaluating 3d medical image segmentation: analysis, selection, and tool.  _BMC medical imaging_ 15, 29 (2015).
  * [25] Zimmerman, B. _et al._ Learning multiparametric biomarkers for assessing mr-guided focused ultrasound treatments.  _IEEE Transactions on Biomedical Engineering_ (2020).
  * [26] Christensen, G. E., Rabbitt, R. D. & Miller, M. I.  Deformable templates using large deformation kinematics.  _IEEE transactions on image processing_ 5, 1435–1447 (1996).
  * [27] Milletari, F., Navab, N. & Ahmadi, S.-A.  V-net: Fully convolutional neural networks for volumetric medical image segmentation.  In _2016 Fourth International Conference on 3D Vision (3DV)_ , 565–571 (IEEE, 2016).
  * [28] Glaunes, J., Trouvé, A. & Younes, L.  Diffeomorphic matching of distributions: A new approach for unlabelled point-sets and sub-manifolds matching.  In _Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004._ , vol. 2, II–II (IEEE, 2004).


# Recent Advances in Software Effort Estimation using Machine Learning

Víctor Uc-Cetina
<EMAIL_ADDRESS>
Facultad de Matemáticas
Universidad Autónoma de Yucatán

###### Abstract

An increasing number of software companies have already realized the
importance of storing project-related data as valuable sources of information
for training prediction models. Such kind of modeling opens the door for the
implementation of tailored strategies to increase the accuracy in effort
estimation of whole teams of engineers. In this article we review the most
recent machine learning approaches used to estimate software development
efforts for both, non-agile and agile methodologies. We analyze the benefits
of adopting an agile methodology in terms of effort estimation possibilities,
such as the modeling of programming patterns and misestimation patterns by
individual engineers. We conclude with an analysis of current and future
trends, regarding software effort estimation through data-driven predictive
models.

Keywords— effort estimation, agile software, stories estimation

A major goal of software project managers and software developers is to
calculate accurate estimates of the effort required to complete a software
development project. Effort prediction in software engineering projects is
considered by some software development managers as a mission impossible kind
of task. Some other managers think it is certainly possible to accomplish it,
at least for specific categories of software applications, if enough
historical information is available and the right methodology is chosen [7].
The common scenario is that each software company ends up developing its own
strategy over time, depending on the kind of applications they develop. In
this article we review the most recently proposed machine learning methods for
software effort prediction, and we categorize them in non-agile and agile
methods.

It is well known that machine learning (ML) has become a very robust data-
driven tool for solving a variety of problems, ranging from face recognition
to chatbots development, from product recommendation to protein structure
prediction. We have been noticing significant impact in different fields such
as computer vision, natural language processing, e-commerce, bioinformatics,
robotics, audio processing, etc. In general applications and fields benefiting
from ML algorithms have been multiplying in the last decade. The same can be
said about ML applications for software engineering problems. In a recent
survey [20], five software engineering activities were clearly identified as
the most investigated using deep ML, namely,

  * •

Software design,

  * •

Software implementation,

  * •

Software testing and debugging,

  * •

Software maintenance, and

  * •

Software management.

Each of these activities include a different number of more specific research
topics. For instance, in software management, two are the most common topics:

  * •

Software repository mining, and

  * •

Effort cost prediction.

Predicting projects efforts is by itself a difficult task. Estimation or
prediction of software development effort is not an exception. In order to
estimate with the highest of prediction ability we need to exploit past
information, we need to develop data-driven strategies. As Jørgensen argues
[6], historical data relevant to the target project improves estimation
accuracy. Moreover, based on recently published well-supported research
results [5], he says that, what we currently know about software effort and
cost estimation is not enough to solve the estimation challenges in the
software industry. Nevertheless, it helps us to identify several actions that
might improve estimation accuracy.

In this article we are interested in the effort prediction problem for whole
projects or large modules of software (non-agile prediction), and also in
agile software effort prediction.

##  1 NON-AGILE EFFORT PREDICTION

In one of the main works on effort prediction carried out two decades ago
Kitchenham et al. [7] studied the effort and duration of 147 projects,
including maintenance and development ones. This data set was obtained from
one software company of outsourcing. They discovered that 63% of the estimates
were calculated within 25% of the actual value, with an average absolute error
of 0.26. In other words, this study supports the idea that estimation of
software development efforts is achievable with a reasonable margin of error,
small enough to consider it as a useful practice in many cases.

One of the earliest studies that used machine learning for building estimators
of software development effort from historical data was presented by
Srinivasan et al. [14]. In this work, regression trees and multilayer
perceptrons algorithms were compared against traditional parametric models for
effort prediction. In both cases the results show that machine learning
algorithms accomplished competitive performances. One of the main advantages
of the methods studied is their adaptability to new data sets.

Jørgensen et al. [6] systematically reviewed 304 software cost estimation
articles, classifying them according to research topic, estimation approach,
research approach, study context, and use of data sets. Based on this
systematic categorization of papers, the authors make some recommendations
over changes they identified as being useful in estimation research.

There are some data-driven models typically studied to solve the effort
prediction software. Two categories that have been compared are linear
regression models and neural networks, being the second type more complex
models with many more parameters to be estimated. Lopez-Martin et al. [8]
performed a comparison of these two categories of algorithms. Specifically
they compared a multiple linear regression methods with multilayer perceptrons
and radial basis function networks. Their conclusion is somewhat expected:
neural networks showed higher accuracy than multiple linear regression.

In a similar approach Bisi et al. [4] tried a methodology based on
evolutionary computation combining swarm optimization for training, principal
component analysis to reduce the dimension of the input data, and finally
genetic algorithms to optimize the weights of a neural network, also reporting
competitive performances.

More recently, BaniMustafa et al. [1] evaluated 3 machine learning models
using the COCOMO NASA data set which includes 93 software projects. The
algorithms used in this study are naïve Bayes, logistic regression, and random
forests. Although these are not by any means state-of-the-art algorithms, the
approach is of general applicability and provide a baseline of methods to
start with and further increase the robustness of our methodology for effort
prediction. The data set includes 24 attributes such as the category of the
application (with 14 possibilities, including avionics, communications,
science, and simulation), development mode (indicating one of the following
options: embedded, organic, semidetached), lines of code measure and the
actual effort in months. The experimental results show a competitive
performance of all 3 methods over the COCOMO model, being the random forest
model the best one. This is not surprising considering that random forest is a
much more robust method compared to the two other tested algorithms.

Effort predictions based on team size has been also approached by Rai et al.
[11]. Recently, they proposed the use of the constructive cost model known as
COCOMO [2] together with support vector regression, to predict effort
estimation.

##  2 AGILE EFFORT PREDICTION

Effort estimation in agile projects focuses on estimating the effort required
for completing user stories. In a recent article, Choetkiertikul et al. [3]
proposed a prediction model for estimating story points based on the
combination of long short-term memory and recurrent highway networks. Their
empirical evaluation with 16 open source projects demonstrates that this
approach consistently outperforms baselines such as random guessing, mean, and
median methods and six other more complex methods such as Doc2Vec and random
forests. Also, Ochodek et al. [9] investigated several neural networks testing
them with 437 use cases from 27 software development projects. The main goal
of this study was to develop an easy-to-train method, which is comparable in
performance to the existing methodologies. The experimental results shows as
the best performer a convolutional neural network used together with a word-
embedding model.

Phannachitta et al. [10] carried out a systematic comparison of software
effort estimators on 13 standard benchmark datasets. Performance metrics and
statistical test methods were used together with different machine learning
algorithms. The experimental results show that an analogy-based model that is
adapted through a gradient boosting machine algorithm and a traditional
adaptation method based on productivity adjustment is the best performer.

Sarro et al. [12] introduced the approach called learning from mistakes. The
main idea is focused on automatically learning from past estimation errors
made by human experts. These errors are used to improve the accuracy of the
predictions by calculating possible future errors in estimates. This empirical
investigation included 402 maintenance and development projects. The main
finding of this study is that the type, severity and magnitude of errors are
all predictable. The algorithms used in this work to predict errors are
classification and regression trees, k-nearest neighbors, naïve Bayes, and
linear programming.

Tawosi et al. [18] replicated and extended a study [13] on multi-objective
software effort estimation in order to increase the confidence in those
previous results. Later on, Tawosi et al. [17] carried out a replication study
of the analysis previously presented by Ochodek et al. [9]. This time the
authors used a larger dataset. However, the results they obtained are not
supporting the hypothesis stating that deep neural networks can significantly
outperform less sophisticated methods for predicting actual efforts in
software development. Moreover, they concluded that semantic analysis of the
texts describing the stories is not enough and it is prone to introduce errors
in the estimations.

There are a couple of datasets usually employed for researching effort
prediction in agile software development projects: the Choetkiertikul dataset
[3] and the TAWOS dataset [16].

The Choetkiertikul dataset was created in 2019 as part of a study that
investigated the use of a deep learning model for estimating story points.
This dataset is stored in a csv format file and contains 23,313 issues from 16
open-source projects. The data can be mainly used for story points-based
estimation.

On the other hand, the TAWOS dataset is a dataset created in 2022 from Jira
repositories and it includes 44 agile open-source software projects with a
total of more than 500,000 issues. It is an easy-to-use dataset stored as a
relational database that can be used with MySql workbench for investigating
different software engineering problems, i. e. agile software effort
estimation. For example, Tawosi et al. [15] used that dataset to investigate
how clustering could be used to estimate story points in agile software
development projects.

##  3 NON-AGILE VS. AGILE

Effort estimation in non-agile project is usually carried out at the level of
the whole project or whole modules, which increases the prediction error,
simply due to the amount of lines of codes that need to be considered. When
the project is developed using an agile methodology, the prediction can be
performed at the level of individual engineers, passing through teams-level,
up to module and project-level, as summarized in Figure 1.

Figure 1: Benefits of performing effort prediction under an agile software
development methodology.

One way of achieving more precise effort estimation in agile development
projects is by keeping track of stories effort estimations and stories actual
efforts needed by engineers of software teams. We argue that when this
information is stored, the following benefits are available:

  1. 1.

Calculation of the current misestimation rate of each engineer.

  2. 2.

Possibility of giving data-driven feedback when new estimations are needed.

  3. 3.

Identification of those engineers with above average skills in you team or
company, for specific tasks or technologies.

Moreover, we believe that:

  1. 1.

As software engineers become more experienced in designing stories, their
stories gradually become more compact and are more easily attainable. We
consider this as a highly desirable design skill.

  2. 2.

As a result of being able to design more compact stories, stories become more
predictable in terms of needed effort and therefore software engineers
gradually reduce their effort misestimation.

Consequently, it is reasonable to expect that over the years, the percentage
of misestimated stories by an engineer will decrease. In other words, there
are patterns inherent to the process that an engineer follows to design
stories and estimate their complexity and those patterns can be learned using
machine learning models.

##  4 CURRENT AND FUTURE TRENDS

We can identify four main families of methods used in recent approaches for
competitively predicting software development effort, all of them with
applications in non-agile and agile software effort prediction (See Figure 2):

Figure 2: Four main families of methods recently used for software effort
prediction in non-agile and agile methodologies.

  * •

Evolutionary algorithms,

  * •

Decision trees algorithms,

  * •

Shallow and deep neural networks, and

  * •

Others.

On the other hand the predictions in all these approaches are calculated
differently, depending on the historical information available. Some of the
features used as inputs for the predictive models are:

  * •

Attributes describing the complexity (e.g. lines of codes),

  * •

Attributes describing the type of software application (e.g. mobile app),

  * •

Names of use cases,

  * •

Features extracted using analogy-based methods,

  * •

Time required to accomplish stories, and

  * •

Complexity points of user stories.

In spite of the fact that deep neural networks are currently a dominating
force in practically all domains where data-driven predictions are needed, we
argue that in the particular case of effort prediction of agile software
stories, even simple methods such as linear regression models can be very
effective and computationally more suitable to be implemented and maintained
in software development management systems.

However, we expect to see more neural network models being applied to the
effort prediction problem. For instance, transformer architectures [19], based
on attention layers, have revolutionized the language processing domain. They
are currently being investigated for different tasks in the computer vision
domain, and we will not have to wait much before transformers start to be
applied as predictive models in software engineering tasks, such as effort
prediction.

In recent years, more and more companies have started using various types of
project management software to keep track of their daily development tasks and
as a result of this practice, larger data sets are being collected. Those data
sets become over time an important asset, allowing companies to generate more
accurate prediction models in the future.

From our point of view, in the coming years we will witness an increment in
tracking and management software systems with a diversity of data-driven
prediction features. The prediction of efforts at different levels of
complexity, ranging from individual software stories up to whole modules or
projects will definitely be one of these features.

##  5 CONCLUSION

Effort prediction of software development is attainable with some degree of
accuracy, making it a viable practice in industry. Modern neural network
models are capable to learn from historical data of software development
projects and extract critical patterns that allow them to predict efforts in
similar projects.

However, when the project to be analyzed for effort prediction is particularly
different to the projects in our available historical data, then the accuracy
in prediction decreases significantly. Therefore more research is still needed
to solve this kind of extrapolation in the estimation of project efforts.

Finally, since an increasing number of software companies are adopting an
agile software development methodology, effort estimation methods focusing on
tracking individual engineers development patterns will be developed.

## References

  * [1] A. BaniMustafa.  Predicting software effort estimation using machine learning techniques.  In 2018 8th International Conference on Computer Science and Information Technology (CSIT), pages 249–256, 2018.
  * [2] B. W. Boehm.  Software Engineering Economics.  Prentice Hall, 1981.
  * [3] M. Choetkiertikul, H. K. Dam, T. Tran, T. Pham, A. Ghose, and T. Menzies.  A deep learning model for estimating story points.  IEEE Transactions on Software Engineering, 45(7):637–656, 2019\.
  * [4] N. Goyal and M. Bisi.  Software development efforts prediction using artificial neural network.  IET Software, 10, 04 2016.
  * [5] M. Jørgensen.  What we do and don’t know about software development effort estimation.  IEEE Software, 31(2):37–40, 2014.
  * [6] M. Jørgensen and M. Shepperd.  A systematic review of software development cost estimation studies.  IEEE Transactions on Software Engineering, 33(1):33–53, 2007.
  * [7] B. Kitchenham, S. Pfleeger, B. McColl, and S. Eagan.  An empirical study of maintenance and development estimation accuracy.  Journal of Systems and Software, 64:57–77, 10 2002.
  * [8] C. Martín and A. Abran.  Neural networks for predicting the duration of new software projects.  Journal of Systems and Software, 101:127–135, 03 2015.
  * [9] M. Ochodek, S. Kopczyńska, and M. Staron.  Deep learning model for end-to-end approximation of cosmic functional size based on use-case names.  Information and Software Technology, 123:106310, 2020.
  * [10] P. Phannachitta.  On an optimal analogy-based software effort estimation.  Information and Software Technology, 125:106330, 05 2020.
  * [11] P. Rai, D. K. Verma, and S. Kumar.  A hybrid model for prediction of software effort based on team size.  IET Software, 15(6):365–375, 2021.
  * [12] F. Sarro, R. Moussa, A. Petrozziello, and M. Harman.  Learning from mistakes: Machine learning enhanced human expert effort estimates.  IEEE Transactions on Software Engineering, pages 1–1, 2020.
  * [13] F. Sarro, A. Petrozziello, and M. Harman.  Multi-objective software effort estimation.  In Proceedings of the 38th International Conference on Software Engineering, ICSE ’16, page 619–630, New York, NY, USA, 2016. Association for Computing Machinery.
  * [14] K. Srinivasan and D. Fisher.  Machine learning approaches to estimating software development effort.  IEEE Transactions on Software Engineering, 21(2):126–137, 1995\.
  * [15] V. Tawosi, A. Al-Subaihin, and F. Sarro.  Investigating the effectiveness of clustering for story point estimation.  In 2018 8th International Conference on Computer Science and Information Technology (CSIT). IEEE, 2022.
  * [16] V. Tawosi, A. A. Al-Subaihin, R. Moussa, and F. Sarro.  A versatile dataset of agile open source software projects.  CoRR, abs/2202.00979, 2022.
  * [17] V. Tawosi, R. Moussa, and F. Sarro.  Deep learning for agile effort estimation have we solved the problem yet?  CoRR, abs/2201.05401, 2022.
  * [18] V. Tawosi, F. Sarro, A. Petrozziello, and M. Harman.  Multi-objective software effort estimation: A replication study.  IEEE Transactions on Software Engineering, pages 1–1, 2021.
  * [19] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin.  Attention is all you need, 2017.
  * [20] Y. Yang, X. Xia, D. Lo, and J. Grundy.  A survey on deep learning for software engineering.  ACM Comput. Surv., dec 2021.