\documentclass[11pt,%
parskip=half,%
paper=a4,%
headings=small,%
DIV15]{scrartcl}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[final,babel,activate=TRUE]{microtype}
%%\usepackage[sc]{mathpazo}
\usepackage{lmodern}
\usepackage{upquote}
%\usepackage{geometry}
%\geometry{verbose,tmargin=2.2cm,bmargin=2.2cm,lmargin=2.2cm,rmargin=2.2cm}
\usepackage{url}
%\usepackage{amsmath}
%\usepackage{bm}
%\usepackage{authblk}
\usepackage{graphicx}
\usepackage[pdftex,%
unicode=true,%
pdfusetitle,%
bookmarks=true,%
plainpages=false,%
colorlinks=true,%
linkcolor=blue,%
citecolor=blue,%
filecolor=blue,%
urlcolor=blue,%
%bookmarksopen=true,%
%bookmarksopenlevel=2,%
breaklinks=false,%
%pdfborder={0 0 1},%
pdfpagelabels=true,%
backref=false,%
pdftitle={Analyzing dose-volume histograms using DVHmetrics for R},%
pdfauthor={Daniel Wollschlaeger, Heiko Karle}]{hyperref}
\usepackage{breakurl}
\usepackage{apacite}                     % after hyperref

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

\ifdefined\hlstd
\renewcommand{\hlstd}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\fi

\ifdefined\hlcom
\renewcommand{\hlcom}[1]{\textcolor[rgb]{0.5,0.4,0.5}{#1}}%
\fi

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

%\VignetteIndexEntry{Analyzing dose-volume histograms using DVHmetrics}
%\VignetteDepends{DVHmetrics}
%\VignetteKeywords{DVHmetrics}
%\VignettePackage{DVHmetrics}
%\VignetteEngine{knitr::knitr}
%%%\VignetteEngine{knitr::rmarkdown}
%%%%\SweaveOpts{engine=R}

\begin{document}

\title{Analyzing dose-volume histograms using \texttt{DVHmetrics} for \textsf{R}}
\author{Daniel Wollschlaeger\\
        \url{wollschlaeger@uni-mainz.de}
        \and
        Heiko Karle\\
        \url{karle@uni-mainz.de}}
\date{University Medical Center Mainz, Germany\\\today}

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

<<setup, echo=FALSE, include=FALSE, cache=FALSE>>=
# set global chunk options
knitr::opts_chunk$set(fig.align='center', fig.show='hold')
knitr::opts_chunk$set(tidy=FALSE, message=FALSE, warning=FALSE)
options(replace.assign=TRUE, useFancyQuotes=FALSE, show.signif.stars=FALSE, digits=4, width=70)
@

\maketitle
\tableofcontents

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{Introduction}
\label{sec:introduction}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

\texttt{DVHmetrics} is an add-on package for the free statistical environment \textsf{R}\footnote{A free short introduction to \textsf{R} can be found at \url{https://www.statmethods.net/}.} \cite{RDevelopmentCoreTeam2008c} with applications in radiation oncology. It provides functionality to read dose-volume-histogram (DVH) text files, to calculate DVH metrics, and to plot DVHs. In addition, it checks and visualizes quality assurance constraints for the DVH.\footnote{For a solution that also reads files in DICOM-RT format, see the \texttt{RadOnc} package for \textsf{R} \cite{Thompson2014}.}

To install \texttt{DVHmetrics}, you need a current version of \textsf{R} and be online. Preferably, a free development environment like \textsf{R}Studio \cite{Allaire2011} should be used.

<<cIntro, eval=FALSE>>=
# install DVHmetrics from the CRAN online package repository
install.packages("DVHmetrics")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{Interfaces}
\label{sec:interfaces}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

\texttt{DVHmetrics} provides two interfaces geared towards users with different levels of familiarity with \textsf{R}: The regular command line functions and a built-in web application.

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{\textsf{R} command line interface}
\label{sec:Rinter}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Users familiar with \textsf{R} can use the \texttt{DVHmetrics} package functions from the \textsf{R} command line. This facilitates statistical post-processing of results with the full capabilities of \textsf{R}. After installing \texttt{DVHmetrics}, you should be able to run (function \texttt{getMetric()} is explained in section \ref{sec:getMetric}):

<<cCmdline>>=
## load DVHmetrics package - required for all following tasks
library(DVHmetrics, verbose=FALSE)

## calculate a DVH metric for built-in data
getMetric(dataMZ, metric="DMEAN", structure="HEART")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Web-based graphical user interface}
\label{sec:webApp}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

For users who are unfamiliar with \textsf{R}, \texttt{DVHmetrics} includes a \texttt{shiny}-based web application \cite{RStudioShiny2014} running locally that eliminates the need to use \textsf{R} syntax. For information on how to use this app, see the documentation by running this from the command line:

<<cWebApp, eval=FALSE>>=
vignette("DVHshiny")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{Read DVH text data}
\label{sec:readData}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

To import DVH data into \textsf{R}, it should be exported as a plain text file from Varian Eclipse$^{\mathrm{TM}}$ (Versions 10--15), CadPlan$^{\mathrm{TM}}$, Pinnacle$^{3\,\mathrm{TM}}$ (version 9\footnote{Pinnacle$^{3}$ files have to be exported such that information from one patient is contained in one directory. The directory layout and required files are explained in \texttt{help(readDVH)}.}), Oncentra MasterPlan$^{\mathrm{TM}}$ (version 4.3), Elekta Monaco$^{\mathrm{TM}}$ (version 5), TomoTheray HiArt$^{\mathrm{TM}}$, RaySearch Labs RayStation$^{\mathrm{TM}}$, Medcom ProSoma$^{\mathrm{TM}}$, or from PRIMO (version 0.3.1.1558). DVH files from different TPSs can be combined into one set of DVHs. Cumulative and differential DVHs are supported, as are sum plans. Eclipse uncertainty plans are supported. The measurement unit for absolute dose can be Gy, cGy, or eV/G for uncalibrated PRIMO files. The measurement unit for volume has to be cm$^{3}$. Multiple DVH text files can be read with \texttt{readDVH()} in one step.

Example: Read one Eclipse file \texttt{dvhFile.txt} from folder \texttt{"c:/folder"} and save the result in object \texttt{res}.\footnote{Note that the way to indicate the path to these files is different from the usual Windows style path: Instead of writing the backslash \texttt{"\textbackslash"} as folder separator, the forward slash \texttt{"/"} must be used.}

<<cReadData1a, eval=FALSE>>=
res <- readDVH("c:/folder/dvhFile.txt", type="Eclipse")
@

Basic information about the files can be displayed with \texttt{print()}, or just by entering the name of a DVH object at the prompt -- here used with built-in DVHs from three patients with radiotherapy, each with seven heart structures.\footnote{Sample data courtesy of Department of Radiation Oncology (Prof. Dr. Schmidberger), University Medical Center Mainz, Germany.}

<<cReadData1b, echo=TRUE>>=
print(dataMZ)
@

Display more information on structures with \texttt{verbose=TRUE}.
<<cReadData1c, echo=TRUE>>=
print(dataMZ, verbose=TRUE)
@

Multiple files with the same name pattern can be specified using wildcards like \texttt{*}. Example: Read all CadPlan files with the file name pattern \texttt{dvhFile*.txt} from folder \texttt{"c:/folder"} and save the result in object \texttt{res}.

<<cReadData2, eval=FALSE>>=
res <- readDVH("c:/folder/dvhFile*.txt", type="Cadplan")
@

When no file pattern is specified, multiple files can be selected using the standard Windows file picker dialogue. On MacOS and Linux, only a single file can be selected interactively.

<<cReadData4, eval=FALSE>>=
res <- readDVH(type="Eclipse")      # opens interactive file picker
@

For DVH files from a sum plan, prescribed dose can be encoded in the plan name like \texttt{name\_70Gy\_etc}. It will then be assumed that ``\% for dose'' is 100.

<<cReadData5, eval=FALSE>>=
res <- readDVH("c:/folder/*", type="Eclipse", planInfo="doseRx")
@

If files contain special characters, it may be necessary to specify the file encoding using options \texttt{encoding="UTF-8"} or \texttt{encoding="UTF-8-BOM"} (when a byte-order-mark is used). If Eclipse uncertainty plans are present, option \texttt{uncertainty=TRUE} is necessary. Structures with uncertainty plans get separate structure names by appending the actual structure name with a suffix derived from the uncertainty plan name.

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{DVH metrics}
\label{sec:metrics}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Calculate DVH metrics}
\label{sec:getMetric}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Function \texttt{getMetric()} calculates freely-defined DVH metrics based on data that has been read in as demonstrated in section \ref{sec:readData}. \texttt{getMetric()} has the following options:

\begin{itemize}
\item Option \texttt{x}: The DVH data.
\item Option \texttt{metric} -- one or many of the following:
\begin{itemize}
\item A \emph{pre-specified} DVH metric is one of the following character strings:
\begin{itemize}
\item \texttt{"DMEAN"}: The volume-weighted mean dose of the structure.
\item \texttt{"DMIN"}: The minimum dose of the non-zero-dose voxels in the structure.
\item \texttt{"DMAX"}: The maximum dose of the non-zero-dose voxels in the structure.
\item \texttt{"DSD"}: The standard deviation of the dose in the structure.
\item \texttt{"DRX"}: The prescription dose.
\item \texttt{"DHI"}: The Homogeneity Index according to ICRU 83: (D2\%-D98\%)/D50\%.
\item \texttt{"DEUD"}: The generalized equivalent uniform dose (gEUD, \citeNP{Niemierko1999,Wu2002}). This can be based on EQD$_{2}$ values if information on the fractionation is provided as well \cite{IAEA2008}.
\item \texttt{"DNTCP"}: The normal tissue complication probability (NTCP) according to the Lyman \citeyear{Lyman1985} probit model, the Niemierko \citeyear{Niemierko1999} logit model, the Poisson model (equations (1) to (3) in \citeNP{Kallman1992} with gEUD plugged in for $D$), or the relative seriality model (equation (18)). This can be based on EQD$_{2}$ values if information on the fractionation is provided as well.
\item \texttt{"DTCP"}: The tumor control probability (TCP) according to the same models as NTCP.
\end{itemize}
\item A \emph{free} DVH metric is a character string which has three mandatory elements and one optional element in the following order:
\begin{itemize}
\item 1$^{\mathrm{st}}$ letter \texttt{"D"} or \texttt{"V"}: \texttt{"D"} If the requested value is a dose, \texttt{"V"} if it is a volume.
\item 2$^{\mathrm{nd}}$ element \texttt{$\langle$number$\rangle$}: If the first letter is \texttt{"D"}, this gives the volume for which the dose value of the cumulative DVH should be reported. If the first letter is \texttt{"V"}, this gives the dose for which the volume value of the cumulative DVH should be reported.
\item 3$^{\mathrm{rd}}$ element \texttt{$\langle$measurement unit$\rangle$}: The measurement unit for the 2$^{\mathrm{nd}}$ element of the metric. Absolute volumes are indicated by \texttt{"CC"} for cm$^{3}$, relative volumes by \texttt{"\%"}. Absolute doses are indicated by \texttt{"Gy"} for Gray or \texttt{"cGy"} for Centigray, relative doses by \texttt{"\%"}.
\item Optional 4$^{\mathrm{th}}$ element \texttt{\_$\langle$measurement unit$\rangle$}: The measurement unit of the output value. Possible units are the same as for the 3$^{\mathrm{rd}}$ element. If missing, dose is reported as absolute dose in the measurement unit used in the DVH. If the measurement unit is missing, volume is reported as relative volume in \%.
\end{itemize}
\item Example metrics are listed in table \ref{tab:metrics}. By default, metrics are calculated using linear interpolation between adjacent supporting points of the cumulative DVH -- without extrapolating beyond the observed volume or dose. The interpolation method can be changed to use monotone Hermite splines \cite{Fritsch1980}, or to local polynomial regression with a Gaussian kernel \cite{Wand1995}. The kernel bandwidth is then determined by the direct plug-in method \cite{Ruppert1995}.
\end{itemize}

\item Option \texttt{patID}: Which patient IDs should be analyzed. With \texttt{fixed=FALSE}, IDs are interpreted as regular expressions matched against those found in the DVH files. By default, IDs are matched exactly. If missing, the metrics are calculated for all patients.

\item Option \texttt{structure}: Which structure should be analyzed. With \texttt{fixed=FALSE}, structure names are interpreted as regular expressions matched against those found in the DVH files. By default, structure names are matched exactly. If missing, the metrics are calculated for all structures.

\item Option \texttt{sortBy}: Results can be sorted according to these variables:
\begin{itemize}
\item \texttt{"observed"}: observed value of the metric
\item \texttt{"structure"}: structure for which the metric is calculated
\item \texttt{"metric"}: type of calculated metric
\item \texttt{"patID"}: patient ID
\end{itemize}

\item Option \texttt{splitBy}: Results can be divided into different tables according to these variables:
\begin{itemize}
\item \texttt{"structure"}: structure for which the metric is calculated
\item \texttt{"metric"}: type of calculated metric
\item \texttt{"patID"}: patient ID
\end{itemize}
\end{itemize}

\begin{table}[!htbp]
\caption{Examples of possible DVH metrics}
\label{tab:metrics}
\centering
\begin{tabular}{lllll}
Metric                & Reference       & Unit reference & Result    & Unit result \\\hline
\texttt{"V10Gy"}      & absolute dose   & Gy       & relative volume & \% \\
\texttt{"V10cGy\_CC"} & absolute dose   & cGy      & absolute volume & cm$^{3}$ \\
\texttt{"V10\%"}      & relative dose   & \%       & relative volume & \% \\
\texttt{"V10\%\_CC"}  & relative dose   & \%       & absolute volume & cm$^{3}$ \\\hline
\texttt{"D10CC"}      & absolute volume & cm$^{3}$ & absolute dose   & as in DVH \\
\texttt{"D10\%\_cGy"} & relative volume & \%       & absolute dose   & cGy \\\hline
\texttt{"DMEAN"}      & ---             & ---      & absolute dose   & as in DVH \\
\texttt{"DEUD"}       & ---             & ---      & absolute dose   & as in DVH \\
\texttt{"DSD"}        & ---             & ---      & absolute dose   & as in DVH \\
\texttt{"DMIN"}       & ---             & ---      & absolute dose   & as in DVH \\
\texttt{"DMAX"}       & ---             & ---      & absolute dose   & as in DVH \\\hline
\texttt{"DHI"}        & ---             & ---      & dose ratio      & --- \\
\texttt{"DNTCP"}      & ---             & ---      & probability     & --- \\
\texttt{"DTCP"}       & ---             & ---      & probability     & --- \\\hline
\end{tabular}
\end{table}

If volume or dose values outside the range of possible values for a structure are requested, it may be that metrics cannot be calculated, and the result will be \texttt{NA} (missing value) with a warning.

In the following examples, we use object \texttt{dataMZ} that is built into the \texttt{DVHmetrics} package. \texttt{dataMZ} was the result from reading three Eclipse DVH files, each with seven structures -- as demonstrated in section \ref{sec:readData}.

Calculate metric \texttt{DMEAN} for all structures for all patients in \texttt{dataMZ}.

<<cMetrics1>>=
getMetric(dataMZ, metric="DMEAN")
@

Calculate metric \texttt{D5cc} just for structure \texttt{HEART} for all patients in \texttt{dataMZ}.

<<cMetrics2>>=
getMetric(dataMZ, metric="D5cc", structure="HEART")
@

Calculate metric \texttt{D5cc} just for structure \texttt{HEART} for all patients in \texttt{dataMZ}, and sort result by the observed value of the metric.

<<cMetrics3>>=
getMetric(dataMZ, metric="D5cc", structure="HEART", sortBy="observed")
@

Calculate metrics \texttt{D10\%} and \texttt{V5Gy} for all structures containing the text \texttt{AMYOC} or \texttt{VALVE}, for patient IDs in \texttt{dataMZ} containing the text \texttt{23}, and sort result by metric and observed value.

<<cMetrics4>>=
getMetric(dataMZ, metric=c("D10%", "V5Gy"),
          structure=c("AMYOC", "VALVE"),
          patID="23",
          sortBy=c("metric", "observed"),
          fixed=FALSE)
@

Calculate metrics \texttt{DMEAN} and \texttt{D5cc} for structure \texttt{HEART} for all patients in \texttt{dataMZ}, sort by the observed value of the metric, and split the output such that one table is generated for each metric.

<<cMetrics5>>=
getMetric(dataMZ, metric=c("DMEAN", "D5cc"), structure="HEART",
          sortBy="observed", splitBy="metric")
@

Calculate metrics \texttt{DMEAN} and \texttt{D5cc} for structures \texttt{HEART} and \texttt{AOVALVE} for all patients in \texttt{dataMZ}, sort by observed value, and split the output such that one table is generated for each combination of structure and metric. Also store the result in object \texttt{met} that can be saved later.

<<cMetrics6>>=
met <- getMetric(dataMZ, metric=c("DMEAN", "D5cc"),
                 structure=c("HEART", "AOVALVE"),
                 sortBy="observed",
                 splitBy=c("structure", "metric"))
met                               # print the calculated results
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Save DVH metrics to file}
\label{sec:saveMetric}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

The calculated DVH metrics can be saved to tab-delimited text files with \texttt{saveMetric()}. These files are easy to import, e.\,g., into spreadsheets like Excel or into other statistics programs.

Assume object \texttt{met} has been calculated before as demonstrated in section \ref{sec:getMetric}. If \texttt{met} is not split into different tables, the following command saves \texttt{met} to the file \texttt{metrics.txt}. If \texttt{met} is divided into multiple tables, this saves \texttt{met} into different files that all have the name pattern \texttt{metrics\_NAME.txt}, where \texttt{NAME} stands, e.\,g., for the names of different structures.

<<cMetricsSave1, eval=FALSE>>=
saveMetric(met, file="c:/folder/metrics.txt")
@

Per default, numbers use the \texttt{.} as decimal separator. This can be changed with option \texttt{dec=","} .

<<cMetricsSave2, eval=FALSE>>=
saveMetric(met, file="c:/folder/metrics.txt", dec=",")
@

If text should be set in quotes in the output file, use \texttt{quote=TRUE} .

<<cMetricsSave3, eval=FALSE>>=
saveMetric(met, file="c:/folder/metrics.txt", quote=TRUE)
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Convenience functions for DMEAN, gEUD, NTCP, TCP}
\label{sec:DMEAN}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

DMEAN, gEUD, NTCP and TCP may be calculated together with other metrics using \texttt{getMetric()}, but there are specialized convenience functions for this task as well. In particular, \texttt{getDMEAN()} calculates the dose mean, median, mode, minimum, and maximum based on the (interpolated) differential DVH instead of relying on the values exported by the TPS.

<<cDmean1>>=
dmean <- getDMEAN(dataMZ[[1]])
subset(dmean, select=c(doseAvg, doseMed, doseMin, doseMax))
@

gEUD is calculated by \texttt{getEUD()}.

<<cDmean2>>=
# note that different tissues should have different parameter values,
# this is just for demonstration purposes
getEUD(dataMZ[[1]], EUDa=2)
@

NTCP and TCP are calculated by \texttt{getNTCP()} and \texttt{getTCP()}, respectively.

<<cDmean3>>=
# note that different tissues should have different parameter values,
# this is just for demonstration purposes
getNTCP(dataMZ[[1]], NTCPtd50=40, NTCPm=0.6, NTCPn=0.5, NTCPtype="probit")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Point-wise mean DVH with standard deviations}
\label{sec:pointWise}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Function \texttt{getMeanDVH()} returns the point-wise mean and median DVH with the point-wise standard deviation for a given list of input DVHs. Other point-wise measures may be calculated as well. Before calculating the point-wise mean and SD, DVHs are first linearly interpolated such that they possess the same set of nodes. This feature can be useful for evaluating different plan options: The DVHs for each plan need to be exported using a different patient ID which thus serves as a plan identifier. See section \ref{sec:showDVH} to show the mean DVH with SD regions.

<<cPointWise1>>=
# point-wise mean and SD for structure HEART over all patients
m1 <- getMeanDVH(dataMZ, fun=list(M=mean, SD=sd), byPat=FALSE, structure="HEART")
head(m1)
@

When using option \texttt{returnDVHObj=TRUE}, the function returns a DVH object that behaves like a regular DVH, and can be used in functions such as \texttt{showDVH()} or \texttt{getMetric()}

<<cPointWise2>>=
# point-wise mean for structure HEART over all patients
m2 <- getMeanDVH(dataMZ, fun=list(mean), byPat=FALSE, structure="HEART",
                 returnDVHObj=TRUE)

getMetric(m2, metric="V5GY")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{DVH diagrams}
\label{sec:plot}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Plot DVH diagrams}
\label{sec:showDVH}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Cumulative as well as differential DVH diagrams can be generated with \texttt{showDVH()}. If you are using \textsf{R}Studio or Architect, all produced diagrams are accessible in the plots tab by clicking on the left and right arrows. Depending on the option \texttt{byPat}, each DVH diagram either shows one patient with multiple structures (\texttt{byPat=TRUE}) or one structure with multiple patients (\texttt{byPat=FALSE}).

<<cPlots1, out.width='3in'>>=
showDVH(dataMZ, byPat=TRUE)
@

Patient IDs and structures can be selected with the \texttt{patID="$\langle$ID$\rangle$"} option and the \texttt{structure="$\langle$NAME$\rangle$"} option. With \texttt{fixed=FALSE}, both accept regular expressions. By default, IDs and structure names are matched exactly. By default, all patients/structures are shown.

<<cPlots2, out.width='3in'>>=
showDVH(dataMZ, byPat=FALSE, patID=c("P123", "P234"))
@

By default, the relative DVH is shown. Absolute volume (if available in the input files) can be plotted with the \texttt{rel=FALSE} option. For differential DVH, set \texttt{cumul=FALSE}.

<<cPlots3, out.width='3in'>>=
# match structures containing "VALVE" and "AMYOC"
showDVH(dataMZ, cumul=FALSE, rel=FALSE,
        structure=c("VALVE", "AMYOC"), fixed=FALSE)
@

Option \texttt{thresh} allows to restrict the range of the $x$-axis such that only relative volumes larger than \texttt{thresh} appear. Use option \texttt{show=FALSE} to prevent the diagrams from being shown if you just need the returned object (here: \texttt{dvhPlot}) to later save the diagrams to file.

<<cPlots4, out.width='3in'>>=
# just save the diagram but don't show it
dvhPlot <- showDVH(dataMZ, structure=c("HEART", "AOVALVE", "AVNODE"),
                   rel=FALSE, thresh=0.001, show=FALSE)
@

With option \texttt{addMSD=TRUE}, the diagram also shows the point-wise mean DVH as a black curve as well as grey shaded areas for the regions defined by the point-wise 1 standard deviation and 2 standard deviations around this mean (see section \ref{sec:pointWise}). With \texttt{byPat=FALSE} and \texttt{addMSD=TRUE}, the average DVH for one structure over many patients can be visually assessed.

<<cPlots5, out.width='3in'>>=
# add point-wise mean DVH and 1 SD/2 SD regions
showDVH(dataMZ, structure="HEART", byPat=FALSE, addMSD=TRUE)
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Save cumulative DVH diagrams to file}
\label{sec:saveDVH}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

DVH diagrams can be saved to file using \texttt{saveDVH()}. A file name pattern can then be supplied to option \texttt{file}. By using different file extensions like \texttt{.pdf}, \texttt{.jpg}, \texttt{.png}, different graphics formats can be automatically selected. In addition, the width and height of the diagram can be specified in inch.

<<cPlotsSave, eval=FALSE>>=
saveDVH(dvhPlot, file="c:/folder/dvh.pdf", width=7, height=5)
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{Quality assurance constraints on the dose-volume relationship}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

For quality assurance, it is possible to define, check, and visualize constraints on the dose-volume relationship for DVHs.

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Define constraints}
\label{sec:defConstraints}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

A DVH constraint is a character string that consists of three parts: The DVH metric (see section \ref{sec:getMetric}), a comparison operator among \texttt{<, >, <=, >=}, and the reference value together with the measurement unit -- one among among \texttt{Gy, cGy, cc, \%}. For constraints involving the relative dose, the DVH must contain the prescription dose.

Some example constraints are \texttt{"V10Gy > 80\%"} (more than 80\% of the structure should have received 10Gy), \texttt{"V20\% < 10CC"} (less than 10cm$^{3}$ of the structure should have received 20\% of the prescription dose), or \texttt{"D10CC > 500cGy"} (The ``hottest'' 10cm$^{3}$ of the structure should have received more than 500cGy). Constraints can also apply to the dose mean, median, and standard deviation as well as to the gEUD and to the (N)TCP.

A DVH constraint can apply to a specific patient or to all patients, and to a specific structure or to all structures.
\begin{itemize}
\item If constraints apply to all patients/structures, the constraint can be a \texttt{character} vector with elements like the examples above.
\item If constraints apply only to some patients/structures, the constraint must be a data frame with variables \texttt{constraint}, \texttt{patID} and \texttt{structure}. Each row then defines one constraint and its scope: \texttt{constraint} must be a character string with one constraint definition as in the examples above. \texttt{patID} must be either a character string with a valid patient ID, or \texttt{"*"} if the the constraint applies to all patients. \texttt{structure} must be either a character string with a valid structure name, or \texttt{"*"} if the the constraint applies to all structures. If variable \texttt{patID} is missing from the data frame, the constraints apply to all available patients. If variable \texttt{structure} is missing from the data frame, the constraints apply to all available structures.
\end{itemize}

Alternatively, it is possible to specify a set of constraints as a table in a text file with one row per constraint and one column for the constraint expression, structure, and patient ID. A table like this can be created in a spreadsheet program like Excel (fig.\ \ref{fig:constrExcel}), be exported to a tab-delimited text-file, and be read in by function \texttt{readConstraints()}. Table \ref{tab:constrPaste} shows some examples.

\begin{table}
\caption{Example for pasted constraints.}
\label{tab:constrPaste}
\begin{tabular}{lll}
\multicolumn{3}{l}{Constraints that apply to all patients and to all structures}\\\hline
\texttt{"D10cc < 20\%"}  & &\\
\texttt{"V5cGy > 100cc"} & &\\
\texttt{"DMEAN < 10Gy"}  & &\\
& & \\
\multicolumn{3}{l}{Constraints that apply to some patients and to all structures}\\\hline
\texttt{"constraint"}    & \texttt{"patID"} &\\
\texttt{"D10cc < 20\%"}  & \texttt{"P123"}  &\\
\texttt{"V5cGy > 100cc"} & \texttt{"*"}     &\\
\texttt{"DMEAN < 10Gy"}  & \texttt{"P234"}  &\\
& & \\
\multicolumn{3}{l}{Constraints that apply to some patients and to some structures}\\\hline
\texttt{"constraint"}    & \texttt{"patID"} & \texttt{"structure"}\\
\texttt{"D10cc < 20\%"}  & \texttt{"P123"}  & \texttt{"*"}\\
\texttt{"V5cGy > 100cc"} & \texttt{"*"}     & \texttt{"HEART"}\\
\texttt{"DMEAN < 10Gy"}  & \texttt{"P234"}  & \texttt{"AOVALVE"}\\\hline
\end{tabular}
\end{table}


\begin{figure}[ht]
\centering
\includegraphics[width=10cm]{constraintsExcel}
\caption{Defining constraints in a spreadsheet program like Excel}
\label{fig:constrExcel}
\end{figure}

<<cConstrDef2, eval=FALSE>>=
dataConstr <- readConstraints("constraints.txt", dec=".", sep="\t")
@

The constraint data frame \texttt{dataConstr} is built into \texttt{DVHmetrics} and applies to the \texttt{dataMZ} DVH data.

<<cConstrDef3, echo=TRUE>>=
dataConstr     # show defined constraints and their scope
@

For checking constraints, and for calculating the difference between the observed DVH and the constraint, the DVH is linearly interpolated, by using monotone Hermite splines, or by local polynomial kernel regression.

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Check constraints}
\label{sec:checkConstraint}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Constraints are checked with \texttt{checkConstraint()}. The output returns information on the observed value of the tested metric, on the compliance with respect to this metric, and on the absolute/relative deviation in volume as well as in dose to the specified constraint value. The units for the absolute deviation are those used in the constraint expression. When the constraint defines a point in dose-volume space, \texttt{checkConstraint()} reports another quantitative measure for the degree of violation: The closest point on the DVH to the constraint as well as its Euclidean distance to the constraint point.

For calculating the minimal Euclidean distance between the constraint point and the DVH, the constraint point is orthogonally projected onto each DVH segment between (interpolated) DVH nodes. The relative Euclidean distance is the minimum of these distances divided by the distance of the constraint point to the closest axis (dose and volume) along the same direction. In doing so, the devation from the expected volume per dose and the devation from the expected dose per volume are condensed in a single metric.

As an example, we use the DVHs and corresponding constraints that are built into the \texttt{DVHmetrics} package.

<<cConstrCheck1, echo=TRUE>>=
## store result in object cc to save to file later
cc <- checkConstraint(dataMZ, constr=dataConstr)
print(cc, digits=2)            # show output with 2 decimal places
@

The result from a constraint check can be saved with function \texttt{saveConstraint()} that works like \texttt{saveMetric()} (see section \ref{sec:saveMetric}).

<<cConstrCheck2, eval=FALSE>>=
saveConstraint(cc, file="c:/folder/constrCheck.txt")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\subsection{Visualize constraints}
\label{sec:showConstraint}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

Constraints that define a point in dose-volume space can be visualized in a DVH with relative volume or absolute volume. The constraints will be converted to match the DVH plot. Only patients and structures within the scope of the defined constraints are shown. The diagram also shows the point on the DVH closest to the constraint. This can be verified visually only if the aspect ratio of the diagram is 1.

As in \texttt{showDVH()} (see section \ref{sec:showDVH}), either one diagram per patient with multiple structures is shown (\texttt{byPat=TRUE}), or one diagram per structure with multiple patients (\texttt{byPat=FALSE}).

<<cConstrShow1, out.width='3in', echo=TRUE>>=
## plot relative volume
showConstraint(dataMZ, constr=dataConstr, byPat=TRUE)
@

<<cConstrShow2, eval=FALSE>>=
## plot absolute volume - store result in sc to save to file later
sc <- showConstraint(dataMZ, constr=dataConstr,
                     byPat=FALSE, rel=FALSE)
@

The result can be saved using \texttt{saveDVH()} as demonstrated in section \ref{sec:saveDVH}.

<<cConstrShow3, eval=FALSE>>=
saveDVH(sc, file="c:/folder/dvhConstraint.pdf")
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section{BED, EQD2, Isoeffective Dose}
\label{sec:BED}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

The linear-quadratic model for the proportion of surviving cells $S$ after dose $d$ is \cite{IAEA2008}:

$$
S = e^{-(\alpha d + \beta d^{2})}
$$

According to the model, the biologically effective dose (BED) with total dose $D$, fraction dose $d$ and a tissue-dependent $\alpha / \beta$ ratio is:

$$
\mathrm{BED} = D \left[1 + \frac{d}{\alpha / \beta}\right]
$$

Given two different fractionation schemes, the total dose $D_{2}$ for the new fraction dose $d_{2}$ that corresponds to the total dose $D_{1}$ from the reference fraction dose $d_{1}$ can be calculated from solving the following equation for the desired measure:

$$
\frac{D_{2}}{D_{1}} = \frac{d_{1} + (\alpha / \beta)}{d_{2} + (\alpha / \beta)}
$$

As a special case, the dose in 2Gy fractions biologically equivalent dose (EQD$_{2}$) is given by:

$$
\mathrm{EQD}_{2} = D_{1} \cdot \frac{d_{1} + (\alpha / \beta)}{2 + (\alpha / \beta)} = \frac{\mathrm{BED}}{1 + \frac{2}{\alpha / \beta}}
$$

The following convenience functions allow for easy calculation of these measures:

<<cBED1>>=
getBED(D=50, fd=2.5, ab=c(2, 3, 4))
getEQD2(D=50, fd=2.5, ab=c(2, 3, 4))
getIsoEffD(D1=70, fd1=2, fd2=3, ab=c(3.5, 10))
@

The same functions can be used to convert complete DVHs to BED, EQD2, or to the iso-effective dose corresponding to some other fraction dose.

<<cBED2>>=
getEQD2(D=dataMZ[[c(1, 1)]], fd=2.5, ab=3)
@

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------
\section*{Acknowledgements}
%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

The authors thank Marcus Stockinger for ideas on checking quality assurance constraints as well as Sandra B\"{u}hrdel, Hannes Rennau, Ulrich Wolf, Bjorne Riis, Nico Banz, and Michael R.\ Young for example DVH files exported from different treatment planning systems.

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

\bibliographystyle{apacite}
\renewcommand{\BAvailFrom}{URL\ }
\renewcommand{\APACrefURL}{URL\ }
\bibliography{lit}

%%%%---------------------------------------------------------------------------
%%%%---------------------------------------------------------------------------

\end{document}
