R/common.R

In bbuchsbaum/neuroim: Data Structures and Handling for Neuroimaging Data

#' matrixToVolumeList
#' converts a matrix to a list of BrainVolumes with values filled at grid coordinates determined by the \code{vox} argument.
#' @param voxmat an N by 3 matrix of voxel coordinates
#' @param mat an N by M matrix of values where M is the number of volumes to create (e.g. one volume per column in \code{mat})
#' @param mask a reference volume defining the geometry of the output volumes. This can either be of type \code{BrainSpace} or \code{BrainVolume}
#' @param default the value that will be used for voxels not contained within voxmat (defualt is \code{NA})
#' @return a \code{list} of \code{BrainVolume} instances, one for each column of \code{mat}
#' @export  
matrixToVolumeList <- function(voxmat, mat, mask, default=NA) {
  if (nrow(voxmat) != nrow(mat)) {
    stop("mismatching dimensions: nrow(voxmat) must equal nrow(mat)")
  }
  
  lapply(1:ncol(mat), function(i) {
    vol <- array(default, dim(mask))   
    vol[voxmat] <- mat[,i]
    
    if (is(mask, "BrainSpace")) {
      BrainVolume(vol, mask)
    } else {
      BrainVolume(vol, space(mask))
    }
  })
}  

matrixToVols <- matrixToVolumeList

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "matrix", fac="integer", FUN="function"),
          def=function(x, fac, FUN) {
            callGeneric(x,as.factor(fac), FUN)
          })

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "matrix", fac="integer", FUN="missing"),
          def=function(x, fac) {
            callGeneric(x,as.factor(fac))
          })

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "matrix", fac="factor", FUN="missing"),
          def=function(x, fac) {
            if (length(fac) != nrow(x)) {
              stop(paste("x must be same length as factor used for splitting rows"))
            }
            
            ind <- split(seq_along(fac), fac)
            out <- do.call(rbind, lapply(levels(fac), function(lev) {
              colMeans(x[ind[[lev]],])
            }))
            
            row.names(out) <- levels(fac)           
            out
          })

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "matrix", fac="factor", FUN="function"),
          def=function(x, fac, FUN) {
            if (length(fac) != nrow(x)) {
              stop(paste("x must be same length as split variable"))
            }        
            
            ind <- split(seq_along(fac), fac)
            out <- do.call(rbind, lapply(names(ind), function(lev) {
              apply(x[ind[[lev]],], 2, FUN)
            }))
            
            row.names(out) <- levels(fac)           
            out
          })

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "BrainVector", fac="factor", FUN="function"),
          def=function(x, fac, FUN) {
            if (length(fac) == prod(dim(x)[1:3])) {
              split_by_voxel <- TRUE
            } else if (length(fac) == dim(x)[4]) {
              split_by_row <- TRUE
            } else {
              stop(paste("length of 'fac' must be equal to number of voxels or to number of volumes"))
            }
            
            if (split_by_voxel) {
            
              ind <- split(seq_along(fac), fac)
              out <- do.call(rbind, lapply(names(ind), function(lev) {
                #idx <- which(fac == lev)
                mat <- series(x, ind[[lev]])
                apply(mat, 1, FUN)
              }))
            
              row.names(out) <- levels(fac)           
              out
            } else {
              m <- as.matrix(x)
              callGeneric(m, fac, FUN)
            }
          })

#' @export
#' @rdname splitReduce-methods
setMethod(f="splitReduce", signature=signature(x = "BrainVector", fac="factor", FUN="missing"),
          def=function(x, fac, FUN) {
            if (length(fac) == prod(dim(x)[1:3])) {
              split_by_voxel <- TRUE
            } else if (length(fac) == dim(x)[4]) {
              split_by_row <- TRUE
            } else {
              stop(paste("length of 'fac' must be equal to number of voxels or to number of volumes"))
            }
            
            if (split_by_voxel) {
            
              ind <- split(seq_along(fac), fac)
              out <- do.call(rbind, lapply(names(ind), function(lev) {
                rowMeans(series(x, ind[[lev]]))
              }))
            
              row.names(out) <- levels(fac)           
              out
            } else {
              m <- as.matrix(x)
              callGeneric(m, fac)
            }
          })



#' @export
#' @rdname splitScale-methods
setMethod(f="splitScale", signature=signature(x = "matrix", f="factor", center="logical", scale="logical"),
          def=function(x, f, center=TRUE, scale=TRUE) {
            if (length(f) != nrow(x)) {
              stop(paste("x must be same length as split variable"))
            }
            
            out <- matrix(0, nrow(x), ncol(x))
            ind <- split(seq_along(f), f)
            
            for (lev in names(ind)) {
              keep <- ind[[lev]]
              xs <- scale(x[keep,,drop=FALSE], center=center, scale=scale)
              out[keep,] <- xs         
            }
            
            out
          })


#' @export
#' @rdname splitScale-methods
#' @importFrom abind abind
setMethod(f="splitScale", signature=signature(x = "matrix", f="factor", center="missing", scale="missing"),
          def=function(x, f) {
            callGeneric(x,f, TRUE, TRUE)
          })


#' .isExtension
#' @rdname internal-methods
#' @keywords internal
.isExtension <- function(fname, extension) {
  last <- substr(fname, nchar(fname)+1 - nchar(extension), nchar(fname))
  return(last==extension)
}

#' .concat4D
#' @rdname internal-methods
#' @keywords internal
.concat4D <- function(x,y, ...) {
	rest <- list(...)
	
	D <- dim(x)[1:3]
	
	rvols <- lapply(rest, function(z) {
			stopifnot(length(dim(z)) >= 3)
			stopifnot(identical(D, dim(z)[1:3]))
			z
	})
	
	out <- cbind(as.matrix(x), as.matrix(y))
	
	if (length(rvols) > 0) {
	  out2 <- do.call(cbind, lapply(rvols, as.matrix))
	  out <- cbind(out, out2)
	}

	nvols <- ncol(out)
	new.dim <- c(D, nvols)
	
	nspace <- BrainSpace(new.dim, origin=origin(x@space), spacing=spacing(x@space),
			axes=axes(x@space), trans=trans(x@space))
	
	DenseBrainVector(out, nspace)
	
}


#' .gridToIndex3D
#' @rdname internal-methods
#' @importFrom assertthat assert_that
#' @keywords internal
.gridToIndex3D <- function(dimensions, voxmat) {
	assert_that(length(dimensions) == 3)
  if (is.vector(voxmat)) {
    assert_that(length(voxmat) == 3)
    voxmat <- matrix(voxmat, 1,3)
  }
  
  assert_that(ncol(voxmat) == 3)
  gridToIndex3DCpp(dimensions, voxmat)

}

#' .gridToIndex
#' @rdname internal-methods
#' @keywords internal
.gridToIndex <- function(dimensions, vmat) {
	D <- Reduce("*", dimensions, accumulate=TRUE)
	apply(vmat, 1, function(vox) {
		sum(sapply(length(D):2, function(i) {
			D[i-1]*(vox[i]-1)
		})) + vox[1]
	})
	
}

#' .indexToGrid
#' @rdname internal-methods
#' @keywords internal
.indexToGrid <- function(idx, array.dim) {
  assert_that(all(idx > 0 & idx <= prod(array.dim)))
  assert_that(length(array.dim) <= 5)
  indexToGridCpp(idx, array.dim)
  
}



#' .getRStorage
#' @rdname internal-methods
#' @keywords internal
.getRStorage <- function(dataType) {
  if (any(toupper(dataType) == c("BINARY", "BYTE", "UBYTE", "SHORT", "INTEGER", "INT", "LONG"))) {
    "integer"
  } else if (any(dataType == c("FLOAT", "DOUBLE"))) {
    "double"
  } else {
	  stop(paste("unrecognized data type", dataType))
  }
}

# @nord
#' @importFrom mmap int8 uint8 int16 int32 real32 real64
#' @importFrom mmap mmap char mmapFlags munmap
.getMMapMode <- function(code) {	
	if (code == "UNKNOWN") {
		stop(paste(".getMMapMode: no memory map mode for UNKNOWN data type: ", code))
	} else if (code == "BINARY") {
		mmap::int8()
	} else if (code == "UBYTE") {
	  mmap::uint8()
	} else if(code == "SHORT") {
	  mmap::int16()
	} else if(code == "INT") {
	  mmap::int32()
	} else if (code == "FLOAT") {
	  mmap::real32()
	} else if (code == "DOUBLE") {
	  mmap::real64()
	} else {
		stop(paste(".getMMapMode: unsupported data type: ", code))
	}
}
	

#' .getDataStorage
#' @rdname internal-methods
#' @keywords internal
.getDataStorage <- function(code) {
  if (code == 0) {
    return("UNKNOWN")
  } else if (code == 1) {
    return("BINARY")
  } else if (code == 2) {
    return("UBYTE")
  } else if(code == 4) {
    return("SHORT")
  } else if(code == 8) {
    return("INT")
  } else if (code == 16) {
    return("FLOAT")
  } else if (code == 64) {
    return("DOUBLE")
  } else {
    stop(paste("nifti(getDataStorage): unsupported data type: ", code))
  }
}

#' .getDataCode
#' @rdname internal-methods
#' @keywords internal
.getDataCode <- function(dataType) {
  if (dataType == "UNKNOWN") {
    return(0)
  }else if (dataType == "BINARY") {
    return(1)
  } else if (dataType == "UBYTE") {
    return(2)
  } else if(dataType == "SHORT") {
    return(4)
  } else if(dataType == "INT") {
    return(8)
  } else if (dataType == "FLOAT") {
    return(16)
  } else if (dataType == "DOUBLE") {
    return(32)
  } else {
    stop(paste("getDataCode: unsupported data type: ", dataType))
  }
}

#' .getDataSize
#' @rdname internal-methods
#' @keywords internal
.getDataSize <- function(dataType) {
  if (dataType == "BINARY") {
    return(1)
  } else if (dataType == "UBYTE") {
	  return(1)
  } else if (dataType == "UBYTE") {
    return(1)
  } else if (dataType == "SHORT") {
    return(2)
  } else if (dataType == "INTEGER") {
    return(4)
  } else if (dataType == "INT") {
    return(4)
  } else if (dataType == "FLOAT") {
    return(4)
  } else if (dataType == "DOUBLE") {
    return(8)
  } else if (dataType == "LONG") {
    return(8)
  }

  stop(paste("unrecognized data type: ", dataType))
}

#' .getEndian
#' @rdname internal-methods
#' @keywords internal
.getEndian <- function(conn) {
  #try little endian
  endian <- "little"

  hsize <- readBin(conn, integer(), 1, endian=endian)
  if (hsize != 348) {
    # might be bigendian
    endian <- "big"
    seek(conn, 0)
    hsize <- readBin(conn, integer(), 1, endian=endian)
    if (hsize != 348) {
      stop("nifti(getEndian): header size is not 348, invalid header.")
    }
  }

  return(endian)
}

#' @rdname internal-methods
#' @name .niftiExt
#' @keywords internal
.niftiExt <- function(filetype) {

  extensions <- list()

  if (filetype == "nifti-single") {
    extensions[["header"]]  <- "nii"
    extensions[["data"]] <- "nii"
  }
  else if (filetype == "nifti-pair") {
    extensions[["header"]]  <- "hdr"
    extensions[["data"]] <- "img"
  }
  else if (filetype == "nifti-gz") {
    extensions[["header"]]  <- "nii.gz"
    extensions[["data"]] <- "nii.gz"
  } else {
    stop(paste("unsupported filetype: ", filetype))
  }

  return(extensions)
}

#' @rdname internal-methods
#' @keywords internal
.matrixToQuatern <- function(mat) {
  xd <- sqrt(drop(crossprod(mat[1:3,1])))
  yd <- sqrt(drop(crossprod(mat[1:3,2])))
  zd <- sqrt(drop(crossprod(mat[1:3,3])))

  if (xd == 0) { mat[1,1] = 1; mat[2:3,1] = 0; xd = 1; }
  if (yd == 0) { mat[2,2] = 1; mat[c(1,3),2] = 0; yd = 1; }
  if (zd == 0) { mat[3,3] = 1; mat[1:2,3] = 0; zd = 1; }

  rmat = mat[1:3, 1:3]
  rmat[,1] = rmat[,1]/xd
  rmat[,2] = rmat[,2]/yd
  rmat[,3] = rmat[,3]/zd

  ####### skipping orthogonalization of columns

  #################################################

  zd = det(rmat)
  qfac = 1

  if (zd > 0) {
    qfac = 1
  } else {
    qfac = -1
    rmat[1:3,3] = -rmat[1:3,3]
  }

  # compute quaternion parameters

  a = rmat[1,1] + rmat[2,2] + rmat[3,3] + 1

  if (a > .5) {
    a = .5 * sqrt(a)
    b = 0.25 * (rmat[3,2]-rmat[2,3]) / a
    c = 0.25 * (rmat[1,3]-rmat[3,1]) / a
    d = 0.25 * (rmat[2,1]-rmat[1,2]) / a
   } else {
     xd = 1.0 + rmat[1,1] - (rmat[2,2]+rmat[3,3])
     yd = 1.0 + rmat[2,2] - (rmat[1,1]+rmat[3,3])
     zd = 1.0 + rmat[3,3] - (rmat[1,1]+rmat[2,2])
     if( xd > 1.0 ){
       b = 0.5 * sqrt(xd)
       c = 0.25* (rmat[1,2]+rmat[2,1]) / b
       d = 0.25* (rmat[1,3]+rmat[3,1]) / b
       a = 0.25* (rmat[3,2]-rmat[2,3]) / b
     } else if( yd > 1.0 ){
       c = 0.5 * sqrt(yd) ;
       b = 0.25* (rmat[1,2]+rmat[2,1]) / c
       d = 0.25* (rmat[2,3]+rmat[3,2]) / c
       a = 0.25* (rmat[1,3]-rmat[3,1]) / c
     } else {
       d = 0.5 * sqrt(zd) ;
       b = 0.25* (rmat[1,3]+rmat[3,1]) / d
       c = 0.25* (rmat[2,3]+rmat[3,2]) / d
       a = 0.25* (rmat[2,1]-rmat[1,2]) / d
     }
     if( a < 0.0 ){ b=-b ; c=-c ; d=-d; a=-a; }
   }

  return(list(quaternion=c(b,c,d), qfac=qfac))

}


.quaternToMatrix <- function(quat, origin, stepSize, qfac) {
  mat <- matrix(0, 4,4)
  mat[4,] <- c(0,0,0,1)

  a <- 1 - sum(quat^2)
  if (a < 1e-07) {
    a <- 1 /(sqrt(sum(quat^2)))
    quat <- quat*a
    a <- 0
  } else {
    a <- sqrt(a)
  }

  stepSize <- ifelse(stepSize > 0, stepSize, 1)
  xd <- stepSize[1]
  yd <- stepSize[2]
  zd <- stepSize[3]

  if (qfac < 0) {
    zd <- -zd
  }

  b <- quat[1]
  c <- quat[2]
  d <- quat[3]


  mat[1,1] <- (a*a+b*b-c*c-d*d) * xd
  mat[1,2] <- 2 * (b*c-a*d)     * yd
  mat[1,3] <- 2 * (b*d+a*c)     * zd
  mat[2,1] <- 2 * (b*c+a*d)     * xd
  mat[2,2] <- (a*a+c*c-b*b-d*d) * yd
  mat[2,3] <- 2 * (c*d-a*b)     * zd
  mat[3,1] <- 2 * (b*d-a*c)     * xd
  mat[3,2] <- 2 * (c*d+a*b)     * yd
  mat[3,3] <- (a*a+d*d-c*c-b*b) * zd

  mat[1:3,4] <- origin

  return(mat)
}
# 
# .byteSwapppedReader <- function(rawmap, meta, offset) {
#   dsize <- .getDataSize(meta@dataType)
#   nels <- prod(meta@Dim)
#   dtype <- .getRStorage(meta@dataType)
#   getElements <- function(indices) {
#     
#   }
# }

.makeMMap <- function(meta) {
  nels <- prod(meta@Dim[1:4]) 	
  
  if (.Platform$endian != meta@endian) {
    ## read raw bytes
    rawbytes <- mmap::mmap(meta@dataFile, mode=mmap::char(), prot=mmap::mmapFlags("PROT_READ"))
    rawbytes <- rawbytes[(meta@dataOffset+1):length(rawbytes)]
    
    mmap::munmap(rawbytes)
    readBin(rawbytes, what=.getRStorage(meta@dataType), size=.getDataSize(meta@dataType), n=nels, endian=meta@endian)
  } else {
    #mmap::mmap(meta@dataFile, mode=.getMMapMode(meta@dataType), off=meta@dataOffset,prot=mmap::mmapFlags("PROT_READ"),flags=mmap::mmapFlags("MAP_PRIVATE"))
    ret <- mmap::mmap(meta@dataFile, mode=.getMMapMode(meta@dataType), prot=mmap::mmapFlags("PROT_READ"))
    offset <- meta@dataOffset/.getDataSize(meta@dataType) + 1
    vals <- ret[offset:nels]
    mmap::munmap(ret)
    vals
  }
  
  
}