## S3 generic function
showConstraint <-
function(x, constr, byPat=TRUE, rel=TRUE, guessX=TRUE, guessY=TRUE,
thresh=1, show=TRUE, visible=FALSE) {
UseMethod("showConstraint")
}
showConstraint.DVHs <-
function(x, constr, byPat=TRUE, rel=TRUE, guessX=TRUE, guessY=TRUE,
thresh=1, show=TRUE, visible=FALSE) {
x <- if(byPat) {
setNames(list(x), x$structure)
} else {
setNames(list(x), x$patID)
}
class(x) <- "DVHLst"
attr(x, which="byPat") <- byPat
#NextMethod("showConstraint")
showConstraint.DVHLst(x, constr=constr, byPat=byPat,
guessX=guessX, guessY=guessY, rel=rel,
thresh=thresh, show=show, visible=visible)
}
## for byPat=TRUE
## x is a list of DVHs (1 per structure)
## for byPat=FALSE
## x is a list of DVHs (1 per patient)
showConstraint.DVHLst <-
function(x, constr, byPat=TRUE, rel=TRUE, guessX=TRUE, guessY=TRUE,
thresh=1, show=TRUE, visible=FALSE) {
## make sure DVH list is organized as required for byPat
if(is.null(attributes(x)$byPat) || attributes(x)$byPat != byPat) {
stop(c("DVH list organization by-patient / by-structure",
" either could not be determined or is different from byPat"))
}
## restrict x and constr to the same IDs/structures
xConstrSub <- harmoConstrDVH(x, constr=constr, byPat=byPat)
## determine dose and volume units in x
xDoseUnits <- vapply(xConstrSub$x, function(y) y$doseUnit, character(1))
xVolUnits <- vapply(xConstrSub$x, function(y) y$volumeUnit, character(1))
## parse constraint and convert dose/volume units if necessary
constrParse <- Map(parseConstraint, xConstrSub$constr,
doseUnit=xDoseUnits,
volUnit=xVolUnits)
## data frame from parsed constraint list
cDFL <- lapply(constrParse, function(y) data.frame(y, stringsAsFactors=FALSE))
## cDF <- melt(cDFL, id.vars=c("constraint", "valid", "metric", "DV",
## "valRef", "unitRef", "cmp", "unitCmp", "metricInv"),
## value.name="valCmp")
## cDF$variable <- NULL
## names(cDF)[names(cDF) == "L1"] <- if(byPat) {
## "structure"
## } else {
## "patID"
## }
cDF <- do.call("rbind", cDFL)
if(byPat) {
cDF[["structure"]] <- names(cDFL)
} else {
cDF[["patID"]] <- names(cDFL)
}
rownames(cDF) <- NULL
## distance constraint to closest point on DVH curve
getMinDstPt <- function(dvh, cnstr) {
Dcoord <- ifelse(cnstr$DV == "D", cnstr$valCmp, cnstr$valRef)
Vcoord <- ifelse(cnstr$DV == "V", cnstr$valCmp, cnstr$valRef)
Dcoord <- ifelse(cnstr$valid, Dcoord, NA_real_)
Vcoord <- ifelse(cnstr$valid, Vcoord, NA_real_)
volRel <- ((cnstr$DV == "D") & (cnstr$unitRef == "%")) |
((cnstr$DV == "V") & (cnstr$unitCmp == "%"))
doseRel <- ((cnstr$DV == "D") & (cnstr$unitCmp == "%")) |
((cnstr$DV == "V") & (cnstr$unitRef == "%"))
dstL <- dvhDistance(dvh,
DV=data.frame(D=Dcoord, V=Vcoord,
volRel=volRel, doseRel=doseRel,
unitRef=cnstr$unitRef, unitCmp=cnstr$unitCmp))
setNames(dstL, cnstr$constraint)
}
minDst <- Map(getMinDstPt, xConstrSub$x, constrParse)
minDstL <- melt_dw(minDst)
## minDstDF <- dcast(minDstL, L1 + L2 ~ L3, value.var="value")
minDstDF <- reshape(minDstL, direction="wide", v.names="value",
timevar="L3", idvar=c("L1", "L2"))
minDstDF <- setNames(minDstDF, gsub("^value\\.(.+)$", "\\1", names(minDstDF)))
names(minDstDF)[names(minDstDF) == "L1"] <- if(byPat) {
"structure"
} else {
"patID"
}
names(minDstDF)[names(minDstDF) == "L2"] <- "constraint"
## merge distance to DVH with remaining data
allDF <- merge(cDF, minDstDF)
## from given info, calculate dose, volume, relative volume
getDoseVol <- function(DV, valRef, unitRef, cmp, unitCmp, valCmp, DVH,
ptMinD, ptMinV, y) {
if(!(DVH %in% names(y))) {
warning(c("DVH for ", DVH, " not found"))
return(NA)
}
if(is.na(ptMinD) || is.na(ptMinV) || is.na(valRef)) {
doseAbs <- NA_real_
volAbs <- NA_real_
volRel <- NA_real_
ptMinDabs <- NA_real_
ptMinVabs <- NA_real_
ptMinVrel <- NA_real_
} else if(DV == "V") {
## input is absolute or relative dose
## output is absolute or relative volume
if(unitRef == "%") {
## input is relative dose -> convert rel to abs
## -> we need doseRx
stopifnot( !is.null(y[[DVH]]$doseRx),
all(!is.na( y[[DVH]]$doseRx)))
doseAbs <- (valRef/100) * y[[DVH]]$doseRx
ptMinDabs <- (ptMinD/100) * y[[DVH]]$doseRx
} else if(unitRef %in% c("GY", "CGY")) {
## input is absolute dose -> nothing to do
doseAbs <- valRef
ptMinDabs <- ptMinD
}
if(unitCmp == "%") {
## output is relative volume -> convert rel to abs
volAbs <- (valCmp/100) * y[[DVH]]$structVol
volRel <- valCmp
ptMinVabs <- (ptMinV/100) * y[[DVH]]$structVol
ptMinVrel <- ptMinV
} else if(unitCmp == "CC") {
## output is absolute volume -> convert abs to rel
volAbs <- valCmp
volRel <- 100*valCmp / y[[DVH]]$structVol[1]
ptMinVabs <- ptMinV
ptMinVrel <- 100*ptMinV / y[[DVH]]$structVol[1]
}
} else if(DV == "D") {
## input is absolute or relative volume
## output is absolute or relative dose
if(unitRef == "%") {
## input is relative volume -> convert rel to abs
volAbs <- (valRef/100) * y[[DVH]]$structVol
volRel <- valRef
ptMinVabs <- (ptMinV/100) * y[[DVH]]$structVol
ptMinVrel <- ptMinV
} else if(unitRef == "CC") {
## input is absolute volume -> convert abs to rel
volAbs <- valRef
volRel <- 100*valRef / y[[DVH]]$structVol
ptMinVabs <- ptMinV
ptMinVrel <- 100*ptMinV / y[[DVH]]$structVol
}
if(unitCmp == "%") {
## output is relative dose -> convert rel to abs
## -> we need doseRx
stopifnot( !is.null(y[[DVH]]$doseRx),
all(!is.na( y[[DVH]]$doseRx)))
doseAbs <- (valCmp/100) * y[[DVH]]$doseRx
ptMinDabs <- (ptMinD/100) * y[[DVH]]$doseRx
} else if(unitCmp %in% c("GY", "CGY")) {
## output is absolute dose -> nothing to do
doseAbs <- valCmp
ptMinDabs <- ptMinD
}
}
return(data.frame(doseAbs=doseAbs, volAbs=volAbs, volRel=volRel,
ptMinDabs=ptMinDabs,
ptMinVabs=ptMinVabs, ptMinVrel=ptMinVrel))
}
## get dose and (abs/rel) volume coords
doseVolL <- Map(getDoseVol,
DV=allDF$DV, valRef=allDF$valRef, unitRef=allDF$unitRef,
cmp=allDF$cmp, unitCmp=allDF$unitCmp, valCmp=allDF$valCmp,
DVH=if(byPat) { allDF$structure } else { allDF$patID },
ptMinD=allDF$ptMinD, ptMinV=allDF$ptMinV,
y=list(xConstrSub$x))
doseVolDF <- do.call("rbind", doseVolL)
## which volume should be plotted?
if(rel) {
doseVolDF$volPlot <- doseVolDF$volRel
doseVolDF$ptMinVplot <- doseVolDF$ptMinVrel
} else {
doseVolDF$volPlot <- doseVolDF$volAbs
doseVolDF$ptMinVplot <- doseVolDF$ptMinVabs
}
## add to existing data and remove invalid constraints
diagDF <- cbind(allDF, doseVolDF)
diagDF <- diagDF[diagDF$valid, ]
## create separate arrow in absolute size - point lower left and upper right
grLen <- 0.7
arrLL <- segmentsGrob(x0=grLen, y0=grLen, x1=0, y1=0, default.units="cm",
arrow=arrow(angle=45, length=unit(grLen, "cm")),
gp=gpar(lwd=3, lineend="butt", linejoin="mitre"))
arrUR <- segmentsGrob(x0=-grLen, y0=-grLen, x1=0, y1=0, default.units="cm",
arrow=arrow(angle=45, length=unit(grLen, "cm")),
gp=gpar(lwd=3, lineend="butt", linejoin="mitre"))
constraintArr <- function(dat) {
getArr <- function(x, y, cmp) {
if(grepl("^<.*", cmp)) {
annotation_custom(grobTree(arrLL), xmin=x, ymin=y, xmax=x+1, ymax=y+1)
} else if(grepl("^>.*", cmp)) {
annotation_custom(grobTree(arrUR), xmin=x, ymin=y, xmax=x+1, ymax=y+1)
}
}
annL <- Map(getArr, dat$doseAbs, dat$volPlot, dat$cmp)
}
## do the actual plotting
## get DVH plot for relevant structures / IDs
## constraint may be more extreme than actual doses
doseMax <- max(vapply(x, function(y) { max(y$dvh[ , "dose"]) }, numeric(1)))
if(isTRUE(any(diagDF$doseAbs > doseMax)) && (length(guessX) == 1L)) {
guessX <- max(diagDF$doseAbs)
}
diag <- showDVH(xConstrSub$x,
cumul=TRUE, byPat=byPat, rel=rel,
guessX=guessX, guessY=guessY,
thresh=thresh, show=FALSE)
## add constraint arrow to plot
diagC <- if(byPat) {
diag + constraintArr(diagDF) +
geom_point(data=diagDF,
aes_string(x="doseAbs", y="volPlot",
color="structure", shape="constraint"),
size=5) +
geom_point(data=diagDF,
aes_string(x="ptMinDabs", y="ptMinVplot",
color="structure", shape="constraint"),
size=5)
} else {
diag + constraintArr(diagDF) +
geom_point(data=diagDF,
aes_string(x="doseAbs", y="volPlot",
color="patID", shape="constraint"),
size=5) +
geom_point(data=diagDF,
aes_string(x="ptMinDabs", y="ptMinVplot",
color="patID", shape="constraint"),
size=5)
}
## for more than 6 constraints, we need more shapes
nConstr <- nlevels(as.factor(diagDF$constraint))
if(nConstr > 6) {
diagC <- diagC + scale_shape_manual(values=seq_len(nConstr))
}
if(show) {
print(diagC)
}
if(visible) {
diagC
} else {
invisible(diagC)
}
}
## x is a DVH list (1 component per id) of lists
showConstraint.DVHLstLst <-
function(x, constr, byPat=TRUE, rel=TRUE, guessX=TRUE, guessY=TRUE,
thresh=1, show=TRUE, visible=FALSE) {
## re-organize x into by-patient or by-structure form if necessary
isByPat <- attributes(x)$byPat
xRO <- if(is.null(isByPat) || (isByPat != byPat)) {
reorgByPat(x, byPat=byPat)
} else {
x
}
## make sure DVH and constraint have the same IDs / structures
xConstrSub <- harmoConstrDVH(xRO, constr=constr, byPat=byPat)
diagL <- Map(showConstraint,
x=xConstrSub$x, constr=xConstrSub$constr,
byPat=byPat, rel=rel, guessX=guessX, guessY=guessY,
thresh=thresh, show=show, visible=visible)
if(visible) {
diagL
} else {
invisible(diagL)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.