R/midplane.reflect.R

In DJFernandes/ABIgeneRMINC: Read, Write, and Analyzing Allen Gene Expression Raw Data

#' Reflect Gene Expression Vector across midplane
#'
#' @param gene.expression gene expression data to reflect (as either vector or array)
#' @param reflect.dimension INTEGER. Which dimension to reflect. This is done automatically if you use vectors and arrays generated by this package. If you want to do it manually, put 1 for the slowest dimension, 3 for the fastest. To reflect Allen orientation vectors across the sagittal midplane, set this value to 1. To reflect MINC orientation vectors across the sagittal midplane, set this value to 3. 
#' @param array.dimensions Vector of length 3 indicating the dimension sizes in order of the SLOWEST dimension to the FASTEST. Must be supplied if gene.expression.vector is missing 'sizes' attributes (execute attr(gene.expression.vector,'sizes')).
#' @param data If "gene" (default), then assume missing data has value -1, which is the syntax chosen by ABI for their gene expression data. If "atlas", then missing data is 0, which is the syntax chosen by ABI`s for their atlases.
#'
#' @return Gene Expression as vector or array
#'
#' @examples
#' # Read Bdnf gene expression (Experiment 75695642) from Allen Brain Institute
#' # Look at a coronal slice
#' # Reflect across the sagittal midplane 
#' # Look at a coronal slice
#'
#'
#' ## Read Bdnf expression (Experiment 75695642)
#' URLtoDownload='http://api.brain-map.org/grid_data/download/75695642?include=energy'
#' gene.expression.vector=read.raw.gene(URLtoDownload,url=TRUE)
#' 
#' ## Coronal slice
#' c.slice=array(gene.expression.vector,rev(attr(gene.expression.vector,'sizes')))[30,,]
#' image(c.slice, ylab='Left-Right' ,xlab='Superior-Inferior')
#' 
#' ## Reflect across the sagittal midplane 
#' reflected.vector=midplane.reflect(gene.expression.vector,1)
#' 
#' ## Coronal slice
#' reflected.c.slice=array(reflected.vector,rev(attr(gene.expression.vector,'sizes')))[30,,]
#' image(reflected.c.slice, ylab='Left-Right' ,xlab='Superior-Inferior')
#'
#' ## Reflect MINC orientation vector across sagittal midplane
#' minc.vector=allenVectorTOmincVector(gene.expression.vector)
#' reflected.minc.vector=midplane.reflect(minc.vector,3)
#' @export

midplane.reflect=function(gene.expression,reflect.dim=NULL,minc.array.sizes=NULL,data='gene'){
     if (!data %in% c('gene','atlas')) {
        stop('data must be either gene or atlas')
     }
     if (is.null(dim(gene.expression))) {
        d1 = T
        gene.expression.vector = gene.expression
     } else {d1 = F}

     if (d1) {
        if (is.null(attr(gene.expression.vector,'sizes'))) {
        if (is.null(array.dimensions) | length(array.dimensions) != 3) { 
           stop(paste("no attribute sizes in gene.expression.vector. Must provide attribute OR",
           "array.dimensions vector of length 3"))
        } else {arrdim=rev(minc.array.sizes)}
        } else {arrdim=rev(attr(gene.expression.vector,'sizes'))}
     } else {
        arrdim = dim(gene.expression)
     }
     
     if (d1) {
        # put into an array. reverse sizes attribute for array dimensions
        gene.expression.array=array(gene.expression.vector,arrdim)
     } else {
        gene.expression.array = gene.expression
     }

     if (attr(gene.expression,'orientation') == 'ABI') {
        if (!is.null(reflect.dim)) {
           if (reflect.dim != 1) {
              warning(
                 'gene.expression is in Allen orientation',
                 'but you are not reflecting across sagittal midplane.')
           }
        } else { reflect.dim = 1 }
     } else if (attr(gene.expression,'orientation') == 'MINC') {
        if (!is.null(reflect.dim)) {
           if (reflect.dim != 3) {
              warning(
                  'gene.expression is in MINC orientation',
                  'but you are not reflecting accross sagittal midplane.')
           }
        } else { reflect.dim = 3 }
     }
     # must reverse dimension to reflect
     # reflect.dim=1 --> rev.reflect.dim=3
     # reflect.dim=2 --> rev.reflect.dim=2
     # reflect.dim=3 --> rev.reflect.dim=1
     rev.reflect.dim=4-reflect.dim
     
     # find midplane
     ref.dim.halflength=floor(median(1:arrdim[rev.reflect.dim]))
     
     #Find indices to split array (depending on which dimension to reflect)
     #define index order to be closest to the mirror
     idx11 = ref.dim.halflength
     idx12 = 1
     idx21 = ifelse(arrdim[rev.reflect.dim] %% 2 == 1, ref.dim.halflength, ref.dim.halflength + 1)
     idx22 = arrdim[rev.reflect.dim]

     if (rev.reflect.dim==3) {
          idx1=paste0(",,",idx11,":",idx12)
          idx2=paste0(",,(",idx21,"):",idx22)
     } else if (rev.reflect.dim==2) {
          idx1=paste0(",",idx11,":",idx12,",")
          idx2=paste0(",(",idx21,"):",idx22,",")
     } else if (rev.reflect.dim==1) {
          idx1=paste0(idx11,":",idx12,",,")
          idx2=paste0("(",idx21,"):",idx22,",,")}
     
     
     genearr1=eval(parse(text=paste0("gene.expression.array[",idx1,"]")))
     genearr2=eval(parse(text=paste0("gene.expression.array[",idx2,"]")))
     
     if (data == 'gene') {
        ## Find where one half is missing data (missing data are elements with -1)
        ## Replace it with the same voxels in the other half
        bool=genearr1==-1
        genearr1[bool]=genearr2[bool]
        
        ## DO the same with the other half
        bool=genearr2==-1
        genearr2[bool]=genearr1[bool]
     } else if (data == 'atlas') {
        ## Find which side has a lower signal
        mn1 = mean(genearr1)
        mn2 = mean(genearr2)
        if (mn1>mn2) {   genearr2 = genearr1
        } else {         genearr1 = genearr2 }
     }
     
     
     ## Put reflected back into a return array
     ret.array=gene.expression.array*0
     eval(parse(text=paste0("ret.array[",idx1,"]=genearr1")))
     eval(parse(text=paste0("ret.array[",idx2,"]=genearr2")))

     ## Output as vector with the old attributes     
     if (d1) { ret = as.vector(ret.array)
     } else  { ret = ret.array }
     
     attributes(ret) = attributes(gene.expression)
     return(ret)
}