## convert DVH to cumulative/differential DVH
## with smoothing the differential DVH
convertDVHsmooth <-
function(x, toType=c("asis", "cumulative", "differential"),
toDoseUnit=c("asis", "GY", "CGY"),
interp=c("asis", "linear", "spline", "ksmooth", "smoothSpl"),
nodes=NULL, rangeD=NULL, perDose=TRUE) {
UseMethod("convertDVHsmooth")
}
## for DVH matrix itself
convertDVHsmooth.matrix <-
function(x, toType=c("asis", "cumulative", "differential"),
toDoseUnit=c("asis", "GY", "CGY"),
interp=c("asis", "linear", "spline", "ksmooth", "smoothSpl"),
nodes=NULL, rangeD=NULL, perDose=TRUE) {
toType <- match.arg(toType)
toDoseUnit <- match.arg(toDoseUnit)
interp <- match.arg(interp)
if(!is.null(nodes)) { stopifnot(nodes > 2L) }
## matrix may include duplicate rows
x <- unique(x)
## when dose entries are duplicated with different volumes (cumulative DVH),
## pick row with max volume
## split matrix into rows with unique dose
## preserve matrix class and keep col names with split.data.frame()
xL <- if(!all(is.na(x[ , "dose"]))) {
split.data.frame(x, x[ , "dose"])
} else {
split.data.frame(x, x[ , "doseRel"])
}
## function to keep row with max volume
keepMaxVol <- function(xSub) {
idx <- if(!all(is.na(xSub[ , "volume"]))) {
xSub[ , "volume"] == max(xSub[ , "volume"])
} else{
xSub[ , "volumeRel"] == max(xSub[ , "volumeRel"])
}
xSub[idx, ]
}
## apply keepMaxVol and re-bind to matrix
x <- do.call("rbind", lapply(xL, keepMaxVol))
rownames(x) <- NULL
## determine number of nodes for interpolation
nodes <- if(!is.null(nodes)) {
if(nodes < 2L) {
warning("nodes is < 2 which is too few, - set to 2")
2L
} else {
nodes
}
} else {
nrow(x)
}
## absolute and relative doses/volumes
dose <- x[ , "dose"]
doseRel <- x[ , "doseRel"]
volume <- x[ , "volume"]
volumeRel <- x[ , "volumeRel"]
## convert dose unit
doseConv <- if(toDoseUnit == "asis") {
dose ## nothing to do
} else if(toDoseUnit == "GY") {
dose / 100 ## cGy to Gy
} else if(toDoseUnit == "CGY") {
dose * 100 ## Gy to cGy
}
## convert DVH type
if(toType == "asis") {
## nothing to convert
## interpolate?
if(interp != "asis") {
## if differential && perDose == FALSE
## normalize -> interpolate -> re-normalize
## check if volume is already sorted -> cumulative DVH
volumeSel <- if(!any(is.na(volumeRel))) {
volumeRel
} else {
volume
}
if(isTRUE(all.equal(volumeSel, sort(volumeSel, decreasing=TRUE)))) {
DVHtype <- "cumulative"
splMethod <- "monoH.FC" # for spline interpolation -> monotonic
} else {
DVHtype <- "differential"
splMethod <- "fmm" # for spline interpolation
}
if((DVHtype == "differential") && !perDose) {
## if perDose == FALSE: normalize -> interpolate -> re-normalize
binW <- diff(c(-doseConv[1], doseConv))
volumeNew <- -diff(volume) / binW
volumeRelNew <- -diff(volumeRel) / binW
}
smFun <- switch(interp,
linear=getInterpLin,
ksmooth=getKSmooth,
smoothSpl=getSmoothSpl,
getInterpSpl) # default
sm <- smFun(doseConv, doseRel, volume, volumeRel, nodes=nodes,
rangeD=rangeD, method=splMethod)
doseNew <- sm[ , "dose"]
doseRelNew <- sm[ , "doseRel"]
if((DVHtype == "differential") && !perDose) {
## re-normalize to non-per-dose volume
## start bin at 0 because no out-of-range interpolation
binW <- diff(c(0, doseConv))
volumeNew <- sm[ , "volume"] * binW
volumeRelNew <- sm[ , "volumeRel"] * binW
} else {
volumeNew <- sm[ , "volume"]
volumeRelNew <- sm[ , "volumeRel"]
}
} else {
doseNew <- doseConv
doseRelNew <- doseRel
volumeNew <- volume
volumeRelNew <- volumeRel
}
} else if(toType == "cumulative") {
## convert from differential to cumulative DVH
## interpolate?
if(interp != "asis") {
if(!perDose) {
## if perDose == FALSE: normalize -> interpolate -> re-normalize
binW <- diff(c(-doseConv[1], doseConv))
volumeNew <- -diff(volume) / binW
volumeRelNew <- -diff(volumeRel) / binW
}
smFun <- switch(interp,
linear=getInterpLin,
ksmooth=getKSmooth,
smoothSpl=getSmoothSpl,
getInterpSpl) # default
sm <- smFun(doseConv, doseRel, volume, volumeRel, nodes=nodes,
rangeD=rangeD, method="monoH.FC")
doseConv <- sm[ , "dose"]
doseRel <- sm[ , "doseRel"]
if(!perDose) {
## re-normalize to non-per-dose volume
## start bin at 0 because no out-of-range interpolation
binW <- diff(c(0, doseConv))
volume <- sm[ , "volume"] * binW
volumeRel <- sm[ , "volumeRel"] * binW
} else {
volume <- sm[ , "volume"]
volumeRel <- sm[ , "volumeRel"]
}
}
## dose category half-widths
doseCatHW <- diff(doseConv)/2
doseRelCatHW <- diff(doseRel)/2
## dose category mid-points, starting at 0
N <- nrow(x)
doseMidPt <- c(0, doseConv[-N] + doseCatHW)
doseRelMidPt <- c(0, doseRel[-N] + doseRelCatHW)
## add one more category beyond max dose
doseNew <- c(doseMidPt, doseConv[N] + doseCatHW[N-1])
doseRelNew <- c(doseRelMidPt, doseRel[N] + doseRelCatHW[N-1])
if(perDose) {
## differential DVH -> volume is per Gy -> mult with bin-width
binW <- diff(c(-doseConv[1], doseConv))
volumeBin <- volume*binW
volumeRelBin <- volumeRel*binW
} else {
volumeBin <- volume
volumeRelBin <- volumeRel
}
## volume "survival" curve starting at max volume
volumeNew <- sum(volumeBin) - c(0, cumsum(volumeBin))
volumeRelNew <- sum(volumeRelBin) - c(0, cumsum(volumeRelBin))
} else if(toType == "differential") {
## convert from cumulative to differential DVH
## dose category half-widths
doseCatHW <- diff(doseConv)/2
doseRelCatHW <- diff(doseRel)/2
doseNew <- doseConv[-length(doseConv)] + doseCatHW
doseRelNew <- doseRel[-length(doseRel)] + doseRelCatHW
## differential DVH -> volume is per Gy -> divide by bin-width
binW <- 2*doseCatHW
volumeNew <- if(perDose) {
-diff(volume) / binW
} else {
-diff(volume)
}
volumeRelNew <- if(perDose) {
-diff(volumeRel) / binW
} else {
-diff(volumeRel)
}
## interpolate?
if(interp != "asis") {
## if perDose == FALSE: normalize -> interpolate -> re-normalize
if(!perDose) {
volumeNew <- -diff(volume) / binW
volumeRelNew <- -diff(volumeRel) / binW
}
## interpolate
smFun <- switch(interp,
linear=getInterpLin,
ksmooth=getKSmooth,
smoothSpl=getSmoothSpl,
getInterpSpl) # default
sm <- smFun(doseNew, doseRelNew, volumeNew, volumeRelNew,
nodes=nodes, rangeD=rangeD, method="fmm")
doseNew <- sm[ , "dose"]
doseRelNew <- sm[ , "doseRel"]
if(!perDose) {
## re-normalize to non-per-dose volume
## start bin at 0 because no out-of-range interpolation
binW <- diff(c(0, doseNew))
volumeNew <- sm[ , "volume"] *binW
volumeRelNew <- sm[ , "volumeRel"]*binW
} else {
volumeNew <- sm[ , "volume"]
volumeRelNew <- sm[ , "volumeRel"]
}
}
}
cbind(dose=doseNew, doseRel=doseRelNew,
volume=volumeNew, volumeRel=volumeRelNew)
}
## for 1 DVH object
convertDVHsmooth.DVHs <-
function(x, toType=c("asis", "cumulative", "differential"),
toDoseUnit=c("asis", "GY", "CGY"),
interp=c("asis", "linear", "spline", "ksmooth", "smoothSpl"),
nodes=NULL, rangeD=NULL, perDose=TRUE) {
toType <- match.arg(toType)
toDoseUnit <- match.arg(toDoseUnit)
if(!is.null(nodes)) { stopifnot(nodes > 2) }
## copy old DVH structure and convert DVH as well as other dose info
DVH <- x
if((toType == "asis") && (toDoseUnit == "asis") && (interp != "asis")) {
## just interpolate existing DVHs
## cumulative
DVH$dvh <- convertDVHsmooth(x$dvh, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
## differential
if(!is.null(x$dvhDiff)) {
DVH$dvhDiff <- convertDVHsmooth(x$dvhDiff, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
}
} else if((toType == "asis") && (toDoseUnit != "asis") && (interp == "asis")) {
## just change dose unit in DVH and remaining dose values
cf <- if( (toupper(x$doseUnit) == "CGY") && (toDoseUnit == "GY")) {
1/100
} else if((toupper(x$doseUnit) == "GY") && (toDoseUnit == "CGY")) {
100
} else {
NA_real_
}
## cumulative
DVH$dvh <- convertDVHsmooth(x$dvh, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
## differential
if(!is.null(x$dvhDiff)) {
DVH$dvhDiff <- convertDVHsmooth(x$dvhDiff, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
}
DVH$dvh[ , "dose"] <- cf*x$dvh[ , "dose"]
DVH$doseMin <- cf*x$doseMin
DVH$doseMax <- cf*x$doseMax
DVH$doseRx <- cf*x$doseRx
DVH$isoDoseRx <- cf*x$isoDoseRx
DVH$doseAvg <- cf*x$doseAvg
DVH$doseMed <- cf*x$doseMed
DVH$doseMode <- cf*x$doseMode
DVH$doseSD <- cf*x$doseSD
DVH$doseUnit <- toDoseUnit
} else if((toType == "asis") && (toDoseUnit != "asis") && (interp != "asis")) {
## change dose unit in DVH and remaining dose values
## and interpolate
cf <- if( (toupper(x$doseUnit) == "CGY") && (toDoseUnit == "GY")) {
1/100
} else if((toupper(x$doseUnit) == "GY") && (toDoseUnit == "CGY")) {
100
} else {
NA_real_
}
DVH$dvh[ , "dose"] <- cf*x$dvh[ , "dose"]
DVH$doseMin <- cf*x$doseMin
DVH$doseMax <- cf*x$doseMax
DVH$doseRx <- cf*x$doseRx
DVH$isoDoseRx <- cf*x$isoDoseRx
DVH$doseAvg <- cf*x$doseAvg
DVH$doseMed <- cf*x$doseMed
DVH$doseMode <- cf*x$doseMode
DVH$doseSD <- cf*x$doseSD
DVH$doseUnit <- toDoseUnit
} else if((toType != "asis") && (toDoseUnit == "asis")) {
## just change DVH type
if(toType == "differential") {
## from cumulative to differential
DVH$dvhDiff <- convertDVHsmooth(x$dvh, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
} else {
## from differential to cumulative
if(!is.null(x$dvhDiff)) {
DVH$dvh <- convertDVHsmooth(x$dvhDiff, toType=toType,
toDoseUnit="asis", interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
} else {
warning("No differential DVH found. Left cumulative DVH as is.")
}
}
} else if((toType != "asis") && (toDoseUnit != "asis")) {
## change DVH type and dose unit
cf <- if( (toupper(x$doseUnit) == "CGY") && (toDoseUnit == "GY")) {
1/100
} else if((toupper(x$doseUnit) == "GY") && (toDoseUnit == "CGY")) {
100
} else {
NA_real_
}
if(toType == "differential") {
## from cumulative to differential
DVH$dvhDiff <- convertDVHsmooth(x$dvh, toType=toType,
toDoseUnit=toDoseUnit, interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
} else {
## from differential to cumulative
if(!is.null(x$dvhDiff)) {
DVH$dvh <- convertDVHsmooth(x$dvhDiff, toType=toType,
toDoseUnit=toDoseUnit, interp=interp,
nodes=nodes, rangeD=rangeD, perDose=perDose)
} else {
warning("No differential DVH found. Left cumulative DVH as is.")
}
}
DVH$doseMin <- cf*x$doseMin
DVH$doseMax <- cf*x$doseMax
DVH$doseRx <- cf*x$doseRx
DVH$isoDoseRx <- cf*x$isoDoseRx
DVH$doseAvg <- cf*x$doseAvg
DVH$doseMed <- cf*x$doseMed
DVH$doseMode <- cf*x$doseMode
DVH$doseSD <- cf*x$doseSD
DVH$doseUnit <- toDoseUnit
}
return(DVH)
}
## for list of DVH objects
convertDVHsmooth.DVHLst <-
function(x, toType=c("asis", "cumulative", "differential"),
toDoseUnit=c("asis", "GY", "CGY"),
interp=c("asis", "linear", "spline", "ksmooth", "smoothSpl"),
nodes=NULL, rangeD=NULL, perDose=TRUE) {
toType <- match.arg(toType)
toDoseUnit <- match.arg(toDoseUnit)
interp <- match.arg(interp)
if(!is.null(nodes)) { stopifnot(nodes > 2) }
dvhL <- Map(convertDVHsmooth, x, toType=toType, toDoseUnit=toDoseUnit,
interp=interp, nodes=list(nodes), rangeD=list(rangeD),
perDose=perDose)
names(dvhL) <- names(x)
class(dvhL) <- "DVHLst"
attr(dvhL, which="byPat") <- attributes(x)$byPat
return(dvhL)
}
## for list of list of DVH objects
convertDVHsmooth.DVHLstLst <-
function(x, toType=c("asis", "cumulative", "differential"),
toDoseUnit=c("asis", "GY", "CGY"),
interp=c("asis", "linear", "spline", "ksmooth", "smoothSpl"),
nodes=NULL, rangeD=NULL, perDose=TRUE) {
toType <- match.arg(toType)
toDoseUnit <- match.arg(toDoseUnit)
interp <- match.arg(interp)
if(!is.null(nodes)) { stopifnot(nodes > 2) }
dvhLL <- Map(convertDVHsmooth, x, toType=toType, toDoseUnit=toDoseUnit,
interp=interp, nodes=list(nodes), rangeD=list(rangeD),
perDose=perDose)
names(dvhLL) <- names(x)
class(dvhLL) <- "DVHLstLst"
attr(dvhLL, which="byPat") <- attributes(x)$byPat
return(dvhLL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.