Using constraint structure and an improved object detection network to detect the 12^{th} Vertebra from CT images with a limited field of view for image-guided radiotherapy
Yunhe Xie Kongbin Kang Gregory Sharp David P. Gierga Theodore S. Hong Thomas Bortfeld
cs.CV physics.med-ph


A constraint feature aided 3D detection network

Y. Xie et al.
Massachusetts General Hospital, Department of Radiation Oncology,
Boston, 02114 USA 
 yxie5@mgh.harvard.edu CVCS, Providence,
RI, 02906, USA
 kongbin.kang@gmail.com
Using constraint structure and an improved object detection network
to detect the 12^th Vertebra from CT images with a limited field
of view for image-guided radiotherapy
    Yunhe Xie1 Kongbin Kang2Gregory Sharp1 David P. Gierga1 Theodore S. Hong1 Thomas
Bortfeld1
    December 30, 2023
==========================================================================================================================================================================


Image guidance has been widely used in radiation therapy. Correctly
identifying the bounding box of the anatomical landmarks from limited
field of views is the key to success. In image-guided radiation therapy
(IGRT), the detection of those landmarks like the 12th vertebra (T12)
still requires tedious manual inspections and annotations; and superior-inferior
misalignment to the wrong vertebral body is still relatively common.
It is necessary to develop an automated approach to detect those landmarks
from images. The challenges of training a model to identify T12 vertebrae
automatically mainly are high shape similarity between T12 and neighboring
vertebrae, limited annotated data, and class imbalance. This study
proposed a novel 3D detection network, requiring only a small amount
of training data. Our approach has the following innovations, including
1) the introduction of an auxiliary network to build constraint feature
map for improving the model's generalization, especially when the
constraint structure is easier to be detected than the main one; 2)
an improved detection head and target functions for accurate bounding
box detection; and 3) an improved loss functions to address the high
class imbalance. Our proposed network was trained, validated and tested
on anotated CT images from 55 patients and demonstrated accurate distinguish
T12 vertebra from its neighboring vertebrae of high shape similarity.
Our proposed algorithm yielded the bounding box center and size errors
of 3.98±2.04mm and 16.83±8.34mm, respectively. Our approach
significantly outperformed state-of-the-arts Retina-Net3D in average
precision (AP) at IoU thresholds of 0.5 and 0.35, with AP increasing
from 0 to 95.4 and 0 to 64.7, respectively. In summary,
our approach has a great potential to be integrated into the clinical
workflow to improve the safety of IGRT. We implemented our approach
in Python and Keras and will make it available under the open-source
MIT License after acceptance.




Although there are approaches (e.g. RetinaNet) utilizing pre-trained
2D models (e.g. VGG, ResNet, etc) to address some of those challenges
for 2D applications. However, 2D features cannot provide sufficient
diagnosis information to distinguish the anaotomy of 3D curved
structure to be detected.






§ INTRODUCTION 


Volumetric data (e.g. computed tomography (CT), cone-beam computed
tomography (CBCT), 4D CT, positron emission tomography, magnetic
resonance imaging, etc.) are commonly used in radiation oncology.
Image-guided radiation therapy (IGRT) utilizes image information to
ensure radiation dose to be delivered precisely and effectively, while
reducing side effects. Correctly identifying anatomical landmarks,
like T12 vertebra, is the key to success. In IGRT, the detection of
such landmarks still requires manual inspections and annotations in
a slice-by-slice manner, which can be tedious. More over, superior-inferior
misalignment to the wrong landmark, particularly in the applications
like CBCT which suffers from high level of noises and a limited field
of view (FOV), is still relatively common. It is necessary to develop
an automated approach to detect those landmarks from images.

Although vertebra structures are highly contracted in CT images, detecting
T12 vertebra from IGRT images remains challenge due to the shape similarity
among the neighboring vertebrae. Automated vertebra detection has
been an active research topic for decades. Previous approaches <cit.>utilized
the hand-crafted features to detect vertebrae, which is can be time-consuming,
labor-intensive, and highly dependent on the expert's knowledge and
experience. Moreover, hand-crafted features may not be able to capture
the complex patterns and variations in the data, especially in high-dimensional
and noisy data. Because of the superior performance of Convolutional
Neural Networks (CNN) and their variant, fully convolutional neural
networks (FCN) to learn image features of the structures directly
from data<cit.>,
it is natural to develop CNN or FCN based approaches. Recently, Cheng
et al. <cit.> presented a FCN based method to detect
vertebrae of the full spine. However, their approach focused on full
spine, using a FCN based approach to detect vertebra from limited
view remains an important research topic that needs to be addressed.

For instance, Herring and Dawant <cit.> employed a
surface-based registration method to automatically identified lumbar
vertebrae in CT images with an unknown field of view. Vrtovec et
al. <cit.> presented a method for spine
curve detection in CT and MR images using a polynomial model to described
curved planar reformation (CPR) of the spine column. While Mastmeyer
et al. <cit.> proposed precise segmentation
of lumbar vertebrae in CT images, the method required manually placement
of landmarks in the vertebral bodies. Schmidt et al. <cit.>
proposed probabilistic graphical models to automatically locate the
vertebral column and labelling the intervertebral disks, but the whole
spine is focus. Klinder et al. <cit.> presented
an automatic vertebra detection and segmentation based on the statistic
model of vertebra. Guo et al. <cit.> employed Haar-like
features for vertebra detection from CT scout images. All those approaches
ultilized the hand-crafted features to detect vertebrae, which is
can be time-consuming, labor-intensive, and highly dependent on the
expert's knowledge and experience. Moreover, hand-crafted features
may not be able to capture the complex patterns and variations in
the data, especially in high-dimensional and noisy data. Cheng et
al. <cit.> presented an automatic vertebra detection method
using the a 2D-Dense-UNet for vertebra localization and a 3D-Dense-UNet
for segmentation of the full spine to incrase detection realiability. 


Convolutional Neural Network (CNN) and its variant, fully convolutional
neural network (FCN), have been used in pioneer studies for detecting
disease structures im medical images <cit.>. Examples
of these including colonic polyps in CT images <cit.>,
cerebral micro-bleeds in MRI scans <cit.>, breast
and lung cancer in ultrasound images <cit.>. Recently,
object detection networks, such as <cit.>,
have shown superior performance in computer vision tasks like self-driving
<cit.>, face/ID recognition
<cit.>, Document Analysis <cit.>. These
developments have led technical advancements for applying it to the
medical image analysis<cit.>.
However, there are two main challenges need to be addressed. First,
those computer vision approaches are designed for 2D applications
and lack supports for 3D tomography images that are crucial in
the medical field. Second, those approaches are trained (or based
on pre-trained models) from millions of annotated 2D images, which
are difficult, if not impossible, to obtain in the medical field.
Successful application of FCN based algorithm to medical images requires
significant modified in order to address the challenges specific to
medical images <cit.>,
including scarcity of available data, high cost of data annotation,
high imbalance of classes, and the high memory demands of 3D networks.

 

Although state-of-the-arts detection networks like Retina-Net <cit.>
have achieved impressive results in computer vision, they often struggle
to distinguish the T12 vertebra from the neighboring vertebrae of
high similarity. To address this challenge, we propose a novel approach
that introduces an auxiliary network to generate deep-learning features
of thoracic rib 12 (T12 rib) to assist the detection network to distinguish
the T12 from its neighbors. T12 rib is the last pair of ribs in the
human ribcage and is relatively easy to detect from CT scans. It articulates
with T12 vertebra at the back of spine and forms a constraint to T12
vertebra. In addition to the constraint feature network, our approach
differs from state-of-the-art methods by attaching the detection head
to both encoder and decoder paths of the U-Net to improve the detection
accuracy, rather than relying merely on the decoder path of the U-Net. 

In summary, our approach has several innovations. Firstly, we introduce
an auxiliary network to generate deep-learning features assisting
the detection network in improving the model's generalization, especially
when the constraint structure is easier than the main structure to
detect. Secondly, we introduce an improved detection head and target
functions for accurate bounding box detection. Finally, we propose
an improved loss functions to address the high class imbalance. To
evaluate our proposed approach, we trained the network to detect T12
vertebra from CT images. This particular vertebra, which is the last
thoracic vertebra in the spine of the human body, was chosen due to
its technical difficulty in being distinguished from its neighboring
vertebrae in the limited FOV. Additionally, the T12 vertebra is widely
used as a landmark in IGRT and various other radiation oncology clinical
applications. By focusing on the detection of T12, we aimed to demonstrate
the efficacy of our approach in a challenging context that has practical
implications in clinical settings.




§ NETWORK ARCHITECTURE


As shown in Fig. <ref>, our proposed architecture
consists of three components: a 3D main feature extraction network
(MFEN), a 3D constraint feature extraction network (CFEN), and a 3D
region detection network (RDN). RDN is then attached to the top of
both MFEN and CFEN to output the bounding box field.



 §.§ 3D FENs for structure to be detected and CFEN for constraint structure


Our MFEN and CFEN are designed to leverage the effectiveness of 3D
U-Net in handling limited annotated training datasets, which has been
shown in various medical research studies. To achieve this, we use
a standard 3D U-Net with a 4-layer up-path and down-path as its
backbones. MFEN uses the full down-path and first two layers in up-path,
and CFEN uses only the down-path of U-Net, as shown in the dash box
of Fig. <ref>a. The input to both networks
of MFEN and CFEN is a volumetric image of dimension N_x× N_y× N_z.
The outputs from CFEN and MFEN down-path feed to a down-path fusion
module. The fusion module first resamples the tensors to a final spatial
resolution M_x× M_y× M_z and then concatenate channels
into a final tensor of dimension M_x× M_y× M_z× d.
The outputs from MFEN up-path also feed into a fusion module to match
the spatial resolution and concatenate the features into a tensor
of M_x× M_y× M_z× u. These two feature fusion
modules combines the multiple layer outputs of the different spatial
resolution, allowing the model to learn the optimal way to select
the correct layer from training data, which is particularly useful
to the structures of anisotropic dimensions. In our experiments, we
used values of N=416×216×128, M=52×36×32
for training, with d and u values of 112 and 96, respectively.
The intuition behind this architecture is that CFEN learns the features
of constraint structure and aids MFEN to learn the major features
to achieve a good generalization. 



 §.§ 3D Region Detection Network 


RDN attaches to both down-path feature fusion and up-path feature
modules. Its output is a tensor of size M_x× M_y× M_z×8l,
where l represents the number of structures to be detected and
8 represents the number of the bounding box parameters (two for probability,
three for locations and three for box sizes). The output tensor represents
a field of box parameters, including the inside-box probability, box
center probability, box location and box size, defined on the coarse-level
spatial grids. In this paper, we use the bounding box aligned to the
image coordinate axis for simplicity. To recover the true box center
from the coarse spatial grids, it is necessary to estimate the the
true center location (or offset). While our network is designed to
detect multiple anatomical structures simultaneously, we only demonstrated
one structure in our experiments.



a) MFEN and CFEN share the same architecture of 3D U-Net. RDN is
attached on two scale fusion modules RDN is attached to two module
of fusion of different scales. b) In RDN network a 8l-vector is
shown for simplicity at the output end of network, its actual output
is a M_x× M_y× M_z×8l tensor, where l is
the number of structures to be detected. The horizontal bars represent
the output vectors from convolution layers with 3×3×3-kernels.
The blue circles, except for p, represent output nodes from convolution
layers with 1×1×1-kernel. The blue circle of center probability,
p_c, is the output from multiplication layer with inputs of p_x,
p_y and p_z.


which consists of three parts: MFEN, CFEN and RDN.


Each blue 3D box represents a 3D U-Net block that contains two
convolution layers, a maximum pooling layer for down-path or a up-sampling
layer for up-path, and a dropout layer.


It's worth noting that FEN is different from FPN <cit.>
since there are no multiple prediction layer in the skip connections.

The architecture of the RDN is as shown in Fig. <ref>b.
The input tensors from both fusion module are firstly convolved with
three small 3×3×3 kernels to generate three low dimensional
feature maps, each consisting of 64-vector at a grid point (x,y,z),
represented by narrow bars in Fig. <ref>b).
Next, the resulting feature vectors are convolved by 8l 1×1×1-kernels
to generate all the x-, y- and z-components of the inside-box
probability p_b, box center probability p_c, location t
and box size s. Finally, the joint probability, p, is a multiplication
from its three axis-components. This architecture was implemented
naturally using 3D convolution and multiplication layers. To mitigate
the effects of image scanning intensity variations, a sigmoid function
is used as the activation function in the convolution layer with 3×3×3-kernels.
Our RDN differs from its counterpart in <cit.>
in that it decomposes the 3D box parameter regressions to multiple
independent regressions on each axis. This approach reduces the number
of weight parameters required in the network, thus helping to mitigate
overfitting with limited training data.



§ TARGET LABELING




 §.§ Augmentation


Besides the network architecture, data augmentation is another important
step to train deep networks with limited training data sets. A properly
designed augmentation teaches the network to focus on robust features
for a good generalization. In our study, each input image and its
corresponding annotated masks were augmented twenty-five times using
a 2D elastic deformation algorithm based on <cit.>.
In each augmentation, a 3×3 grid of random displacements were
drawn from a Gaussian distribution (μ=(0,0) pixels and σ=(10,10)
pixels). Then the displacements were interpolated to the pixel level.
All the 2D slices of the volumetric image were deformed using the
same displacements and a spline interpolation. 



 §.§ RDN and FEN Output Target


CFEN was trained using standard U-Net setup and only the down path
was kept after training. We also devised a customized ground truth
approach to train RDN and FEN networks. Although standard target generating
methods, such as the intersection over union (IoU), are commonly used
in computer vision <cit.>, we
found that they were not effective in our experiments due to the optimization
plateau issues<cit.>. Instead, we developed a
hand-crafted probability function that has a global peak value at
the box center to overcome these problems.

To create the ground truth, we started by obtaining annotated images
and use the structure contours drawn by experienced radiation oncologist
or physicists to compute the target box parameters, including the
box center locations and box sizes. We then used these parameters
to calculate an 8-vector target, (p̂_c,p̂_bt̂_x,t̂_y,t̂_z,ŝ_x,ŝ_y,ŝ_z),
for every spatial grid location (x,y,z). The target was designed
as follows: when inside a target box b, the inside-box probability
target p̂_b=1 and 0 otherwise. The center probability p_c
is sigmoid decay from the box center, (x_bc,y_bc,z_bc), along
the three axes; the center offset targets, t̂_b, linearly
increases from the box center; and the box size target, ŝ_b,
is a constant vector. When outside the box b, p̂_b, t̂_b,
and ŝ_b are set to zero. This approach teaches the network
to focus on the features of the structure and ignore those outside
it. All those targets are scaled by the down sampling rate, (r_x,r_y,r_z),
caused by pooling operators inside FEN to match the coarse spatial
resolution of RDN output. 


    p̂_c(x,y,z)     =    σ(4[1-|Δ x|/r_x])σ(4[1-|Δ y|/r_y])σ(4[1-|Δ z|/r_z])
    
    p_b(x,y,z)     =    
    1     (x,y,z)∈ b
    
    0     otherwise
    
    t̂_b(x,y,z)     =    (Δ x/r_x,Δ y/r_y,Δ z/r_z)     (x,y,z)∈ b
    (0,0,0)     otherwise
    
    ŝ_b(x,y,z)     =    (w_b/r_x,l_b/r_y,h_b/r_z)     (x,y,z)∈ b
    (0,0,0)     otherwise

where σ(·) is the sigmoid function. Δ x, Δ y
and Δ z are x-x_bc, y-y_bc and z-z_bc, respectively,
and w_b, l_b, and h_b are the box width, length and
height, respectively. 



§ RDN LOSS FUNCTIONS


We assign a total loss to measure the discrepancy between the RDN
output tensor and the target tensor. The total loss, L_total,
consists of the losses of the probability p, the location offset
t, and the box size s. Here p, t, and s are the estimated
probability, center offset and box size, and p̂, t̂
and ŝ are their ground-truth counterparts. The loss functions
are defined as

    L_c    =    ∑_b∑_i,j,k1_(i,j,k)∈ b‖ t_b(x_i,y_j,z_k)-t̂_b(x_i,y_j,z_k)‖/Σ_i,j,k1_(i,j,k)∈ b
    
    L_s    =    ∑_b∑_i,j,k1_(i,j,k)∈ b‖ s_b(x_i,y_j,z_k)-ŝ_b(x_i,y_j,z_k)‖/Σ_i,j,k1_(i,j,k)∈ b
    
    L_p    =     1-2∑_i,j,k∑_bp_b(x_i,y_j,z_k)p̂_b(x_i,y_j,z_k)/∑_i,j,k∑_b(p_b^2(x_i,y_j,z_k)+p̂_b^2(x_i,y_j,z_k))+∑_b∑_i,j,k1_(i,j,k)∈ b‖ p_c-p̂_c‖/Σ_i,j,k1_(i,j,k)∈ b

where i, j and k are spatial grid index of the RDN output
tensor and ‖ .‖ is the l_1-norm. The total
loss is defined as L_total=L_p+L_c+L_s. The summation over
all the spatial grid index (i,j,k) is needed in order to measure
the overall discrepancy of the tensor.

With the above definitions, RDN and FEN networks can be trained jointly
end-to-end. The trained network then computes the output tensor when
a tomography image is given. The detection was performed by finding
the maximum probability inside the output tensor and extracting its
associated box parameters. 



§ IMPLEMENTATION DETAILS 


We implemented our method using Python (ver 3.7.10) with Keras (ver
2.4.3) and TensorFlow (ver 2.4.1) as back-end. We cropped all training
images to 416x288x128 to fit into 24GB GPU memory, and used (2,2,2)
pooling size for all pooling layers. Our model was trained end-to-end
using Adam optimizer with a learning rate of 10^-4. It took about
16 hours to train the model on a single NVidia Titan RTX GPU card
with CUDA driver ver 11.4. During the prediction phase, our trained
model generates proposal regions of the top five scores and no structure
is detected if the maximum probability is smaller than 0.5. The training
was run for 150 epochs and the model with the minimum fitting loss
was kept.



§ EXPERIMENT RESULTS


This study analyzed 55 CT scans from patients who received proton
or photon treatments for various sites including esophagus, lung,
liver, and pancreas at  Masschusset General Hospital from 2012 to
2020. The patients were aged between 34 and 87 years, with 30 males
and 25 females. The CT images were acquired using a GE force CT scanner
at 140 kVp with a field of view ranging from 41 cm to 70 cm. The pixel
sizes ranged from 0.98mm to 1.37mm, and the slice thickness was 2.5mm.
Patients were immobilized on a flat couch in the same posture as for
the radiation treatments during the scanning, and T12 vertebra and
ribs were contoured by experienced radiation oncologists and board-certified
medical physicists. This study was approved by the research IRB.

Those datasets were randomly plit into training, validation and test
datasets consisting of 36, 5 and 14 patients, respectively.
To ensure performance comparability, we implement a 3D RetinaNet with
the same backbone as our proposed algorithm as the baseline. For the
purpose of IGRT safety check, we measured the distance error and the
size error between the ground truth box center and the model predicted
center. Using our trained model, the detection took approximately
3.30 seconds in a CT dataset of 416×288×128 voxels. 


On the test CT images, our proposed approach achieved a mean error
of 0.58, 0.74, and 1.00mm in the left-right, anterior-posterior
and superior-inferior directions, respectively, for the detected center.
The corresponding standard deviations of the center location error
were 2.90, 2.78, and 3.36mm, respectively. For the box sizes,
our approach yielded a mean error of 0.64, -2.38, and -1.40mm
with std. of 12.52, 17.76, and 13.23mm on the direction left-right,
anterior-posterior and superior-inferior directions, respectively.
We summarized these results in Table <ref>.
Specifically, we accurately identified the center of T12s with a mean
overall detection error of 3.98±2.04mm, and the superior-inferior
mean detection error is 1.00mm. Given a typical slice thickness
is 2.5mm, our approach's accuracy is sufficient enough to be integrated
into the IGRT workflow.

We conducted ablation studies to compare our approach to the baseline
RetinaNet3D and investigated the detection performance changes with
each proposed improvement. In the study, the matching IoU threshold
for positive proposals was selected as 0.35 and 0.5, which is the
same as the popular choice of 2D <cit.> or an extrapolating
the 3D IoU threshold<cit.> to 3D by (√(0.5))^3=0.35,
respectively. In practice, we found that the detected bounding box
with threshold of 0.35 already closely matched the ground-truth
. As shown in Table <ref>, our approach significantly outperformed
the baseline RetinaNet3D in T12 detection AP. Introducing the constraint
feature map improved performance but was still insufficient. However,
incorporating rib constraints and improved RDN significantly enhanced
performance. In the experiments, we found that the RetinaNet3D had
poor T12 detection performance due to detecting the wrong center or
the detection size way off. With the introduction of constraint features
and an improved RDN, our approach consistently and accurately identified
the bounding box, as demonstrated in Fig. <ref>. 

Fig. <ref> shows an example of T12 was correctly detected
in CT by our proposed approach, even this patient's T12 vertebra was
cracked due to compression.










§ CONCLUSION


Our study has presented a novel 3D detection network that incorporates
a constraint feature network, an improved detection head and target
functions, and an improved loss function to accurately detect bounding
boxes of anatomical structures of high similarities in 3D volumetric
data. Our approach outperforms existing methods, with significant
improvement in average precision for different IoU thresholds. With
the ability to train on a small amount of data, our approach has a
great potential to be integrated into the clinical workflow for more
accurate, efficient and safe image guided therapy. Future research
could explore the generalization of the proposed approach to other
anatomical structures and imaging modalities, as well as evaluate
its performance in a larger and more diverse dataset. Furthermore,
we plan to investigate 3D data augmentation techniques to simulate
patients' movements, bending, stretching, and compression etc.



§ ACKNOWLEDGMENT


We would also like to thank Dr. Yuling Guo and Dr. Yi Dong for the
technical discussions. 

splncs04


