R/plotBasicNetwork.R

In neuroconductor-devel/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Simple plotBasicNetwork function.
#'
#' takes an object output from renderSurfaceFunction and a list of centroids
#' and plots the centroid network over the rendering object
#'
#' If \code{edgecolors} is not specified, a heat-like color palette is used. Weights
#' can be quantile transformed or clipped at a given quantile in order to improve
#' contrast.
#'
#' If weights are not specified, only the nodes are plotted.
#'
#' @param centroids input matrix of size N 3D points ( in rows ) by 3 (
#' in columns ), for N nodes.
#' @param brain input rendering object which is output of renderSurfaceFunction.
#' or a function derived from renderSurfaceFunction.
#' @param weights edge weights, a symmetric matrix of size N. Weights should be non-negative.
#' @param nodecolors a color or color vector for nodes.
#' @param edgecolors a color(map) for edges. If a color map function, weights will be transformed
#' to the range [0,1], which is compatible with functions returned by \code{colorRamp}.
#' @param backgroundColor background color.
#' @param nodetype sphere or other node type supported by RGL.
#' @param edgeContrast a vector of length 2, specifying the contrast range for edge colors.
#' Weights are normalized to the range [0,1].
#' The normalized weights can be rescaled with this parameter, eg \code{c(0.05,0.95)}
#' would stretch the contrast over the middle 90\% of normalized weight values.
#' @param quantileTransformWeights quantile transform the weights.
#' @param lwd line width for drawing edges.
#' @param minRadius minimum node radius. Ignored if the radius is specified explicitly with \code{radius}.
#' @param maxRadius maximum node radius. The node radius between \code{minRadius} and \code{maxRadius} is
#' determined from the sum of the edge weights connecting the node. Ignored if the radius is specified
#' explicitly with \code{radius}.
#' @param radius a constant radius or vector of length nrow(centroids). If not specified, node radius is
#' determined by the sum of edge weights connected to the node.
#' @param showOnlyConnectedNodes boolean, if \code{TRUE}, only nodes with non-zero edge weights are plotted.
#' 
#' @return None
#' @author Avants BB, Duda JT, Cook PA
#' @examples
#'
#' \dontrun{ 
#' # more complete example
#'   mnit<-getANTsRData("mni")
#'   mnit<-antsImageRead(mnit)
#'   mnia<-getANTsRData("mnia")
#'   mnia<-antsImageRead(mnia)
#'   mnit<-thresholdImage( mnit, 1, max(mnit) )
#'   mnit<-iMath(mnit,"FillHoles") 
#'   mnit = thresholdImage(mnit,  1, 2 )
#'   cnt<-getCentroids( mnia, clustparam = 0 )
#'   aalcnt<-cnt[1:90,1:3]
#'   brain<-renderSurfaceFunction( surfimg =list( mnit ) , 
#'   alphasurf=0.1 ,smoothsval = 1.5 )
#'   testweights<-matrix( rep( 0, 90*90 ) ,nrow=90)
#'   testweights[31,37]<-1  # ant cingulate to hipp
#'   testweights[31,36]<-2  # ant cingulate to post cingulate
#'   testweights[11,65]<-3  # broca to angular
#'   plotBasicNetwork( centroids = aalcnt , brain , weights=testweights, edgecolors = "red" )
#'   id<-rgl::par3d('userMatrix')
#'   rid<-rotate3d( id , -pi/2, 1, 0, 0 )
#'   rid2<-rotate3d( id , pi/2, 0, 0, 1 )
#'   rid3<-rotate3d( id , -pi/2, 0, 0, 1 )
#'   rgl::par3d(userMatrix = id )
#'   dd<-make3ViewPNG(  rid, id, rid2,  paste('network1',sep='') )
#'   rgl::par3d(userMatrix = id )
#' # another example 
#' mni<-getANTsRData("mni")
#' mni<-antsImageRead(mni)
#' mnit<-thresholdImage( mni, 1, max(mni) )
#' mnit<-iMath(mnit,"FillHoles")
#' mniseg = kmeansSegmentation( mni, 3 )$segmentation
#' wmbkgd = thresholdImage( mniseg, 3, 3 ) %>% 
#' iMath( "GetLargestComponent" ) %>% 
#' iMath( "FillHoles" )
#' wmbkgd = smoothImage( iMath( wmbkgd, "MD", 1 ), 2.0 )
#' brain<-renderSurfaceFunction( surfimg =list( wmbkgd ) , 
#' alphasurf=0.8 ,smoothsval = 1.0 )
#' data( powers_areal_mni_itk )
#' coords = powers_areal_mni_itk[,1:3]
#' id<-rgl::par3d('userMatrix')
#' rid<-rotate3d( id , -pi/2, 1, 0, 0 )
#' rid2<-rotate3d( id , pi/2, 0, 0, 1 )
#' rid3<-rotate3d( id , -pi/2, 0, 0, 1 )
#' rgl::par3d(userMatrix = id )
#' handMat2 = t( matrix(  c(-0.9998656511 , 0.01626961,  0.00198165 ,   
#' 0 ,-0.0163816363, -0.99584705 ,-0.08955579   , 0, 0.0005163439, 
#' -0.08957647,  0.99597979 ,   0,  0.0000000000,  0.00000000,  
#' 0.00000000  ,  1),  ncol = 4) )
#' loccolor = as.character( powers_areal_mni_itk$Color )
#' loccolor[ loccolor == "Peach" ] = "sienna1"
#' loccolor[ loccolor == "Cyan" ] = "cyan"
#' loccolor[ loccolor == "Orange" ] = "orange"
#' loccolor[ loccolor == "Purple" ] = "darkorchid1"
#' loccolor[ loccolor == "Pink" ] = "deeppink"
#' loccolor[ loccolor == "Red" ] = "red"
#' loccolor[ loccolor == "Gray" ] = "gray74"
#' loccolor[ loccolor == "Teal" ] = "turquoise4"
#' loccolor[ loccolor == "Blue" ] = "blue"
#' loccolor[ loccolor == "Yellow" ] = "yellow"
#' loccolor[ loccolor == "Black" ] = "black"
#' loccolor[ loccolor == "Brown" ] = "brown"
#' loccolor[ loccolor == "Pale blue" ] = "steelblue1"
#' loccolor[ loccolor == "Green" ] = "green"
#' tt = plotBasicNetwork( centroids = coords, brain, 
#' nodecolors = loccolor, radius=3 )
#' dd<-make3ViewPNG(  handMat2, id, rid2, tempfile( fileext='.png' ) )
#' rgl::par3d(userMatrix = id )
#' }
#'
#' @export plotBasicNetwork
plotBasicNetwork <- function(
                             centroids,
                             brain,
                             weights = NA,
                             edgecolors = -1,
                             backgroundColor = "white",
                             nodecolors = "blue",
                             nodetype = "s",
                             edgeContrast = c(0, 1),
                             quantileTransformWeights = FALSE,
                             lwd = 2,
                             minRadius=0,
                             maxRadius=3,
                             radius = NA,
                             showOnlyConnectedNodes = TRUE ) {
    if (missing(centroids) | missing(brain)) {
        print(args(plotBasicNetwork))
        return(1)
    }
    if(!usePkg('rgl')){
        print("rgl is necessary for this function.")
        return(NULL)
    }

    edgeColorsIsFunction <- is.function(edgecolors)

    # Test if weights are all NA here. Then use this boolean in place of (is.na(weights), which
    # causes warnings in scalar contexts like if statements
    weightsNA <- all(is.na(weights))

    numNodes <- nrow(centroids)
    
    rgl::rgl.bg(color = backgroundColor)
    rgl::par3d(windowRect = c(100, 100, 600, 600))
    mesh <- .getvertices(brain[[1]])
    nSurfaceVerts <- dim(mesh$vertices)[1]
    mesh$vertices <- rbind(mesh$vertices, as.matrix(centroids))
    labelVerts <- c(1:nrow(centroids)) + nSurfaceVerts
    
    if (weightsNA) {
        if ( all(is.na( radius )) ) radius = 3
        rgl::spheres3d(mesh$vertices[labelVerts, ], color = nodecolors, type = nodetype, radius = radius)
        return(1)
    }
  
    if ( !(length(dim(weights) == 2 && dim(weights)[1] == numNodes && dim(weights)[2] == numNodes)) ) {
        stop("Weights must be a 2D symmetric matrix with one entry for each pair of nodes")
    }
    
    if ( all(is.na(radius)) ) { # node scaled by strength
        radiusw <- rep(0, nrow(centroids))
        radiusscale <- as.numeric(apply(weights, FUN = sum, MARGIN = 1, na.rm = T) +
                                  apply(weights, FUN = sum, MARGIN = 2, na.rm = T))
        radiusscale <- ( radiusscale/max(radiusscale) )
        nodeMask <- 1 * (radiusscale > 0) # unconnected nodes stay at 0 if showOnlyConnectedNodes
        radiusw <- (minRadius + (maxRadius - minRadius) * radiusscale)

        if (showOnlyConnectedNodes) {
            radiusq <- radiusw * nodeMask
        }
      
        radius <- radiusw
    }
   
    ggg <- weights
    luniqvals <- length( unique( ggg[ggg > 0] ) )
    
    binaryWeights <- FALSE
    
    if (luniqvals == 1) {
        binaryWeights <- TRUE
        # If some constant that is not binary, convert to binary so that we don't mess up
        # color mapping
        weights <- weights / max(weights)
    }
    else if (quantileTransformWeights) {
        if ( luniqvals > 250 ) ncuts = 1.0 / 250.0 else ncuts = 1.0 / ( luniqvals*0.5 )
        
        qqq = quantile(ggg[ ggg > 0 ], probs=seq(0, 1, by=ncuts), na.rm=TRUE, names=FALSE, type=7 )
        qqq = unique( qqq )
        myquartile <- cut(ggg, breaks = qqq, include.lowest=TRUE )
        myquartile <- as.numeric( myquartile )
        myquartile[ is.na( myquartile) ] = 0
        myquartile[ is.nan( myquartile) ] = 0
        weights <- matrix( myquartile, nrow=nrow(weights) )
    }
    
    rgl::spheres3d(mesh$vertices[labelVerts, ], color = nodecolors, type = nodetype, radius = radius)

    edgelocations <- c()
    edgeweights <- c()
    
    
    for (i in c(1:nrow(weights))) {
        for (j in c(1:ncol(weights))) {
            if (weights[i, j] > 0 & weights[i, j] < Inf) {
                # Draw each edge once only
                if (i < j) { 
                    edgelocations <- c(edgelocations, nSurfaceVerts + c(i, j))
                    edgeweights <- c(edgeweights, weights[i, j])
                }
            }
        }
    }

    # normalize weights to range [0-1]
    edgeweightsNorm <- edgeweights
    
    if (!binaryWeights) {
        minWeight <- min(edgeweights)
        maxWeight <- max(edgeweights)
        
        edgeweightsNorm <- (edgeweights - minWeight) / (maxWeight - minWeight)
        
        # user defined contrast stretch
        # Want to map (0.05, 0.95) to (0,1)
        # then clip values outside this range
        edgeweightsNorm <- (edgeweightsNorm - edgeContrast[1]) / (edgeContrast[2] - edgeContrast[1])
        
        edgeweightsNorm[edgeweightsNorm > 1] <- 1
        edgeweightsNorm[edgeweightsNorm < 0] <- 0
    }
    
    # Map colors unless color map or other explicit color scheme provided.
    if (!edgeColorsIsFunction && (length(edgecolors) == 1) && (edgecolors[1] == -1)) {
        # heat.colors makes high end white, which is invisible with default background
        
        colormap <- colorRampPalette(c("#650000","dark red", "red", "darkorange", "orange", "yellow", "#FAFAD0"))
        
        edgecolors <- edgeweightsNorm
        
        colorArr <- colormap(256)
        
        for (i in c(1:length(edgeweightsNorm))) {
            edgecolors[i] <- colorArr[1 + floor(edgeweightsNorm[i] * 255)]
        }
        
    }
    if (edgeColorsIsFunction) {
        # Convert function to color array
        colSeq <- seq(0,1,1/255)
        
        colorArr <- rgb(edgecolors(colSeq), maxColorValue = 255)

        edgecolors <- edgeweightsNorm
        
        for (i in c(1:length(edgeweightsNorm))) {
            edgecolors[i] <- colorArr[1 + floor(edgeweightsNorm[i] * 255)]
        }
        
    }
    
    rgl::segments3d(mesh$vertices[edgelocations, ], col = edgecolors, lwd = lwd)
    
}


.getvertices <- function(inrglmesh) {
    cter <- nrow(inrglmesh[[1]])
    vertices <- matrix(NA, nrow = 3 * cter, ncol = 3)
    inds <- c(1:cter)
    vertices[(3 * inds - 2), ] <- inrglmesh[[1]]
    vertices[(3 * inds - 1), ] <- inrglmesh[[2]]
    vertices[(3 * inds - 0), ] <- inrglmesh[[3]]
    indices <- rep(NA, nrow(vertices))
    return(list(vertices = vertices, indices = indices))
}