R/plot_regist_fda.R

In snandi/Registration: Contains functions for registration algorithms, simulation and testing

#########################################################################################
## Plot Registration output from fda, curve by curve
#########################################################################################
# geom_line linetypes
# 0 = blank, 1 = solid, 2 = dashed, 3 = dotted, 4 = dotdash, 5 = longdash, 6 = twodash

# registOutput <- Regfd.M
# TrueWarp = NULL
# Lambda = NULL
# Ylabel = NULL
# Xlabel = NULL
#' @import fda
#' @importFrom gridExtra grid.arrange
#' @import ggplot2
#' @importFrom RFunctionsSN rowSE
#' @importFrom graphics plot par
#' 
#@import fdasrvf

plotbyCurve_regist_fda <- function(
  registOutput,
  TrueWarp = NULL,
  Lambda = NULL,
  Ylabel = NULL,
  Xlabel = NULL
) {
  #   require(ggplot2)
  #   require(gridExtra)
  #   require(fdasrvf)
  #  plots the output of registration function register.fd, one curve at a time
  #
  #  Argument:
  #  registOutput ... The named list generated by function REGISTER.FD
  #              The members of registOutput that are required are:
  #              REGFD  ... the registered functions
  #              WFD    ... the functions W(t) defining the warping functions h(t)
  #              YFD    ... the unregistered functions
  #              Y0FD   ... the target functions
  #  If required objects are missing, registOutput was probably generated by
  #  an older verson of REGISTERFD, and the registration should be redone.
  
  # TrueWarp: is a matrix of the true warping function. This is mainly for 
  #           comparisons of true warping functions vs estimated warping
  #           functions and used in simulations only
  
  #  Last modified  by S. Nandi on July 2, 2015
  
  #  check the argument registOutput
  if (!inherits(registOutput, "list")) stop("registOutput is not a list object.")
  
  #  check the argument TrueWarp
  if (!is.null(TrueWarp)) {
    if (!(is.matrix(TrueWarp))) TrueWarp <- as.matrix(TrueWarp)
    if (!(is.matrix(TrueWarp))) stop('TrueWarp needs to be in matrix format')
  }
  
  #  extract the required members of registOutput
  missinginfo <- FALSE
  if (is.null(registOutput$regfd)) {
    missinginfo <- TRUE
  } else {
    yregfd  <- registOutput$regfd
  }
  if (is.null(registOutput$Wfd)) {
    missinginfo <- TRUE
  } else {
    Wfd     <- registOutput$Wfd
  }
  if (is.null(registOutput$yfd)) {
    missinginfo <- TRUE
  } else {
    yfd     <- registOutput$yfd
  }
  if (is.null(registOutput$y0fd)) {
    missinginfo <- TRUE
  } else {
    y0fd    <- registOutput$y0fd
  }
  if (missinginfo) stop("registOutput does not containing required objects.")
  
  #  generate a fine mesh of argument values
  ybasis  <- yfd$basis
  yrange  <- ybasis$rangeval
  ynbasis <- ybasis$nbasis
  nfine   <- 201
  argfine <- seq(yrange[1],yrange[2],len=nfine)
  
  if (!is.null(TrueWarp)) argOrig <- seq(yrange[1],yrange[2],len = nrow(TrueWarp))  
  
  #  evaluate the required functional on over this fine mesh
  ymat    <- eval.fd(argfine, yfd)
  y0mat   <- eval.fd(argfine, y0fd)
  yregmat <- eval.fd(argfine, yregfd)
  warpmat <- eval.monfd(argfine, Wfd)
  warpmat <- yrange[1] + (yrange[2]-yrange[1])*warpmat/
    (matrix(1,nfine,1)%*%warpmat[nfine,])
  
  #  extract the number of functions NCURVE and the number of variables NVAR
  ydim     <- dim(yfd$coef)
  ncurve   <- ydim[2]
  
  if (!is.null(names(yfd$fdname)[[3]])) {
    ylabel = names(yfd$fdname)[[3]]
  }  else {
    ylabel = "Function value"
  }
  casename <- yfd$fdname[[2]]
  if (length(casename) != ncurve) casename = as.character(1:ncurve)
  
  if (dim(y0mat)[2] == 1) y0mat = y0mat %*% matrix(1,1,ncurve)
  ylimit   <- c(min(ymat),max(ymat))
  
  graphics::par(mfrow=c(1,2),pty = "s")
  i <- 1
  for (i in 1:ncurve) {
    ## ymat: original; 
    ## y0mat: mean to be registered to
    ## yregmat: registered curve
    
    Data.Orig <- as.data.frame(cbind(argfine, ymat = ymat[,i]))
    Data.Mean <- as.data.frame(cbind(argfine, y0mat = y0mat[,i]))
    Data.Reg <- as.data.frame(cbind(argfine, yregmat = yregmat[,i]))
    if(is.null(Xlabel)) Xlabel <- 'Pixel Position'
    if(is.null(Ylabel)) Ylabel <- 'Intensity'
    MainTitle <- paste("Curve Number",casename[i])
    Col.Orig <- 'deepskyblue'
    Col.Mean <- 'orange'
    Col.Reg <- 'springgreen'
    
    PlotCurve <- ggplot()+
      ylim(ylimit) +
      geom_line(data = Data.Orig, aes(x = argfine, y = ymat, col = 'Original'), linetype = 2, size = 1) +
      geom_line(data = Data.Mean, aes(x = argfine, y = y0mat, col = 'Mean'), linetype = 1, size = 1) + 
      geom_line(data = Data.Reg, aes(x = argfine, y = yregmat, col = 'Registered'), linetype = 1, size = 1) +
      xlab(label = Xlabel) + ylab(label = Ylabel) +
      ggtitle(MainTitle) +
      scale_colour_manual(name = '', values = c('Original' = Col.Orig, 'Mean' = Col.Mean, 'Registered' = Col.Reg),
                          breaks = c("Original", "Mean", "Registered")) +
      theme(plot.title = element_text(face = "bold", size = 12, colour = "white"),
            #             panel.background = element_rect(fill = 'black'), 
            #             plot.background = element_rect(color = 'black', fill = "gray10"), 
            panel.background = element_rect(fill = 'gray60'), 
            plot.background = element_rect(color = 'gray60', fill = "gray60"), 
            axis.text = element_text(colour = "white", size = 10), 
            axis.title.x = element_text(colour = "white", size = 10), 
            axis.title.y = element_text(colour = "white", size = 10), 
            panel.grid.major = element_line(colour = "gray30", size = 0.35), 
            panel.grid.minor = element_line(colour = "gray20", size = 0.25),
            legend.position = c(0.5,0.9), 
            legend.direction = 'horizontal', 
            legend.text = element_text(size = 10, colour = 'white'),
            #             legend.background = element_rect(fill = 'black')
            legend.background = element_rect(fill = 'gray60')
      )
    
    #PlotCurve
    
    Data.Warp <- as.data.frame(cbind(argfine, Warp = warpmat[,i]))
    if (is.null(TrueWarp)){
      PlotWarp <- ggplot()+
        geom_line(data = Data.Warp, aes(x = argfine, y = Warp, col = 'Warping fn'), linetype = 1, 
                  size = 1, colour = Col.Reg) +
        geom_abline(intercept = 0, slope = 1, col = 'white') + 
        xlab(label = '') + ylab(label = '') +
        ggtitle('Warping Functions, Using Min Eig Value') +
        theme(
          plot.title = element_text(face = "bold", size = 12, colour = "white"),
          panel.background = element_rect(fill = 'black'), 
          plot.background = element_rect(color = 'black', fill  =  "gray10"), 
          axis.text = element_text(colour = "white", size = 10), 
          panel.grid.major = element_line(colour = "gray30", size = 0.35), 
          panel.grid.minor = element_line(colour = "gray20", size = 0.25)
        )
    } else{
      Data.Warp <- cbind(Data.Warp, TrueWarp = TrueWarp[,i])
      ###???!!! Add the code with TrueWarp !!!???###
    }
    grid.arrange(PlotCurve, PlotWarp, ncol = 1)  
  }
  return(yregmat)
}

#########################################################################################
## Plot Registration output from fda, all curves together
#########################################################################################
# registOutput <- Regfd.M
# TrueWarp = NULL
# Lambda = NULL
# Ylabel = NULL
# Xlabel = NULL
# This function currently prints the plots it generates. It should be wrapped inside a
# pdf object call. 
plotAll_regist_fda <- function(
  registOutput,
  TrueWarp = NULL,
  Lambda = NULL,
  TitleText = "",
  Ylabel = NULL,
  Xlabel = NULL,
  BeforeAfterDistance = TRUE,
  PlotBeforeRegist = TRUE,
  Xarg_fine = NULL,
  saveToPDF = FALSE
) {
  #   require(ggplot2)
  #   require(gridExtra)
  #   require(reshape2)
  #  plots the output of registration function register.fd, one curve at a time
  #
  #  Argument:
  #  registOutput ... The named list generated by function REGISTER.FD
  #              The members of registOutput that are required are:
  #              REGFD  ... the registered functions
  #              WFD    ... the functions W(t) defining the warping functions h(t)
  #              YFD    ... the unregistered functions
  #              Y0FD   ... the target functions
  #  If required objects are missing, registOutput was probably generated by
  #  an older verson of REGISTERFD, and the registration should be redone.
  
  # TrueWarp: is a matrix of the true warping function. This is mainly for 
  #           comparisons of true warping functions vs estimated warping
  #           functions and used in simulations only
  
  #  Last modified  by S. Nandi on July 2, 2015
  
  #  check the argument registOutput
  if (!inherits(registOutput, "list")) stop("registOutput is not a list object.")
  
  #  check the argument TrueWarp
  if (!is.null(TrueWarp)) {
    if (!(is.matrix(TrueWarp))) TrueWarp <- as.matrix(TrueWarp)
    if (!(is.matrix(TrueWarp))) stop('TrueWarp needs to be in matrix format')
  }
  
  #  extract the required members of registOutput
  missinginfo <- FALSE
  if (is.null(registOutput$regfd)) {
    missinginfo <- TRUE
  } else {
    yregfd  <- registOutput$regfd
  }
  if (is.null(registOutput$Wfd)) {
    missinginfo <- TRUE
  } else {
    Wfd     <- registOutput$Wfd
  }
  if (is.null(registOutput$yfd)) {
    missinginfo <- TRUE
  } else {
    yfd     <- registOutput$yfd
  }
  if (is.null(registOutput$y0fd)) {
    missinginfo <- TRUE
  } else {
    y0fd    <- registOutput$y0fd
  }
  if (missinginfo) stop("registOutput does not containing required objects.")
  
  #  generate a fine mesh of argument values
  ybasis  <- yfd$basis
  yrange  <- ybasis$rangeval
  ynbasis <- ybasis$nbasis
  if(is.null(Xarg_fine)){
    nfine   <- 201
    argfine <- seq(yrange[1],yrange[2],len = nfine)
  } else{
    nfine <- length(Xarg_fine)
    argfine <- Xarg_fine
  }
  
  if (!is.null(TrueWarp)) argOrig <- seq(yrange[1],yrange[2],len=nrow(TrueWarp))  
  
  #  evaluate the required functional on over this fine mesh
  ymat    <- eval.fd(argfine, yfd)
  y0mat   <- eval.fd(argfine, y0fd)
  yregmat <- eval.fd(argfine, yregfd)
  warpmat <- eval.monfd(argfine, Wfd)
  warpmat <- yrange[1] + (yrange[2]-yrange[1])*warpmat/
    (matrix(1,nfine,1)%*%warpmat[nfine,])
  
  #  extract the number of functions NCURVE and the number of variables NVAR
  ydim     <- dim(yfd$coef)
  ncurve   <- ydim[2]
  
  if (!is.null(names(yfd$fdname)[[3]])) {
    ylabel = names(yfd$fdname)[[3]]
  }  else {
    ylabel = "Function value"
  }
  casename <- yfd$fdname[[2]]
  if (length(casename) != ncurve) casename = as.character(1:ncurve)
  
  if (dim(y0mat)[2] == 1) y0mat = y0mat %*% matrix(1,1,ncurve)
  ylimit   <- c(min(ymat),max(ymat))
  
  Data.Orig <- melt(data = ymat, id = "Curve")
  MainTitle <- paste('Before Registration', '\n', TitleText)
  colnames(Data.Orig) <- c('Pixel', 'Curve', 'Intensity')
  Data.Orig$Pixel <- argfine
  Median_toRegist <- as.data.frame(cbind(Pixel = argfine, Intensity = L1median(t(ymat))$estimate, 
                                         Curve = 'Median'))
  Median_toRegist <- within(data=Median_toRegist,{
    Pixel <- as.numeric(as.vector(Pixel))
    Intensity <- as.numeric(as.vector(Intensity))
  })
  if(is.null(Xlabel))  Xlabel <- 'Pixel'
  if(is.null(Ylabel))  Xlabel <- 'Intensity'
  
  Plot.Orig  <- ggplot(data = Data.Orig, aes_string(x = "Pixel", y = "Intensity", 
                                                    colour = "Curve", group = "Curve")) + 
    geom_line() + 
    ggtitle(MainTitle) + 
    xlab(label = Xlabel) +
    ylab(label = Ylabel) + 
    #    geom_line(aes(x = Pixel, y = Intensity), data = Median_toRegist, size = 2, col = 'white') +
    theme(plot.title = element_text(face = "bold", size = 12, colour = "white"),
          panel.background = element_rect(fill = 'black'), 
          plot.background = element_rect(color = 'black', fill = "gray10"), 
          ## panel.background = element_rect(fill = 'gray60'), 
          ## plot.background = element_rect(color = 'gray60', fill = "gray60"), 
          axis.text = element_text(colour = "white", size = 10), 
          axis.title.x = element_text(colour = "white", size = 12), 
          axis.title.y = element_text(colour = "white", size = 12), 
          panel.grid.major = element_line(colour="gray30", size = 0.35), 
          panel.grid.minor = element_line(colour = "gray20", size = 0.25),
          legend.position = '' 
    )
  
  Data.Regist <- melt(data = yregmat, id = "Curve")
  MainTitle <- paste('After Registration, Using Min Eig Value', '\n', TitleText)
  colnames(Data.Regist) <- c('Pixel', 'Curve', 'Intensity')
  Data.Regist$Pixel <- argfine
  if(is.null(Xlabel))  Xlabel <- 'Pixel'
  if(is.null(Ylabel))  Xlabel <- 'Intensity'
  Median_toRegist <- as.data.frame(cbind(Pixel = argfine, Intensity = L1median(t(yregmat))$estimate, 
                                         Curve = 'Median'))
  Median_toRegist <- within(data=Median_toRegist,{
    Pixel <- as.numeric(as.vector(Pixel))
    Intensity <- as.numeric(as.vector(Intensity))
  })
  
  
  Plot.Regist  <- ggplot(data = Data.Regist, aes_string(x = "Pixel", y = "Intensity", 
                                                        colour = "Curve", group = "Curve")) + 
    geom_line() + 
    ggtitle(MainTitle) + 
    geom_line(aes(x = Pixel, y = Intensity), data = Median_toRegist, size = 2, col = 'white') +
    xlab(label = Xlabel) +
    ylab(label = Ylabel) + 
    theme(plot.title = element_text(face = "bold", size = 12, colour = "white"),
          panel.background = element_rect(fill = 'black'), 
          plot.background = element_rect(color = 'black', fill = "gray10"), 
          ## panel.background = element_rect(fill = 'gray60'), 
          ## plot.background = element_rect(color = 'gray60', fill = "gray60"), 
          axis.text = element_text(colour = "white", size = 10), 
          axis.title.x = element_text(colour = "white", size = 12), 
          axis.title.y = element_text(colour = "white", size = 12), 
          panel.grid.major = element_line(colour = "gray30", size = 0.35), 
          panel.grid.minor = element_line(colour = "gray20", size = 0.25),
          legend.position = '' 
    )
  
  Consensus_Mean <- rowMeans(yregmat)
  Consensus_SE <- rowSE( Matrix = yregmat )
  Consensus_MeanSE <- as.data.frame(cbind(
    Intensity = Consensus_Mean, 
    LL = Consensus_Mean - 2*Consensus_SE,
    UU = Consensus_Mean + 2*Consensus_SE, 
    Pixel = argfine))
  Consensus_MeanSE$Curve <- 'Mean'
  
  Plot.Regist_wSE <- Plot.Regist +
    geom_smooth( aes(y = Intensity, x = Pixel, ymin = LL, ymax = UU), data = Consensus_MeanSE, stat = 'identity', 
                 color = 'gray60', size = 1, fill = 'white', alpha = 0.9)
  
  Data.Warp <- melt(data = warpmat, id = "Curve")
  MainTitle <- paste('Warping functions, Using Min Eig Value', '\n', TitleText)
  colnames(Data.Warp) <- c('Pixel', 'Curve', 'Warp')
  yrangediff <- yrange[2] - yrange[1] + 1
  Data.Warp$Pixel <- Data.Warp$Pixel * yrangediff/ nfine + yrange[1] - 1
  Data.Warp$Curve <- as.factor(Data.Warp$Curve)
  
  Plot.Warp  <- ggplot(data = Data.Warp, aes_string(x = "Pixel", y = "Warp", 
                                                    colour = "Curve", group = "Curve")) + 
    geom_line() + 
    ggtitle(MainTitle) +
    geom_abline(intercept = 0, slope = 1, col = 'white', size = 1) + 
    xlab(label = Xlabel) +
    theme(plot.title = element_text(face = "bold", size = 12, colour = "white"),
          panel.background = element_rect(fill = 'black'), 
          plot.background = element_rect(color = 'black', fill = "gray10"), 
          ## panel.background = element_rect(fill = 'gray60'), 
          ## plot.background = element_rect(color = 'gray60', fill = "gray60"), 
          axis.text = element_text(colour = "white", size = 10), 
          axis.title.x = element_text(colour = "white", size = 12), 
          axis.title.y = element_text(colour = "white", size = 12), 
          panel.grid.major = element_line(colour="gray30", size = 0.35), 
          panel.grid.minor = element_line(colour="gray20", size = 0.25),
          legend.position = '' 
    )
  
  ## ymat: original; 
  ## y0mat: mean to be registered to
  ## yregmat: registered curve
  if( BeforeAfterDistance ){
    PhaseDist.Before <- fn_pairwiseDistance_fdasrvf(Mat = ymat, Xaxis = argfine)
    PhaseDist.After <- fn_pairwiseDistance_fdasrvf(Mat = yregmat, Xaxis = argfine)
    PhaseDistPlot.Before <- qplot(Dx, data = PhaseDist.Before, geom = "histogram", binwidth = 0.02) + 
      ggtitle('Pairwise elastic distance, Before registration') +
      xlab('Phase Distance')
    PhaseDistPlot.After <- qplot(Dx, data = PhaseDist.After, geom = "histogram", binwidth = 0.02) + 
      ggtitle('Pairwise elastic distance, After registration') +
      xlab('Phase Distance')
    colnames(Data.Warp) <- c('OriginalTime', 'Curve', 'WarpedTime')
  }
  
  if(saveToPDF){
    if(!is.null(Xlabel)){
      Iteration <- as.numeric(substr(x=Xlabel, start=nchar('After Iteration  '), stop=nchar(Xlabel)))
      if(PlotBeforeRegist | Iteration == 1) plot(Plot.Orig)
    } else{
      if( PlotBeforeRegist) plot(Plot.Orig)
    }
    plot(Plot.Regist)
    plot(Plot.Regist_wSE)
    plot(Plot.Warp)
    if( BeforeAfterDistance ) grid.arrange(PhaseDistPlot.Before, PhaseDistPlot.After, ncol = 1)
    return(Data.Warp)
  } else{
    AllPlots <- vector(mode='list', length=5) 
    if(!is.null(Xlabel)){
      Iteration <- as.numeric(substr(x=Xlabel, start=nchar('After Iteration  '), stop=nchar(Xlabel)))
      if(PlotBeforeRegist | Iteration == 1) AllPlots[[1]] <- Plot.Orig
    } else{
      if( PlotBeforeRegist) AllPlots[[1]] <- Plot.Orig
    }
    AllPlots[[2]] <- Plot.Regist
    AllPlots[[3]] <- Plot.Regist_wSE
    AllPlots[[4]] <- Plot.Warp
    if( BeforeAfterDistance ) AllPlots[[5]] <- grid.arrange(PhaseDistPlot.Before, PhaseDistPlot.After, ncol = 1)
    names(AllPlots) <- c('Plot.Orig', 'Plot.Regist', 'Plot.Regist_wSE', 'Plot.Warp', 'Plot.PhaseDist')
    return(list(AllPlots = AllPlots, Data.Warp = Data.Warp))
  }
}