espadon overview

Easy Study of Patient DICOM Data in Oncology

Espadon is a package simplifying the use and exploitation of DICOM data in radiotherapy. Espadon works in a native way (user friendly):

on images (CT, MR, PT) and fully manages changes of frame of reference (REG)
on dosimetry (rt-dose, rt-plan)
on structures (rt-struct)

It is also able to use any DICOM format file, as long as the user knows where to look for the information.In addition to the simplified use of the above-mentioned files, espadon contains many functions that allow the user to rework documents to produce new features that can be integrated into machine learning.

Among other features, espadon is suitable for pseudonymizing DICOM, decoding them, visualizing them and producing dosimetry data in relation to clinical data. The package can export espadon object in DICOM files.

Also see https://espadon.cnrs.fr.

Input data

The espadon package acts directly on the DICOM files or the folders containing them, thanks to espadon functions whose names contain the word “dicom”. For example, to load an imaging volume from a CT or MR, simply execute the following instruction:

library(espadon)
vol <- load.obj.from.dicom (image_path)

image_path represents either the image folder path, or a vector of all the image DICOM file path.

To speed up the loading of data when studying a patient’s images in depth, espadon offers the dicom.to.Rdcm.converter data pre-formatting function. This function creates files that are smaller than the original files, without loss of information, and prepares the data for use by espadon’s analysis and display functions. The new data format is a Rdcm format.

library(espadon)
dicom.to.Rdcm.converter("D:/dcm/patient001", "D:/Rdcm/patient001“)

DICOM Handling

The espadon package allows direct analysis of DICOM files. The following instructions, for example, create a table (R dataframe), and thus to have a synthetic view of the DICOM file content.

library(espadon)
dcm.filename <- file.choose ()
df <- dicom.parser (dicom.raw.data.loader (dcm.filename), as.txt = TRUE, 
                    try.parse = TRUE)

A file in .xlsx format can also be created with the xlsx.from.dcm function.

Thanks to the dicom.set.tag.value function, the content of the various tags can be modified or deleted, as far as they are in ascii format.

The dicom.raw.data.anonymizer function provides pseudonymization by shifting dates, changing the patient’s PIN and name, removing all patient data except age, weight, and height, and also removing any other name, phone, operator, service, and undocumented tag content.

Patient overview

Espadon provides a unified view of the different imaging objects, by assigning them a name that takes into account their relationship. The following functions (depending on whether the patient’s directory contains DICOM or .Rdcm files) create these links.

The load.patient.from.dicom and load.patient.from.Rdcm functions (depending on whether the patient’s folder contains DICOM or .Rdcm files) create an R-list, in which all the DICOM objects of the patient are identified. In this patient list, there is, among others, the element T.MAT, computed from the DICOM reg objects. This one provides all the necessary information to overlay two imaging objects acquired on different machines, or at different times. The display.obj.links function allows a nice and convenient display of these links.

library(espadon)
pat.dir <- choose.dir () # patient folder containing mr, ct, rt-struct, rt-dose, rt-plan and reg…
pat <- load.patient.from.dicom (pat.dir)
display.obj.links (pat)

When the patient object list is loaded into memory, the load.obj.data function loads the full object information, i.e. voxels content for MR, CT or PT, rt-dose or contour coordinates for rt-struct.

S <- load.obj.data (pat$rtstruct[[1]])
CT <- load.obj.data (pat$ct[[1]]) 
MR <- load.obj.data (pat$mr[[1]])
D <- load.obj.data (pat$rtdose[[1]])

A toy patient has been created for testing purposes, and can be loaded with the instruction : toy.load.patient. The examples below will use the toy patient, but you can also use the codes with your own data.

library (espadon)
pat <- toy.load.patient(modality = c("ct", "mr", "rtdose", "rtstruct"),
                        dxyz = c(2, 2, 2),beam.nb = 3)
S <- pat$rtstruct[[1]]
CT <- pat$ct[[1]]
MR <- pat$mr[[1]]
D <- pat$rtdose[[1]]
display.obj.links (pat)

2D Display

Visualization of imaging volumes is essential for controlling the loaded DICOM objects. The espadon package allows this display, with the possibility of using the color palettes of your choice.

The display.kplane function displays a section of the imaging volume and contours in the requested reference frame. The display.plane function displays transverse, frontal or sagittal views in the patient’s frame of reference. It allows the superposition of images and contours. The espadon plot function displays either a section of the imaging volume, or contours in the reference frame of the displayed object, and is more flexible than display.plane.

par(mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.9)
display.kplane (CT, k =10, col = pal.RVV(1000), 
                flop = TRUE, interpolate = TRUE)


ptv.center <- as.numeric(S$roi.info[S$roi.info$name == "ptv", c("Gx","Gy","Gz")])


display.plane(CT, D, S, main ="transverse", view.coord = ptv.center[3], 
              legend.shift = -100)

display.plane(CT, D, S, main ="frontal", view.type = "front", 
              view.coord = ptv.center[2], legend.shift = -100)

display.plane(CT, D, S, main ="sagittal", view.type = "sagi",
              view.coord = ptv.center[1], legend.shift = -100)

par(mar = c(4, 4, 2, 1), cex.axis = 0.8, cex.lab = 0.8, cex.main = 0.9)
plot(CT,view.type = "yz", view.coord = ptv.center, col = pal.RVV(100), las = 1)
plot(D, add = TRUE, view.type = "yz", view.coord = ptv.center, 
     col= pal.rainbow((100)))
Scut <- plot(S, add = TRUE, view.type = "yz", view.coord = ptv.center, lwd = 2)
legend("topleft",legend = Scut$roi.info$name, text.col = Scut$roi.info$color,
       cex = 0.8, bty="n")

3D Display

The espadon package has several 3D visualization features:

the display.3D.contour function displays the contours of the selected region of interest as described in the rt-struct.
the display.3D.stack function displays selected slices of the imaging volume.
the display.3D.sections function displays the frontal, sagittal and transverse view at a point. These last 2 functions can use the same color palettes as the 2D visualization functions, except that a color can only be totally transparent or visible.

Thanks to the internal management of geometrical frames of reference, the different displays can be mixed, as shown in the example below:


library(rgl)
bg3d ("black")
par3d (userMatrix = matrix (c(0, 0, -1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1), ncol = 4), 
       windowRect = c(0, 50, 300, 300), zoom = 0.7)

display.3D.contour (S, roi.sname = "eye", display.ref = CT$ref.pseudo, 
                    T.MAT = pat$T.MAT)
display.3D.stack (MR, k.idx = round(seq(3, length(MR$k.idx)-3,length.out=10)), 
                  display.ref = CT$ref.pseudo, 
                  T.MAT = pat$T.MAT,
                  ktext = FALSE)
display.3D.sections (CT, cross.pt = c(0, 110, 0), 
                     col = pal.RVV (200, alpha = c(rep (0, 90), rep (1, 110))),
                     breaks = seq(-1000, 1000, length.out = 201), 
                     border.col = "#844A39")

rglwidget()
close3d()

Geometry tools

Several espadon functions have been created to transform imaging volumes.

Functions such as add.margin, nesting.cube, nesting.roi, nesting.bin increase or decrease the number of voxels in the volume.

nesting.MR <- nesting.roi (MR, S, roi.sname = "brain", T.MAT = pat$T.MAT, 
                           vol.restrict = TRUE, xyz.margin = c(2, 2, 2))
display.3D.stack (nesting.MR, ktext = FALSE, line.lwd= 2)
display.3D.contour (S, roi.sname = "brain", display.ref = MR$ref.pseudo, 
                    T.MAT = pat$T.MAT)

The vol.regrid function resamples the volume to the grid of another imaging volume.

MR.on.CT <- vol.regrid (MR, CT, T.MAT=pat$T.MAT, interpolate = TRUE)
display.3D.stack (MR.on.CT, display.ref = CT$ref.pseudo , T.MAT = pat$T.MAT, 
                  ktext = FALSE, line.lwd = 2, line.col = "#844A39")

More generally the get.plane function creates a cut of the volume in any direction. The get.direction function creates the direction vector from 3 points in space.

# Points determination : centers of the eyes and the chiasm found in rtstruct data
roi.idx <- select.names (S$roi.info$roi.pseudo, 
                         roi.name = c("left eye" ,"right eye","ptv"))
pt <- S$roi.info[roi.idx, c("Gx", "Gy", "Gz")]
origin.pt <- as.numeric (apply (pt, 2, mean))

# Plane creation in the CT frame of reference 
new.plane <- get.plane (CT, origin = origin.pt, 
                        plane.orientation =  orientation.create (pt), 
                        interpolate = TRUE)

# Display of the new plane, in the cut plane frame of reference
tmat <- ref.cutplane.add (new.plane, origin = origin.pt, 
                          ref.cutplane = "ref.plane", T.MAT = pat$T.MAT)
plan <- vol.in.new.ref (new.plane, "ref.plane", tmat)
display.plane (plan, struct= S, roi.idx = roi.idx, T.MAT = tmat, 
               view.coord = plan$xyz0[3],
               bottom.col = pal.RVV(200), 
               bottom.breaks = seq(-1000, 1000, length.out = 201), 
               main = "cut plane", bg = "#379DA2", sat.transp = T, 
               legend.plot = FALSE, interpolate = TRUE)

Binary and weight objects

Binary objects are volumes whose voxel values are either TRUE, FALSE or NA (i.e. undetermined). They act as a mask, and are used to select only the desired part of a volume. The example below show how to select the brain voxels described in the rtstruct object.

nesting.MR <- nesting.roi (MR, S, roi.sname = "patient", T.MAT = pat$T.MAT, 
                           vol.restrict = TRUE, 
                           xyz.margin = c(1, 1, 1)) # to reduce the computation time 
bin.patient <- bin.from.roi (nesting.MR, struct = S, roi.sname = "patient", 
                           T.MAT = pat$T.MAT)
MR.patient <- vol.from.bin (nesting.MR, bin.patient)
plot(MR.patient)

Weight objects are volumes whose voxel values represent the proportion of the selection area contained within the voxels. As a result, voxel values range from 0 to 1.

bin.patient_ <- bin.from.roi (nesting.MR, struct = S, roi.sname = "patient", 
                           T.MAT = pat$T.MAT)
MR.patient_ <- vol.from.bin (nesting.MR, bin.patient_)
par(mfrow=c(1,2))
plot(bin.patient, main = "binary", cut.interpolate  = FALSE)
plot(bin.patient_, main="weight", cut.interpolate  = FALSE)

par(mfrow=c(1,1))

The bin.from.vol function creates binary volumes by selecting voxels that are inside or outside an interval [min, max]. The example below selects voxels with a value greater than 150.

bin.min.150 <- bin.from.vol (MR.patient, min = 70)

display.plane (bin.min.150, view.coord = 0, view.type = "sagi",
               main = "bin.min.150", bg = "#379DA2", legend.plot = FALSE, 
               sat.transp = TRUE, interpolate = TRUE)

Other image processing functions are available:

bin.erosion decreases the volume of a radius r
bin.dilation increases the volume of a radius r
bin.closing is usefull for filling holes that are smaller than the radius and merging two shapes close to each other.
bin.opening is usefull for removing volumes that are smaller than the radius and smoothing shapes.
bin.clustering groups and labels TRUE voxels that have a 6-connectivity (i.e. sharing a common side).

These functions can be used to separate multiple volumes in the binary selection, as shown below:

bin.smooth <- bin.closing (bin.min.150, radius = 2 * max(bin.min.150$dxyz))
bin.cluster <- bin.clustering (bin.smooth)

n=5
col <- c ("#00000000", rainbow (n))
breaks <- seq (-0.5, n + 0.5, length.out = n+2)
display.plane (MR, top = bin.cluster, main = "After clustering", 
               view.coord = 0, view.type = "sagi",
               top.col = col, top.breaks = breaks, 
               interpolate = FALSE)

# Creation of new patient contours
new.bin.patient <- bin.from.vol (bin.cluster, min = 1, max = 1)
S.patient <- struct.from.bin (new.bin.patient , roi.name = "patient", 
                                roi.color = "#FF0000" )

# Display of the first 3 largest sub-volumes, and ventricle contours
library (rgl)
bg3d ("black")
par3d (userMatrix = matrix (c(0, 0, -1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1), 
                            ncol = 4), 
       windowRect = c(0, 50, 300, 300), zoom = 0.7)
display.3D.stack(bin.cluster, k.idx = bin.cluster$k.idx, 
                 col = c ("#00000000", rainbow (3, alpha = 1)), 
                 breaks = seq (0, 3, length.out = 5), border = FALSE, 
                 ktext = FALSE)
play3d (spin3d (rpm = 4), duration = 15)

bg3d ("black")
par3d (userMatrix = matrix (c(0, 0, -1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1), 
                            ncol = 4), 
       windowRect = c(0, 50, 300, 300), zoom = 1)
display.3D.contour (S.patient, roi.col="yellow")
play3d (spin3d (rpm = 4), duration = 15)

Meshes

The espadon package has the mesh.from.bin function to create a triangle mesh from a binary volume and the display.3D.mesh function to display it in any frame of reference.

In the example below, the patient folder contains a CT-scan files and a DICOM rtstruct file in which lung and heart are outlined :


bin.leye <- bin.sum (bin.from.roi (pat$ct[[1]], struct = pat$rtstruct[[1]], 
                                   roi.name = "left eye"),
                     bin.from.roi (pat$ct[[1]], struct = pat$rtstruct[[1]], 
                                   roi.name = "left optical nerve"))

bin.reye <- bin.sum (bin.from.roi (pat$ct[[1]], struct = pat$rtstruct[[1]], 
                                   roi.name = "right eye"),
                     bin.from.roi (pat$ct[[1]], struct = pat$rtstruct[[1]], 
                                   roi.name = "right optical nerve"))

mesh.leye <- mesh.from.bin (bin.leye)
mesh.reye <- mesh.from.bin (bin.reye)

display.3D.mesh (mesh.eye, color = "burlywood2", specular = "black", 
                 alpha = 0.8)
display.3D.mesh (mesh.proc, color = "red", specular = "black", 
                 alpha = 1)

Sometimes, the mesh contains holes. The mesh.repair function can be used to repair the mesh by removing the degenerated triangles.

Histograms

Histograms computing is a basic function proposed by R. The espadon package adds the feature to compute these differential histograms over an imaging volume such as CT, MR, PT, rt-dose, and within a predefined area:

either by a region of interest (RoI) of an rt-struct file : histo.from.roi function
or by a selection by binary volume : histo.from.bin function.

To be able to simulate random shifts of the selection, or variations of contours, the histo.from.roi function includes Monte Carlo calculations, according to a normal distribution. If the imaging volume is a dose distribution (from an rtdose file), it is then easy to calculate a cumulative histogram, called in this case dose-volume histograms or DVH, using the histo.DVH function.

The functions display.dV_dx, display.DVH, display.DVH.pc allow to display pretty graphs, including the quantile zones 100%, 95% and 50% of the Monte-Carlo data.

# resample D on the cut planes on which the chiasm has been outlined
D.on.CT <- vol.regrid (D, CT)

# Calculation of histogram with simulation of random displacements according to 
# a normal distribution of standard deviation 1 mm in all directions:
h <- histo.from.roi (D.on.CT, S, 
                     roi.name="ptv", 
                     breaks = seq(0, 80, 0.1),
                     alias = "ptv", MC = 10)

DVH <- histo.DVH (h, alias = "PTV")

display.DVH (DVH, MC.plot = TRUE, col = "#ff0000", lwd = 2)

Radiotherapy indices

The espadon package offers the possibility of computing many radiotherapy indices. After specifying the target volumes to be treated (if available), the healthy volumes to be protected (if available), the dose distribution from rt-dose file, the functions rt.indices.from.roi and rt.indices.from.bin can calculate, on demand :

dosimetry indices : minimum, maximum, average dose, standard deviation, requested D.x% and requested D.xcc.
volume indices : total volume, surface area and requested V.xGy or V.x%.
conformity indices : prescription isodose target volume (PITV), prescription dose spillage (PDS), conformity index (CI), conformation number (CN), new conformity index (NCI), Dice similarity coefficient (DSC), conformity index based on distance (CIdistance), conformity distance index (CDI), triple point conformity scale (CS3), underdosed lesion factor (ULF), overdosed healthy tissues factor (OHTF), geometric conformity index (gCI), COIN , critical organ scoring index (COSI), generalized COSI (gCOSI).
homogeneity indices from RTOG or ICRU.
gradient indices : gradient index based on volumes ratio, modified gradient index (mGI).

indices <- rt.indices.from.roi (D.on.CT, S, target.roi.name = "ptv", 
                     healthy.roi.sname = c("eye", 'proc', 'ghost'),
                     presc.dose = round(0.9 * D$max.pixel),
                     conformity.indice = "", homogeneity.indices="",
                     gradient.indices = "",
                     V.xpc = c(5, 50, 75), DVH = FALSE)
head(indices$dosimetry)
#>                            D.min D.max D.mean   STD
#> ptv                       46.422    52 49.983 0.041
#> left eye                   0.000     0  0.000 0.000
#> right eye                  0.000     0  0.000 0.000
#> labyrinth processing unit  0.000     0  0.000 0.000
#> ghost container            0.000     0  0.000 0.000
head(indices$volume)
#>                             V_tot      area V_47Gy    V_5%  V_50% V_75%
#> isodose                        NA        NA  5.656 292.648 56.032 7.624
#> ptv                         4.160  11.92563  4.159   4.160  4.160 4.160
#> left eye                   23.128  41.49084  0.000   0.000  0.000 0.000
#> right eye                  23.132  41.49890  0.000   0.000  0.000 0.000
#> labyrinth processing unit  16.108  69.00278  0.000   0.000  0.000 0.000
#> ghost container           162.416 162.98173  0.000   0.000  0.000 0.000

The espadon package can also compute the Gamma index and the Chi index in 2D or 3D, using the rt.gamma.index and rt.chi.index.

Spatial similarity metrics

The espadon package calculates several spatial similarity metrics, such as :

the volumetric DICE similarity defined by DSC = 2 (VA ∩ VB) / (VA+VB)
the Dice-jaccard coefficient defined by DJC = (VA ∩ VB) / (VA ∪ VB)
the mean distance to conformity MDC, over-contouring mean distance over.MDC and under-contouring mean distance under.MDC, defined by Jena et al [1]
the maximum , mean and quantiles of Hausdorff distances
the surface DICE metric defined by Nikolov et al [2]

Metrics using the volumes of regions of interest (ROIs) are calculated using the sp.similarity.from.bin function, and those using surfaces are calculated using the sp.similarity.from.mesh function.

Here is the calculation of the distances between the left ocular system and the brain.



# The binary object for brain is :
bin.brain <- bin.from.roi (CT, S, roi.sname = "brain", modality ="weight")

# The mesh objects for brain:
# FYI: the smooth.iteration argument avoids staircase meshing, due to z-axis binning.
#      As a result, ROI contours are also smoothed in XY.
mesh.brain <- mesh.from.bin (bin.brain, smooth.iteration = 10, alias = "brain")


# spatial similarity between lest ocular system and brain
metrics.le <- c(list(sp.similarity.from.bin (bin.leye , bin.brain)),
                sp.similarity.from.mesh (mesh.leye , mesh.brain, 
                                         hausdorff.quantile = c (0.5, 0.95),
                                         surface.tol = 1:3))
#> Warning in sp.similarity.from.bin(bin.leye, bin.brain): both volumes must be of
#> binary modality.

metrics.le[[1]]
#> NULL
metrics.le[[2]]
#>     label        HD
#> 1  HD.max 103.86927
#> 2 HD.mean  49.04248
#> 3   HD.50  47.97817
#> 4   HD.95  90.58114
metrics.le[[3]]
#>   tol       sDSC     sAPL
#> 1   1 0.00938282 8661.953
#> 2   2 0.01974961 8501.652
#> 3   3 0.02752938 8335.481

References

[1] Jena R, et al. (2010). “A novel algorithm for the morphometric assessment of radiotherapy treatment planning volumes.”Br J Radiol., 83(985), 44-51.doi:10.1259/bjr/27674581.

[2] Nikolov S, et al. (2018). “Deep learning to achieve clinically applicable segmentation of head and neck anatomy for radiotherapy.” ArXiv,abs/1809.04430.