---
title: "espadon overview"
author: "Cathy & Jean-Marc Fontbonne"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{espadon overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6, fig.height = 4, out.width = "90%", fig.align = 'center'
)
```
```{r}
par(cex = 0.5, title.cex=0.8)
```
## 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:

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = TRUE}
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.


```{r, echo = TRUE, eval = TRUE}
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)
```
```{r, echo = TRUE, eval = TRUE}
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:

```{r, eval = FALSE}

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.

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
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.

```{r, eval = TRUE}
# 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.

```{r, eval = TRUE}
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.

```{r, eval = TRUE}
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.

```{r, eval = TRUE}
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:

```{r, eval = TRUE}
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)
```
```{r, eval = 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 :

```{r, eval = TRUE}

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)
```

```{r, eval = FALSE}
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.

```{r, eval = TRUE}
# 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).

```{r, eval = TRUE}
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)
head(indices$volume)
```

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.
```{r, eval = TRUE}


# 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))

metrics.le[[1]]
metrics.le[[2]]
metrics.le[[3]]
```

### 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.
