R/plot_regist_fdasrvf.R

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

# registOutput <- Out_Male
# TrueWarp=NULL
# Lambda=NULL
# Ylabel=NULL
# Xlabel=NULL
# time=Age[-1]

#' @import fda
#' @importFrom gridExtra grid.arrange
#' @import ggplot2
#@import fdasrvf
#' @import reshape2

plotbyCurve_regist_fdasrv <- function(registOutput, TrueWarp=NULL, 
                                      time=Age[-1], 
                                      Lambda=NULL,
                                      Ylabel=NULL,
                                      Xlabel=NULL
) {
  # plots the output of registration function time_warping, one curve at a time
  #
  # Argument:
  # registOutput ... The named list generated by function time_warping
  #     f0    original functions
  #     fn    aligned functions - matrix (N x M ) of M functions with N samples
  #     qn    aligned srvfs - similar structure to fn
  #     q0    original srvf - similar structure to fn
  #     fmean function mean or median - vector of length N
  #     mqn   srvf mean or median - vector of length N
  #     gam   warping functions - similar structure to fn
  
  # 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$f0)) {
    missinginfo <- TRUE
  } else {
    f0  <- registOutput$f0
  }
  if (is.null(registOutput$fn)) {
    missinginfo <- TRUE
  } else {
    fn  <- registOutput$fn
  }
  if (is.null(registOutput$gam)) {
    missinginfo <- TRUE
  } else {
    gam  <- registOutput$gam
  }
  
  if (missinginfo) stop("registOutput does not containing required objects.")
  
  #  generate a fine mesh of argument values
  yrange  <- range(fn)
  
  if (!is.null(TrueWarp)) argOrig <- seq(yrange[1],yrange[2],len=nrow(TrueWarp))  
  
  #  evaluate the required functional on over this fine mesh
  ymat   <- f0
  yregmat <- fn
  warpmat <- gam
  
  #   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
  ncurve   <- ncol(fn)
  
  ylimit   <- c(min(yregmat),max(yregmat))
  
  i <- 1
  for (i in 1:ncurve) {
    ## ymat: original; 
    ## y0mat: mean to be registered to
    ## yregmat: registered curve
    Data <- as.data.frame(cbind(time=time, ymat=ymat[,i], y0mat=rowMeans(f0), yregmat=yregmat[,i]))
    if(is.null(Xlabel)) Xlabel <- 'Pixel Position'
    if(is.null(Ylabel)) Ylabel <- 'Intensity'
    MainTitle <- paste("Curve Number", i)
    Col.Orig <- 'deepskyblue'
    Col.Mean <- 'orange'
    Col.Reg <- 'springgreen'
    
    PlotCurve <- ggplot()+
      ylim(ylimit) +
      geom_line(data=Data, aes(x=time, y=ymat, col='Original'), linetype=2, size=1) +
      geom_line(data=Data, aes(x=time, y=y0mat, col='Mean'), linetype=1, size=1) + 
      geom_line(data=Data, aes(x=time, 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"), 
            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')
      )
    
    #PlotCurve
    
    Data.Warp <- as.data.frame(cbind(time=time, Warp=warpmat[,i]))
    Data.Warp$Warp <- Data.Warp$Warp * 18
    
    if (is.null(TrueWarp)){
      PlotWarp <- ggplot()+
        geom_line(data=Data.Warp, aes(x=time, 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, SRVF') +
        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)  
  }
}

plotAll_regist_fdasrv <- function(registOutput, TrueWarp=NULL, time,
                                  Lambda=NULL,
                                  Ylabel=NULL,
                                  Xlabel=NULL
) {
#   require(ggplot2)
#   require(gridExtra)
#   require(reshape2)
  # plots the output of registration function time_warping, one curve at a time
  #
  # Argument:
  # registOutput ... The named list generated by function time_warping
  #     f0    original functions
  #     fn    aligned functions - matrix (N x M ) of M functions with N samples
  #     qn    aligned srvfs - similar structure to fn
  #     q0    original srvf - similar structure to fn
  #     fmean function mean or median - vector of length N
  #     mqn   srvf mean or median - vector of length N
  #     gam   warping functions - similar structure to fn
  
  # 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$f0)) {
    missinginfo <- TRUE
  } else {
    f0  <- registOutput$f0
  }
  if (is.null(registOutput$fn)) {
    missinginfo <- TRUE
  } else {
    fn  <- registOutput$fn
  }
  if (is.null(registOutput$gam)) {
    missinginfo <- TRUE
  } else {
    gam  <- registOutput$gam
  }
  
  if (missinginfo) stop("registOutput does not containing required objects.")
  
  #  generate a fine mesh of argument values
  yrange  <- range(fn)
  
  if (!is.null(TrueWarp)) argOrig <- seq(yrange[1],yrange[2],len=nrow(TrueWarp))  
  
  #  evaluate the required functional on over this fine mesh
  ymat   <- f0
  yregmat <- fn
  warpmat <- gam
  
  #   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
  ncurve   <- ncol(fn)
  
  ylimit   <- c(min(yregmat),max(yregmat))
  
  Data.Orig <- melt(data=ymat, id="Curve")
  MainTitle <- paste('Original Curves, Using SRVF')
  colnames(Data.Orig) <- c('Pixel', 'Curve', 'Intensity')
  Data.Orig$Curve <- as.factor(Data.Orig$Curve)
  Plot.Orig  <- ggplot(data=Data.Orig, aes_string(x="Pixel", y="Intensity", 
                                                  colour="Curve", group="Curve")) + 
    geom_line() + 
    ggtitle(MainTitle) +
    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), 
          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 = '' 
    )
  plot(Plot.Orig)
  
  Data.Regist <- melt(data=yregmat, id="Curve")
  MainTitle <- paste('Registered Curves, Using SRVF')
  colnames(Data.Regist) <- c('Pixel', 'Curve', 'Intensity')
  Data.Regist$Curve <- as.factor(Data.Regist$Curve)
  Plot.Regist  <- ggplot(data=Data.Regist, aes_string(x="Pixel", y="Intensity", 
                                                      colour="Curve", group="Curve")) + 
    geom_line() + 
    ggtitle(MainTitle) +
    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), 
          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 = '' 
    )
  plot(Plot.Regist)
  
  Data.Warp <- melt(data=warpmat, id="Curve")
  MainTitle <- paste('Warping functions, Using SRVF')
  colnames(Data.Warp) <- c('Pixel', 'Curve', 'Warp')
  Data.Warp$Pixel <- (Data.Warp$Pixel - min(Data.Warp$Pixel))/
    (max(Data.Warp$Pixel) - min(Data.Warp$Pixel))
  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) + 
    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), 
          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 = '' 
    )
  plot(Plot.Warp)
  
  PhaseDist.Before <- fn_pairwiseDistance_fdasrvf(Mat=f0, Xaxis=time)
  PhaseDist.After <- fn_pairwiseDistance_fdasrvf(Mat=fn, Xaxis=time)
  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')
  plot(grid.arrange(PhaseDistPlot.Before, PhaseDistPlot.After, ncol=1))
  return(Data.Warp)
}