app.R
In DVHmetrics: Analyze Dose-Volume Histograms and Check Constraints

#####---------------------------------------------------------------------------
#####---------------------------------------------------------------------------
## DVHmetrics shiny App
## Daniel Wollschlaeger <wollschlaeger@uni-mainz.de>
#####---------------------------------------------------------------------------
#####---------------------------------------------------------------------------

library(shiny)
library(DVHmetrics)
library(bs4Dash)

source("global.R", encoding="UTF8")

## max number of IDs (=DVH files) we can make plots for
maxNID <- 100

#####---------------------------------------------------------------------------
## app code
#####---------------------------------------------------------------------------

shiny::shinyApp(
    #####-----------------------------------------------------------------------
    ## UI
    #####-----------------------------------------------------------------------
    ui=bs4DashPage(
        # theme = "custom.css",
        title="Analyze dose-volume histograms using DVHmetrics",
        sidebar_mini=TRUE,
        sidebar_collapsed=FALSE,
        sidebar=source("app_ui_sidebar.R", encoding="UTF8")$value,
        navbar=bs4DashNavbar(skin="light",
                             border=FALSE,
                             fixed=TRUE,
                             sidebarIcon="bars",
                             # controlbarIcon="th",
                             leftUi=tagList(
                                 tags$code(tags$h3("Analyze dose-volume histograms using 'DVHmetrics'"))
                             )
        ),
        body=bs4DashBody(
            bs4TabItems(
                bs4TabItem(
                    tabName="tab_data",
                    source("app_ui_tab_data.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_metrics",
                    source("app_ui_tab_metrics.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_show_dvh",
                    source("app_ui_tab_show_dvh.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_check_constraints",
                    source("app_ui_tab_check_constraints.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_show_constraints",
                    source("app_ui_tab_show_constraints.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_bed_eqd2",
                    source("app_ui_tab_bed_eqd2.R", local=TRUE, encoding="UTF8")$value
                ),
                bs4TabItem(
                    tabName="tab_about",
                    source("app_ui_tab_about.R", local=TRUE, encoding="UTF8")$value
                )
            )
        )
    ),
    #####-----------------------------------------------------------------------
    ## server
    #####-----------------------------------------------------------------------
    server=function(input, output, session) {
        session$onSessionEnded(stopApp)

        ## directory where current DVH data is saved
        DVHdir <- tempdir()

        ## reactive conductor
        DVH <- reactive({
            input$applyData
            isolate({
                DVHdata <- if(input$DVHin == '1') {
                    dataMZ
                } else if(input$DVHin == '2') {
                    if(!is.null(input$DVHupload)) {
                        types <- c('1'="Eclipse", '2'="Cadplan", '3'="Masterplan",
                                   '4'="Pinnacle", '5'="Monaco", '6'="HiArt",
                                   '7'="RayStation", '8'="ProSoma", '9'="PRIMO")
                        plans <- c('1'="none", '2'="doseRx")
                        # encodings <- c('1'="Default", '2'="UTF-8", '3'="UTF-8-BOM")
                        # enc <- encodings[input$fileenc]
                        # if(enc == "Default") { enc <- NULL }
                        argL  <- if(("DVHadd" %in% input$DVHreadOpts) && file.exists("DVHprev.rds")) {
                            setwd(DVHdir)
                            add <- readRDS("DVHprev.Rds")
                            list(x=input$DVHupload$datapath,
                                 type=types[input$DVHtype],
                                 # encoding=enc,
                                 courseAsID=("DVHcourse" %in% input$DVHreadOpts),
                                 # volume_from_dvh=("volume_from_dvh" %in% input$DVHreadOpts),
                                 uncertainty=("uncertainty" %in% input$DVHreadOpts),
                                 planInfo=plans[input$DVHplanInfo],
                                 add=add)
                        } else {
                            list(x=input$DVHupload$datapath,
                                 type=types[input$DVHtype],
                                 # encoding=enc,
                                 courseAsID=("DVHcourse" %in% input$DVHreadOpts),
                                 # volume_from_dvh=("volume_from_dvh" %in% input$DVHreadOpts),
                                 uncertainty=("uncertainty" %in% input$DVHreadOpts),
                                 planInfo=plans[input$DVHplanInfo])
                        }
                        # argL <- Filter(Negate(is.null), argL)
                        do.call(readDVH, argL)
                    } else {
                        NULL
                    }
                }

                ## add representation by structure
                DVHdataByStruct <- if(!is.null(DVHdata)) {
                    DVHmetrics:::reorgByPat(DVHdata, byPat=FALSE)
                } else {
                    NULL
                }

                DvHcurr <- list(DVH=DVHdata, DVHbyStruct=DVHdataByStruct)
                setwd(DVHdir)
                saveRDS(DvHcurr$DVH, file="DVHprev.rds")
                return(DvHcurr)
            })
        })

        output$DVHinfo <- renderPrint({
            dvh <- DVH()$DVH
            if(!is.null(dvh)) {
                if(input$DVHverbose == '1') {
                    print(dvh, verbose=FALSE)
                } else {
                    print(dvh, verbose=TRUE)
                }
            } else {
                invisible(NULL)
            }
        })

        output$metrSelStruct <- renderUI({
            dvh <- DVH()$DVH
            if(!is.null(dvh)) {
                checkboxGroupInput("metrSelStruct",
                                   label=h5("Select structures"),
                                   choices=getStrIDs(dvh, what="structure", choices=TRUE),
                                   selected=1)
            } else {
                NULL
            }
        })

        observeEvent(
            eventExpr=input$metrSelStructAll,
            handlerExpr= {
                if(input$metrSelStructAll == 0) { return(NULL) }
                dvh <- DVH()$DVH
                if(!is.null(dvh)) {
                    # browser()
                    choices <- getStrIDs(dvh, what="structure", choices=TRUE)
                    if((input$metrSelStructAll %% 2) == 0) {
                        updateCheckboxGroupInput(session, "metrSelStruct", "Select structures",
                                                 choices=choices)
                    } else {
                        updateCheckboxGroupInput(session, "metrSelStruct", "Select structures",
                                                 choices=choices,
                                                 selected=choices)
                    }
                }
            }
        )

        output$metrSelPat <- renderUI({
            dvh <- DVH()$DVH
            if(!is.null(dvh)) {
                checkboxGroupInput("metrSelPat",
                                   label=h5("Select patients"),
                                   choices=getStrIDs(dvh, what="patient", choices=TRUE),
                                   selected=1)
            } else {
                NULL
            }
        })

        ## https://stackoverflow.com/a/35044158
        observeEvent(
            eventExpr=input$metrSelPatAll,
            handlerExpr= {
                if(input$metrSelPatAll == 0) { return(NULL) }
                dvh <- DVH()$DVH
                if(!is.null(dvh)) {
                    # browser()
                    choices <- getStrIDs(dvh, what="patient", choices=TRUE)
                    if((input$metrSelPatAll %% 2) == 0) {
                        updateCheckboxGroupInput(session, "metrSelPat", "Select patients",
                                                 choices=choices)
                    } else {
                        updateCheckboxGroupInput(session, "metrSelPat", "Select patients",
                                                 choices=choices,
                                                 selected=choices)
                    }
                }
            }
        )

        output$plotSelStruct <- renderUI({
            dvh <- DVH()$DVH
            if(!is.null(dvh)) {
                checkboxGroupInput("plotSelStruct",
                                   label=h5("Select structures"),
                                   choices=getStrIDs(dvh, what="structure", choices=TRUE),
                                   selected=1)
            } else {
                NULL
            }
        })

        observeEvent(
            eventExpr=input$plotSelStructAll,
            handlerExpr= {
                if(input$plotSelStructAll == 0) { return(NULL) }
                dvh <- DVH()$DVH
                if(!is.null(dvh)) {
                    # browser()
                    choices <- getStrIDs(dvh, what="structure", choices=TRUE)
                    if((input$plotSelStructAll %% 2) == 0) {
                        updateCheckboxGroupInput(session, "plotSelStruct", "Select structures",
                                                 choices=choices)
                    } else {
                        updateCheckboxGroupInput(session, "plotSelStruct", "Select structures",
                                                 choices=choices,
                                                 selected=choices)
                    }
                }
            }
        )

        output$plotSelPat <- renderUI({
            dvh <- DVH()$DVH
            if(!is.null(dvh)) {
                checkboxGroupInput("plotSelPat",
                                   label=h5("Select patients"),
                                   choices=getStrIDs(dvh, what="patient", choices=TRUE),
                                   selected=1)
            } else {
                NULL
            }
        })

        observeEvent(
            eventExpr=input$plotSelPatAll,
            handlerExpr= {
                if(input$plotSelPatAll == 0) { return(NULL) }
                dvh <- DVH()$DVH
                if(!is.null(dvh)) {
                    # browser()
                    choices <- getStrIDs(dvh, what="patient", choices=TRUE)
                    if((input$plotSelPatAll %% 2) == 0) {
                        updateCheckboxGroupInput(session, "plotSelPat", "Select patients",
                                                 choices=choices)
                    } else {
                        updateCheckboxGroupInput(session, "plotSelPat", "Select patients",
                                                 choices=choices,
                                                 selected=choices)
                    }
                }
            }
        )

        output$metrics <- DT::renderDataTable({
            dvh        <- DVH()$DVH
            sortOpts   <- c('1'="observed", '2'="structure", '3'="metric", '4'="patID")
            splitOpts  <- c('1'="structure", '2'="metric", '3'="patID")
            selMetrics <- if(length(input$metrInput) > 0) {
                metrRaw <- unlist(strsplit(input$metrInput, "[[:blank:],]"))
                metrRaw[nzchar(metrRaw)]
            } else {
                c("DMEAN", "D1CC", "V10%")
            }
            selStruct <- if(length(input$metrSelStruct) > 0) {
                if(!is.null(dvh)) {
                    getStrIDs(dvh, what="structure")[as.numeric(input$metrSelStruct)]
                } else {
                    NULL
                }
            } else {
                NULL
            }
            selPat <- if(length(input$metrSelPat) > 0) {
                if(!is.null(dvh)) {
                    getStrIDs(dvh, what="patient")[as.numeric(input$metrSelPat)]
                } else {
                    NULL
                }
            } else {
                NULL
            }
            interp  <- "linear" # c("linear", "spline", "ksmooth")[as.numeric(input$metrInterp)]
            EUDa    <- if(input$metrEUDa  != "") { as.numeric(input$metrEUDa)  } else { NULL }
            EUDfd   <- if(input$metrEUDfd != "") { as.numeric(input$metrEUDfd) } else { NULL }
            EUDab   <- if(input$metrEUDab != "") { as.numeric(input$metrEUDab) } else { NULL }

            NTCPtype    <- c("probit", "logit", "poisson")[as.numeric(input$metrNTCPtype)]
            NTCPtd50    <- if(input$metrNTCPtd50    != "") { as.numeric(input$metrNTCPtd50)    } else { NULL }
            NTCPn       <- if(input$metrNTCPn       != "") { as.numeric(input$metrNTCPn)       } else { NULL }
            NTCPm       <- if(input$metrNTCPm       != "") { as.numeric(input$metrNTCPm)       } else { NULL }
            NTCPgamma50 <- if(input$metrNTCPgamma50 != "") { as.numeric(input$metrNTCPgamma50) } else { NULL }

            sortSel <- input$metrSortBy
            sortBy <- if(length(sortSel) > 0) {
                sortOpts[sortSel]
            } else {
                NULL
            }
            if(!is.null(dvh)) {
                argL <- list(x=dvh,
                             metric=selMetrics,
                             patID=selPat,
                             structure=selStruct,
                             sortBy=sortBy,
                             interp=interp,
                             EUDa=EUDa, EUDfd=EUDfd, EUDab=EUDab,
                             NTCPtype=NTCPtype, NTCPtd50=NTCPtd50, NTCPn=NTCPn, NTCPm=NTCPm, NTCPgamma50=NTCPgamma50,
                             TCPtype=NTCPtype, TCPtcd50=NTCPtd50,  TCPn=NTCPn,  TCPm=NTCPm,  TCPgamma50=NTCPgamma50)
                argL <- Filter(Negate(is.null), argL)
                metr <- do.call(getMetric, argL)
                metr$observed <- round(metr$observed, 2)
                metr
            } else {
                NULL
            }
        })#, options=list(pageLength=25))

        output$saveMetrics <- downloadHandler(
            filename=function() { "metrics.txt" },
            content=function(file) {
                dvh       <- DVH()$DVH
                sortOpts  <- c('1'="observed", '2'="structure", '3'="metric", '4'="patID")
                splitOpts <- c('1'="metric", '2'="structure", '3'="patID")
                selMetrics <- if(length(input$metrInput) > 0) {
                    metrRaw <- unlist(strsplit(input$metrInput, "[[:blank:],]"))
                    metrRaw[nzchar(metrRaw)]
                } else {
                    "DMEAN"
                }
                selPat <- if(length(input$metrSelPat) > 0) {
                    if(!is.null(dvh)) {
                        getStrIDs(dvh, what="patient")[as.numeric(input$metrSelPat)]
                    } else {
                        NULL
                    }
                } else {
                    NULL
                }
                selStruct <- if(length(input$metrSelStruct) > 0) {
                    if(!is.null(dvh)) {
                        getStrIDs(dvh, what="structure")[as.numeric(input$metrSelStruct)]
                    } else {
                        NULL
                    }
                } else {
                    NULL
                }

                interp  <- "linear"
                EUDa    <- if(input$metrEUDa  != "") { as.numeric(input$metrEUDa)  } else { NULL }
                EUDfd   <- if(input$metrEUDfd != "") { as.numeric(input$metrEUDfd) } else { NULL }
                EUDab   <- if(input$metrEUDab != "") { as.numeric(input$metrEUDab) } else { NULL }

                NTCPtype    <- c("probit", "logit", "poisson")[as.numeric(input$metrNTCPtype)]
                NTCPtd50    <- if(input$metrNTCPtd50    != "") { as.numeric(input$metrNTCPtd50)    } else { NULL }
                NTCPn       <- if(input$metrNTCPn       != "") { as.numeric(input$metrNTCPn)       } else { NULL }
                NTCPm       <- if(input$metrNTCPm       != "") { as.numeric(input$metrNTCPm)       } else { NULL }
                NTCPgamma50 <- if(input$metrNTCPgamma50 != "") { as.numeric(input$metrNTCPgamma50) } else { NULL }

                sortSel <- input$metrSortBy
                sortBy <- if(length(sortSel) > 0) {
                    sortOpts[sortSel]
                } else {
                    NULL
                }
                argL <- list(x=dvh,
                             metric=selMetrics,
                             patID=selPat,
                             structure=selStruct,
                             sortBy=sortBy,
                             interp=interp,
                             EUDa=EUDa, EUDfd=EUDfd, EUDab=EUDab,
                             NTCPtype=NTCPtype, NTCPtd50=NTCPtd50, NTCPn=NTCPn, NTCPm=NTCPm, NTCPgamma50=NTCPgamma50,
                             TCPtype=NTCPtype, TCPtcd50=NTCPtd50,  TCPn=NTCPn,  TCPm=NTCPm,  TCPgamma50=NTCPgamma50)
                argL <- Filter(Negate(is.null), argL)
                metr <- do.call(getMetric, argL)
                dec <- c('1'=".",  '2'=",")[input$saveMetrDec]
                sep <- c('1'="\t", '2'=" ", '3'=",", '4'=";")[input$saveMetrSep]
                write.table(metr, file=file, dec=dec, sep=sep, row.names=FALSE)
            },
            contentType='text/plain' # MIME type
        )

        #     output$DVHplotOrg <- renderPlot({
        #         dvh <- DVH()$DVH
        #         if(!is.null(dvh)) {
        #             argL <- list(x=dvh,
        #                          patID=    getStrIDs(dvh, what="patient")[  as.numeric(input$plotSelPat)],
        #                          structure=getStrIDs(dvh, what="structure")[as.numeric(input$plotSelStruct)],
        #                          thresh=input$plotThreshVol,
        #                          rel=input$plotPlotVol == '1')
        #             do.call(showDVH, argL)
        #         } else {
        #             NULL
        #         }
        #     })

        output$DVHplot <- renderUI({
            dvh       <- DVH()
            byPat     <- input$plotByPat == '1'
            selPat    <- getStrIDs(dvh$DVH, what="patient")[  as.numeric(input$plotSelPat)]
            selStruct <- getStrIDs(dvh$DVH, what="structure")[as.numeric(input$plotSelStruct)]
            plotOutputL <- if(byPat) {
                ## 1 diagram per patient
                lapply(seq_along(dvh$DVH), function(i) {
                    iMax <- sum(selPat %in% names(dvh$DVH))
                    if(i <= iMax) {    # avoid creating empty plots
                        plotName <- paste0("DVHplot", i)
                        plotOutput(plotName,
                                   dblclick=paste0(plotName, "_dblclick"),
                                   brush=brushOpts(id=paste0(plotName, "_brush"),
                                                   clip=TRUE,
                                                   resetOnNew=FALSE))
                    } else {
                        NULL
                    }
                })
            } else {
                ## 1 diagram per structure
                lapply(seq_along(dvh$DVHbyStruct), function(i) {
                    iMax <- sum(selStruct %in% names(dvh$DVHbyStruct))
                    if(i <= iMax) {    # avoid creating empty plots
                        plotName <- paste0("DVHplot", i)
                        plotOutput(plotName,
                                   dblclick=paste0(plotName, "_dblclick"),
                                   brush=brushOpts(id=paste0(plotName, "_brush"),
                                                   clip=TRUE,
                                                   resetOnNew=FALSE))
                    } else {
                        NULL
                    }
                })
            }

            ## weed out NULL components, convert the list to a tagList and return
            plotOutputL <- Filter(Negate(is.null), plotOutputL)
            do.call(tagList, plotOutputL)
        })

        ## adapted from Winston Chang: https://gist.github.com/wch/5436415
        ## call renderPlot maxNID times and discard those not required (return NULL)
        ## problem: number of plots should be dynamic -> reactive context required
        for(i in seq_len(maxNID)) {
            ## need local so that each item gets its own number. Without it, the value
            ## of i in the renderPlot() will be the same across all instances, because
            ## of when the expression is evaluated.
            local({
                localI <- i
                ## DVH plot
                output[[paste0("DVHplot", localI)]] <- renderPlot({
                    dvh   <- DVH()
                    rel   <- input$plotPlotVol == '1'
                    cumul <- input$plotType    == '1'
                    byPat <- input$plotByPat   == '1'
                    selPat    <- getStrIDs(dvh$DVH, what="patient")[  as.numeric(input$plotSelPat)]
                    selStruct <- getStrIDs(dvh$DVH, what="structure")[as.numeric(input$plotSelStruct)]
                    if(( byPat && (sum(selPat    %in% names(dvh$DVH))         < localI)) ||
                       (!byPat && (sum(selStruct %in% names(dvh$DVHbyStruct)) < localI))) {
                        NULL
                    } else {
                        ## 1 diagram per patient/structure
                        x <- if(byPat) {
                            sharedNames <- intersect(selPat, names(dvh$DVH))
                            dvh$DVH[sharedNames][[localI]]
                        } else {
                            sharedNames <- intersect(selStruct, names(dvh$DVHbyStruct))
                            dvh$DVHbyStruct[sharedNames][[localI]]
                        }

                        brush <- input[[paste0("DVHplot", localI, "_brush")]]
                        if(!is.null(brush)) {
                            showDVH(x=x,
                                    cumul=cumul,
                                    byPat=byPat,
                                    patID=selPat,
                                    structure=selStruct,
                                    guessX=c(brush$xmin, brush$xmax),
                                    guessY=c(brush$ymin, brush$ymax),
                                    thresh=input$plotThreshVol,
                                    rel=rel,
                                    addMSD=input$plotMSD,
                                    show=FALSE,
                                    visible=TRUE)
                        } else {
                            showDVH(x=x,
                                    cumul=cumul,
                                    byPat=byPat,
                                    patID=selPat,
                                    structure=selStruct,
                                    thresh=input$plotThreshVol,
                                    rel=rel,
                                    addMSD=input$plotMSD,
                                    show=FALSE,
                                    visible=TRUE)
                        }
                    }
                })

                ## constraint plot
                output[[paste0("constraintPlot", localI)]] <- renderPlot({
                    dvh    <- DVH()
                    constr <- DVHconstr()
                    rel    <- input$constrPlotVol == "1"
                    byPat  <- input$constrByPat   == "1"
                    ## 1 diagram per patient/structure
                    ## restrict DVH and constr to the same IDs/structures
                    xConstrSub <- if(byPat) {
                        DVHmetrics:::harmoConstrDVH.DVHLstLst(dvh$DVH,         constr=constr, byPat=byPat)
                    } else {
                        DVHmetrics:::harmoConstrDVH.DVHLstLst(dvh$DVHbyStruct, constr=constr, byPat=byPat)
                    }

                    if(is.null(constr) || (length(xConstrSub$x) < localI)) {
                        NULL
                    } else {
                        brush <- input[[paste0("constraintPlot", localI, "_brush")]]
                        if(!is.null(brush)) {
                            showConstraint(x=xConstrSub$x[[localI]],
                                           constr=xConstrSub$constr[[localI]],
                                           byPat=byPat,
                                           thresh=input$constrThreshVol,
                                           guessX=c(brush$xmin, brush$xmax),
                                           guessY=c(brush$ymin, brush$ymax),
                                           rel=rel,
                                           show=FALSE,
                                           visible=TRUE)
                        } else {
                            showConstraint(x=xConstrSub$x[[localI]],
                                           constr=xConstrSub$constr[[localI]],
                                           byPat=byPat,
                                           thresh=input$constrThreshVol,
                                           rel=rel,
                                           show=FALSE,
                                           visible=TRUE)
                        }
                    }
                })
            })
        }

        output$saveDVHPDF <- downloadHandler(
            filename=function() { "DVH.pdf" },
            content=function(file) {
                dvh   <- DVH()
                byPat <- input$plotByPat   == '1'
                rel   <- input$plotPlotVol == "1"
                cumul <- input$plotType    == "1"
                selPat    <- getStrIDs(dvh$DVH, what="patient")[  as.numeric(input$plotSelPat)]
                selStruct <- getStrIDs(dvh$DVH, what="structure")[as.numeric(input$plotSelStruct)]
                x <- if(byPat) {
                    dvh$DVH
                } else {
                    dvh$DVHbyStruct
                }

                argL <- list(x=x,
                             cumul=cumul,
                             byPat=byPat,
                             patID=selPat,
                             structure=selStruct,
                             thresh=input$plotThreshVol,
                             rel=rel,
                             addMSD=input$plotMSD)
                pdf(file, width=7, height=5)
                do.call(showDVH, argL)
                dev.off()
            },
            contentType='application/pdf' # MIME type
        )

        makeDVHJPG <- function(fName, dvh, byPat, rel, cumul, selPat, selStruct, thresh) {
            argL <- list(x=dvh,
                         cumul=cumul,
                         byPat=byPat,
                         patID=selPat,
                         structure=selStruct,
                         thresh=thresh,
                         rel=rel)
            jpeg(fName, quality=100, width=700, height=500)
            do.call(showDVH, argL)
            dev.off()
        }

        output$saveDVHJPG <- downloadHandler(
            filename=function() { "DVH_JPGs.zip" },
            content=function(file) {
                dvh    <- DVH()
                byPat  <- input$plotByPat   == '1'
                rel    <- input$plotPlotVol == "1"
                cumul  <- input$plotType    == "1"
                thresh <- input$plotThreshVol
                addMSD <- input$plotMSD
                x <- if(byPat) {
                    dvh$DVH
                } else {
                    dvh$DVHbyStruct
                }

                selPat    <- getStrIDs(dvh$DVH, what="patient")[  as.numeric(input$plotSelPat)]
                selStruct <- getStrIDs(dvh$DVH, what="structure")[as.numeric(input$plotSelStruct)]
                nFiles <- if(byPat) {
                    sum(selPat    %in% names(x))
                } else {
                    sum(selStruct %in% names(x))
                }

                ## write all JPEGs into temporary directory and zip them
                fNames <- paste0("DVH_", sprintf("%02d", seq_len(nFiles)), ".jpg")
                tmpdir <- tempdir()
                setwd(tmpdir)

                for(i in seq_along(fNames)) {
                    fN         <- fNames[i]
                    selPatI    <- if(byPat) { selPat[i] } else { selPat }
                    selStructI <- if(byPat) { selStruct } else { selStruct[i] }
                    makeDVHJPG(fN, dvh=x, byPat=byPat, rel=rel, cumul=cumul,
                               selPat=selPatI, selStruct=selStructI,
                               addMSD=addMSD, thresh=thresh)
                }

                zip(zipfile=file, files=fNames)
                if(file.exists(paste0(file, ".zip"))) {
                    file.rename(paste0(file, ".zip"), file)
                }
            },
            contentType = "application/zip"
        )

        output$saveDVHMSD <- downloadHandler(
            filename=function() { "DVH_M_SD.txt" },
            content=function(file) {
                dvh    <- DVH()
                byPat  <- input$plotByPat   == '1'
                cumul  <- input$plotType    == '1'
                x <- if(byPat) {
                    dvh$DVH
                } else {
                    dvh$DVHbyStruct
                }

                selPat    <- getStrIDs(dvh$DVH, what="patient")[  as.numeric(input$plotSelPat)]
                selStruct <- getStrIDs(dvh$DVH, what="structure")[as.numeric(input$plotSelStruct)]

                argL <- list(x=x,
                             fun=list(mean=mean, median=median, min=min, max=max, sd=sd),
                             byPat=byPat,
                             thin=1,
                             cumul=cumul,
                             patID=selPat,
                             structure=selStruct,
                             fixed=TRUE)
                argL   <- Filter(Negate(is.null), argL)
                DVHMSD <- do.call(getMeanDVH, argL)
                dec    <- c('1'=".",  '2'=",")[input$saveDVHDec]
                sep    <- c('1'="\t", '2'=" ", '3'=",", '4'=";")[input$saveDVHSep]
                write.table(DVHMSD, file=file, dec=dec, sep=sep, row.names=FALSE)
            },
            contentType='text/plain' # MIME type
        )

        ## reactive conductor
        DVHconstr <- reactive({
            input$applyConstraints
            isolate({
                constr <- if(input$constrIn == '1') {
                    dataConstr
                } else if(input$constrIn == '2') {
                    if(is.null(input$constrUpload)) {
                        NULL
                    } else {
                        dec <- c('1'=".",  '2'=",")[input$constrDec]
                        sep <- c('1'="\t", '2'=" ", '3'=",", '4'=";")[input$constrSep]
                        ## peek at data to see if we have a header
                        line1 <- readLines(input$constrUpload$datapath, n=1)
                        if(grepl("constraint", line1, ignore.case=TRUE)) {
                            read.table(file=input$constrUpload$datapath,
                                       header=TRUE, sep=sep, dec=dec, stringsAsFactors=FALSE)
                        } else {
                            read.table(file=input$constrUpload$datapath,
                                       header=FALSE, sep=sep, dec=dec, stringsAsFactors=FALSE)
                        }
                    }
                } else if(input$constrIn == '3') {
                    dec <- c('1'=".",  '2'=",")[input$constrPasteDec]
                    sep <- c('1'="\t", '2'=" ", '3'=",", '4'=";")[input$constrPasteSep]
                    ## peek at data to see if we have a header
                    line1 <- readLines(textConnection(input$constrPaste), n=1)
                    if(grepl("constraint", line1, ignore.case=TRUE)) {
                        read.table(file=textConnection(input$constrPaste),
                                   header=TRUE, sep=sep, dec=dec, stringsAsFactors=FALSE)
                    } else {
                        read.table(file=textConnection(input$constrPaste),
                                   header=FALSE, sep=sep, dec=dec, stringsAsFactors=FALSE)
                    }
                } else {
                    NULL
                }

                if(!is.null(constr)) {
                    cNames <- tolower(names(constr))
                    cNames[cNames == "patid"] <- "patID"
                    constr <- setNames(constr, cNames)
                }

                return(constr)
            })
        })

        output$constraints <- DT::renderDataTable({
            constr <- DVHconstr()
            dvh    <- DVH()$DVH
            outSel <- constrOutInv[input$constrOut]
            interp <- "linear" #c("linear", "spline", "ksmooth")[as.numeric(input$constrInterp)]
            EUDa   <- if(input$constrEUDa  != "") { as.numeric(input$constrEUDa)  } else { NULL }
            EUDfd  <- if(input$constrEUDfd != "") { as.numeric(input$constrEUDfd) } else { NULL }
            EUDab  <- if(input$constrEUDab != "") { as.numeric(input$constrEUDab) } else { NULL }

            NTCPtype    <- c("probit", "logit", "poisson")[as.numeric(input$constrNTCPtype)]
            NTCPtd50    <- if(input$constrNTCPtd50    != "") { as.numeric(input$constrNTCPtd50)    } else { NULL }
            NTCPn       <- if(input$constrNTCPn       != "") { as.numeric(input$constrNTCPn)       } else { NULL }
            NTCPm       <- if(input$constrNTCPm       != "") { as.numeric(input$constrNTCPm)       } else { NULL }
            NTCPgamma50 <- if(input$constrNTCPgamma50 != "") { as.numeric(input$constrNTCPgamma50) } else { NULL }

            sortOpts <- c('1'="compliance", '2'="dstMin", '3'="deltaV", '4'="deltaD",
                          '5'="observed", '6'="patID", '7'="structure", '8'="constraint")
            sortSel  <- input$constrSortBy
            sortBy   <- if(length(sortSel) > 0) {
                sortOpts[sortSel]
            } else {
                "none"
            }

            if(!is.null(constr) && !is.null(dvh)) {
                argL <- list(x=dvh,
                             constr=constr, byPat=TRUE, interp=interp,
                             semSign=input$constrSemSign,
                             sortBy=sortBy,
                             EUDa=EUDa, EUDfd=EUDfd, EUDab=EUDab,
                             NTCPtype=NTCPtype, NTCPtd50=NTCPtd50, NTCPn=NTCPn, NTCPm=NTCPm, NTCPgamma50=NTCPgamma50,
                             TCPtype=NTCPtype, TCPtcd50=NTCPtd50,  TCPn=NTCPn,  TCPm=NTCPm,  TCPgamma50=NTCPgamma50)
                argL <- Filter(Negate(is.null), argL)
                x <- do.call(checkConstraint, argL)
                x$observed  <- round(x$observed,  2)
                x$deltaVpc  <- round(x$deltaVpc,  2)
                x$deltaDpc  <- round(x$deltaDpc,  2)
                x$dstMin    <- round(x$dstMin,    2)
                x$dstMinRel <- round(x$dstMinRel, 2)
                x[ , outSel]
            } else {
                NULL
            }

        })#, options=list(pageLength=25))

        #     output$constraintPlotOrg <- renderPlot({
        #         constr <- DVHconstr()
        #         dvh    <- DVH()$DVH
        #         if(!is.null(constr) && !is.null(dvh)) {
        #             argL <- list(x=dvh,
        #                          thresh=input$constrThreshVol,
        #                          rel=input$constrPlotVol == '1',
        #                          constr=constr)
        #             do.call(showConstraint, argL)
        #         } else {
        #             NULL
        #         }
        #     })

        output$constraintPlot <- renderUI({
            dvh    <- DVH()
            constr <- DVHconstr()
            byPat  <- input$constrByPat == "1"
            xConstrSub <- if(byPat) {
                DVHmetrics:::harmoConstrDVH.DVHLstLst(dvh$DVH,         constr=constr, byPat=byPat)
            } else {
                DVHmetrics:::harmoConstrDVH.DVHLstLst(dvh$DVHbyStruct, constr=constr, byPat=byPat)
            }

            plotOutputL <- lapply(seq_along(xConstrSub$x), function(i) {
                plotName <- paste0("constraintPlot", i)
                plotOutput(plotName,
                           dblclick=paste0(plotName, "_dblclick"),
                           brush=brushOpts(id=paste0(plotName, "_brush"),
                                           clip=TRUE,
                                           resetOnNew=FALSE))
            })

            ## convert the list to a tagList and return
            do.call(tagList, plotOutputL)
        })

        output$saveConstrTxt <- downloadHandler(
            filename=function() { "constraints.txt" },
            content=function(file) {
                interp <-"linear"
                EUDa   <- if(input$constrEUDa  != "") { as.numeric(input$constrEUDa)  } else { NULL }
                EUDfd  <- if(input$constrEUDfd != "") { as.numeric(input$constrEUDfd) } else { NULL }
                EUDab  <- if(input$constrEUDab != "") { as.numeric(input$constrEUDab) } else { NULL }

                NTCPtype    <- c("probit", "logit", "poisson")[as.numeric(input$constrNTCPtype)]
                NTCPtd50    <- if(input$constrNTCPtd50    != "") { as.numeric(input$constrNTCPtd50)    } else { NULL }
                NTCPn       <- if(input$constrNTCPn       != "") { as.numeric(input$constrNTCPn)       } else { NULL }
                NTCPm       <- if(input$constrNTCPm       != "") { as.numeric(input$constrNTCPm)       } else { NULL }
                NTCPgamma50 <- if(input$constrNTCPgamma50 != "") { as.numeric(input$constrNTCPgamma50) } else { NULL }

                argL <- list(x=DVH()$DVH,
                             constr=DVHconstr(), interp=interp,
                             EUDa=EUDa, EUDfd=EUDfd, EUDab=EUDab,
                             NTCPtype=NTCPtype, NTCPtd50=NTCPtd50, NTCPn=NTCPn, NTCPm=NTCPm, NTCPgamma50=NTCPgamma50,
                             TCPtype=NTCPtype, TCPtcd50=NTCPtd50,  TCPn=NTCPn,  TCPm=NTCPm,  TCPgamma50=NTCPgamma50)
                argL   <- Filter(Negate(is.null), argL)
                constr <- do.call(checkConstraint, argL)
                dec    <- c('1'=".",  '2'=",")[input$saveConstrDec]
                sep    <- c('1'="\t", '2'=" ", '3'=",", '4'=";")[input$saveConstrSep]
                write.table(constr, file=file, dec=dec, sep=sep, row.names=FALSE)
            },
            contentType='text/plain' # MIME type
        )

        output$saveConstrPDF <- downloadHandler(
            filename=function() { "constraints.pdf" },
            content=function(file) {
                rel   <- input$constrPlotVol == '1'
                byPat <- input$constrByPat   == '1'
                pdf(file, width=7, height=5)
                showConstraint(x=DVH()$DVH,
                               constr=DVHconstr(),
                               byPat=byPat,
                               thresh=input$constrThreshVol,
                               rel=rel)
                dev.off()
            },
            contentType='application/pdf' # MIME type
        )

        makeConstrJPG <- function(fName, dvh, constr, byPat, thresh, rel) {
            argL <- list(x=dvh,
                         constr=constr,
                         byPat=byPat,
                         thresh=thresh,
                         rel=rel)
            jpeg(fName, quality=100, width=700, height=500)
            do.call(showConstraint, argL)
            dev.off()
        }

        output$saveConstrJPG <- downloadHandler(
            filename=function() { "constraints.zip" },
            content=function(file) {
                dvh    <- DVH()
                constr <- DVHconstr()
                rel    <- input$constrPlotVol == '1'
                byPat  <- input$constrByPat   == '1'
                thresh <- input$constrThreshVol
                x <- if(byPat) {
                    dvh$DVH
                } else {
                    dvh$DVHbyStruct
                }

                xConstrSub <- DVHmetrics:::harmoConstrDVH.DVHLstLst(x, constr=constr, byPat=byPat)
                nFiles     <- length(xConstrSub$x)

                ## write all JPEGs into temporary directory and zip them
                fNames <- paste0("DVH_", sprintf("%02d", seq_len(nFiles)), ".jpg")
                tmpdir <- tempdir()
                setwd(tmpdir)

                for(i in seq_along(fNames)) {
                    fN <- fNames[i]
                    makeConstrJPG(fN, dvh=xConstrSub$x, constr=xConstrSub$constr,
                                  byPat=byPat, rel=rel, thresh=thresh)
                }

                zip(zipfile=file, files=fNames)
                if(file.exists(paste0(file, ".zip"))) {
                    file.rename(paste0(file, ".zip"), file)
                }
            },
            contentType = "application/zip"
        )

        output$BED <- renderPrint({
            if(input$BEDtype %in% c('1', '2')) {
                D <- if(input$BED_BED_D  != "") {
                    as.numeric(strsplit(input$BED_BED_D, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fd <- if(input$BED_BED_FD != "") {
                    as.numeric(strsplit(input$BED_BED_FD, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fn <- if(input$BED_BED_FN != "") {
                    as.numeric(strsplit(input$BED_BED_FN, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                ab <- if(input$BED_BED_AB != "") {
                    as.numeric(strsplit(input$BED_BED_AB, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
            } else if(input$BEDtype == '3') {
                D1 <- if(input$BED_IED_D1  != "") {
                    as.numeric(strsplit(input$BED_IED_D1, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                D2 <- if(input$BED_IED_D2  != "") {
                    as.numeric(strsplit(input$BED_IED_D2, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fd1 <- if(input$BED_IED_FD1 != "") {
                    as.numeric(strsplit(input$BED_IED_FD1, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fd2 <- if(input$BED_IED_FD2 != "") {
                    as.numeric(strsplit(input$BED_IED_FD2, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fn1 <- if(input$BED_IED_FN1 != "") {
                    as.numeric(strsplit(input$BED_IED_FN1, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                fn2 <- if(input$BED_IED_FN2 != "") {
                    as.numeric(strsplit(input$BED_IED_FN2, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
                ab <- if(input$BED_IED_AB != "") {
                    as.numeric(strsplit(input$BED_IED_AB, "[[:blank:]]")[[1]])
                } else {
                    NULL
                }
            }

            if(input$BEDtype == '1') {
                getBED(D=D, fd=fd, fn=fn, ab=ab)
            } else if(input$BEDtype == '2') {
                getEQD2(D=D, fd=fd, fn=fn, ab=ab)
            } else if(input$BEDtype == '3') {
                getIsoEffD(D1=D1, D2=D2, fd1=fd1, fd2=fd2, ab=ab)
            }
        })
    }
)