R/fmri_3dvisual.R

In TCIU: Spacekime Analytics, Time Complexity and Inferential Uncertainty

#' @title visualization of the 3D brain with the activated areas
#' @description a visualization method, using \code{plotly} to draw the 3D plot of the brain 
#' with the activated areas determined by p-values, which is generated from fMRI data
#'
#' @param pval a 3D array of p-values used to plot activated area of the brain
#' @param mask a 3D nifti or 3D array of data to show the shell of the brain
#' @param p_threshold NULL or a numeric value that can be selected randomly below 0.05 to 
#' drop all p-values above the threshold. If 'low5_percent' method is used, 
#' make 'p_threshold' as NULL. The default is 0.05.
#' @param method a string that represents method for the plot.
#' There are 2 options: 'scale_p' and 'low5_percent'. The default is 'scale_p'. 
#' 'scale_p' is to draw the plot with fixed color scale for fixed range of p value.
#' 'low5_percent' is to draw the plot for the smallest 5 percent of p value when all the p values are not significant.
#' @param color_pal the name of the color palettes provided by \code{RColorBrewer}. The default is "YlOrRd".
#' @param multi_pranges an option under 'scale_p' method to decide whether there are at most 9 colors 
#' in the legend for the ranges of p value, or at most 4 colors. 
#' The default is TRUE, choosing the larger number of colors for the plot.
#' @param title the title of the plot. The default is NULL.
#' 
#' @details
#' The function \code{fmri_3dvisual} is used to visualize the 3D plot of the brain 
#' with its activated parts based on provided p values. The p values are generated by
#' applying statistical test on fMRI data. When providing input of a 3D p-values data, 
#' a 3D interactive plot will be generated with surface for the brain shell 
#' and scatter points in different colors and size representing different stimulated levels.
#'
#' @author SOCR team <\url{http://socr.umich.edu/people/}>
#'
#' @return a list of two elements
#' \itemize{
#'   \item plot - the 3d plot of the fMRI data drawn by \code{plotly}
#'   \item pval_df - data.frame with the p value for each voxel and the specified color for it
#' }
#'
#' @examples
#' # sample 3D data of mask provided by the package
#' dim(mask)
#' # sample 3D p value provided by the package
#' dim(phase2_pval)
#' 
#' # make the 3D plot
#' fmri_3dvisual(phase2_pval, mask, p_threshold = 0.05, method="scale_p")$plot
#' 
#' @export
#'
#' @import plotly scales RColorBrewer fancycut geometry dplyr
#' @importFrom grDevices colorRampPalette contourLines

fmri_3dvisual = function(pval,
                         mask,
                         p_threshold = 0.05, 
                         method = "scale_p",
                         color_pal = "YlOrRd",
                         multi_pranges = TRUE,
                         title = NULL){

  floor_dec = function(x, level=1) round(x - 5*10^(-level-1), level)
  ceiling_dec = function(x, level=1) round(x + 5*10^(-level-1), level)
  # "boundry_pt_func()" is a small function, specialized to generate boundary point for the 2d z-plane 
  # of a 3d structure based on the the location of z
  # the input of this function are the location of z and the 3d strucure as a 3d array 
  # the output is a dataframe which contains the index of the boundary points for 
  # 2d plane of the 3d structure at location z
  boundry_pt_func = function(z, struc_3d=mask){
    dim_maskz = dim(struc_3d[,,z])
    x = seq(1, dim_maskz[1])
    y = seq(1, dim_maskz[2])
    contour_pt = contourLines(x = x, y = y, z = struc_3d[,,z])
    if (length(contour_pt) != 0){
      boundry_pt_df = cbind(contour_pt[[1]]$x, contour_pt[[1]]$y, z)
      return(boundry_pt_df)
    }
  }
  
  # generate cutomized error message based on different wrong input
  if((is.array(pval)!=TRUE) |
      (length(dim(pval))!=3)){
    stop("'pval' should be a 3d array.")
  }else if(((class(mask) %in% c("array", "nifti")) != TRUE) |
            (length(dim(mask))!=3)){
    stop("'mask' should be a 3d array.")
  }else if(is.null(p_threshold) != TRUE){
    if((is.numeric(p_threshold) != TRUE) | 
        (p_threshold > 0.05) | (p_threshold <= 0)){
      stop("'p_threshold should be a numeric value in range of (0, 0.05].'")
    }
  }else if((method %in% c("scale_p", "low5_percent")) != TRUE){
    stop("'method' should only choose from 'scale_p' or 'low5_percent'.")
  }
  
  # this part is for the mesh plot of the surface of the 3d brain
  # construct a dataframe for the index of all the boundary points for the 3d brain
  boundry_pt_df = 
    data.frame(do.call(rbind, lapply(1:dim(mask)[3], boundry_pt_func)))
  colnames(boundry_pt_df) = c("x", "y", "z")
  
  # get the delaunay triangles for the surface of the 3d brain to have the mesh plot
  dela_contour = delaunayn(boundry_pt_df)
  zmean = apply(dela_contour, MARGIN=1, function(row){mean(boundry_pt_df[row, 3])})
  # assign color to the mesh plot of the surface of the 3d brain
  color_sel = colour_ramp(colorRampPalette(c("#E0E0E0", "#D9D9D9"))(10))(scales::rescale(x=zmean))
  
  p = plot_ly()%>%
    add_trace(boundry_pt_df, 
              x=boundry_pt_df[, 1], y=boundry_pt_df[, 2], z=boundry_pt_df[, 3],
              i = dela_contour[, 1]-1, j = dela_contour[, 2]-1, k = dela_contour[, 3]-1,
              type = "mesh3d", opacity = 0.01, name = "Brain Shell",
              # facecolor = color_sel,
              contour = list(show = TRUE, color="#000", width=15))
  if (!is.null(title)){
  p = p%>%
      layout(title = title) 
  }
  
  # this part is to process p value data to add the points of the motor area in the mesh brain plot
  # this plot can only be made on 1 - p value
  one_minus_p_3d = (1- pval)
  dim_pval = dim(one_minus_p_3d)
  # convert 3d array 1 - p value to a dataframe to be the input for plotly
  pval_df = 
    data.frame(x = rep(c(1:dim_pval[1]), dim_pval[2]*dim_pval[3]),
                y = rep(rep(c(1:dim_pval[2]), each=dim_pval[1]), dim_pval[3]),
                z = rep(c(1:dim_pval[3]), each=dim_pval[1]*dim_pval[2]))
  names(pval_df) = c("x", "y", "z")
  pval_df$one_minus_p_3dn = as.vector(one_minus_p_3d)
  pval_df$p_val = 1 - pval_df$one_minus_p_3dn
  
  
  # pval_df_filter = filter(pval_df, p_val <= p_threshold)
  # dim_pval_df_filter = dim(pval_df_filter)
  # if(p_threshold==0.05 & dim_pval_df_filter == 0){
  #   stop("There is no p value survived below 0.05. Change the method to 'low5_percent'!")
  # }else if(dim_pval_df_filter == 0){
  #   stop("There is no p value survived below the selected p_threshold. 
  #        Change one p_threshold or the plot method!
  #        If no p value survive under 0.05, try 'low5_percent'!")
  # }
  # add message to tell if all data >0.05 or p_threshold > 0.05, suggest low5percent
  switch (method,
          "scale_p" = {
            
            pval_df = dplyr::filter(pval_df, pval_df$p_val <= p_threshold)
            pval_df$cut_invs = wafflecut(pval_df$p_val, 
                                          c('[0,1e-8]', '(1e-8, 1e-7]', '(1e-7, 1e-6]', '(1e-6, 1e-5]', 
                                            '(1e-5, 1e-4]', '(1e-4, 1e-3]', '(1e-3, 1e-2]', '(1e-2, 5e-2]'))
            # consider whether we can decide to change the interval cut based on the condition of p calue provided
            pval_df$colorgrp = as.numeric(pval_df$cut_invs)
            if (length(unique(pval_df$colorgrp)) <= 4){
              message("The number of p ranges is originally less than 4, 
                    it is recommended to turn 'multi_pranges' as FALSE.")
            }
            if(multi_pranges == FALSE){
              pval_df$cut_invs = wafflecut(pval_df$p_val, 
                                            c('[0,1e-7]', '(1e-7, 1e-5]', 
                                              '(1e-5, 1e-3]', '(1e-3, 5e-2]'))
              pval_df$colorgrp = as.numeric(pval_df$cut_invs)
            }
            
            label_p = levels(unique(pval_df$cut_invs))
            
            },
          
          "low5_percent" = {
            
            quantile5_pt = quantile(pval_df$p_val, 0.05)
            pval_df = dplyr::filter(pval_df, pval_df$p_val <= quantile5_pt)
            p_val_max_minus_min = max(pval_df$p_val)-min(pval_df$p_val)
            if (p_val_max_minus_min != 0){
              
              cut_pts = min(pval_df$p_val) + 0:9 * (p_val_max_minus_min) / 9
              cut_pts[length(cut_pts)] = ceiling_dec(cut_pts[length(cut_pts)], 2)
              cut_pts[1] = floor_dec(cut_pts[1], 2)
              # add lapply to cut pts[2:-1] to round 3
              cut_pts = unique(sapply(cut_pts, function(x) round(x, 3)))
              invs_total = levels(cut(pval_df$p_val, cut_pts, right=TRUE))
              invs_total[1] = sub("\\(", "[", invs_total[1])
              pval_df$cut_invs = wafflecut(pval_df$p_val, invs_total)
              pval_df$colorgrp = as.numeric(pval_df$cut_invs)
              label_p = levels(unique(pval_df$cut_invs))
            }else{
              
              pval_df$cut_invs = as.character(round(min(pval_df$p_val), 3))
              label_p = unique(pval_df$cut_invs)
              pval_df$colorgrp = 1
            }
          }
  )

  if (multi_pranges == FALSE){
    color_choice = rev(brewer.pal(9, color_pal))[c(FALSE, TRUE)]
  }else{
    color_choice = rev(brewer.pal(9, color_pal))
  }
  pval_df$corresp_color = color_choice[pval_df$colorgrp]

  for (i in sort(unique(pval_df$colorgrp))){
    pts_grp = pval_df[pval_df$colorgrp == i, ]
    p = add_markers(p, data=pts_grp, 
                    x=~x, y=~y, z=~z, mode="markers", 
                    name=paste0("p value in ", label_p)[i],
                    marker = list(opacity=0.6, symbol="circle", 
                                  size=rev(c(1:length(color_choice)))[i]+2*as.numeric(I(multi_pranges==FALSE)),
                                  color=~corresp_color,
                                  line = list(color = ~corresp_color,
                                              width = 0.3))
    )
  }
  
  return(list(plot=p, pval_df = pval_df))
}