R/multichannelPCA.R

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

#' Shape analysis of multichannel images based on PCA
#'
#' Converts a set of multichannel images (e.g. deformation fields ) to a matrix
#' to enable various forms of PCA.  Returns the components of shape variability
#' and variance explained.  May employ different decomposition methods (WIP).
#'
#' @param x list containing multichannel images all of the same size
#' @param mask mask to apply to the multichannel images
#' @param k rank to use
#' @param pcaOption currently only PCA and randPCA, latter being much faster.
#' We also allow fastICA.
#' @param auxiliaryModality if you pass this matrix, then will do CCA.  This
#' will only work with one option.
#' @param center subtract the mean column vector.
#' @param sigma parameter for kernel PCA.
#' @param verbose produces more explanatory output.
#' @return list of the pca output and conversion to multichannel images
#' @author Avants BB
#' @examples
#'
#' img1 = antsImageRead( getANTsRData( "r16" ) ) %>%
#'   resampleImage( c(4,4) )
#' img2 = antsImageRead( getANTsRData( "r64" ) ) %>%
#'   resampleImage( c(4,4) )
#' img3 = antsImageRead( getANTsRData( "r27" ) ) %>%
#'   resampleImage( c(4,4) )
#' img4 = antsImageRead( getANTsRData( "r30" ) ) %>%
#'   resampleImage( c(4,4) )
#' reg1 = antsRegistration( img1, img2, 'SyN' )
#' reg2 = antsRegistration( img1, img3, 'SyN' )
#' reg3 = antsRegistration( img1, img4, 'SyN' )
#' w1 = antsImageRead( reg1$fwdtransforms[1] )
#' w2 = antsImageRead( reg2$fwdtransforms[1] )
#' w3 = antsImageRead( reg3$fwdtransforms[1] )
#' mask = getMask( img1 )
#' x = list( w1, w2, w3 )
#' dpca = multichannelPCA( x, mask )
#' warpTx = antsrTransformFromDisplacementField( dpca$pcaWarps[[1]] )
#' warped = applyAntsrTransform( warpTx, data = img1, reference = img1)
#'
#' @export multichannelPCA
multichannelPCA <- function(
  x,
  mask,
  k = NA,
  pcaOption = "PCA",
  auxiliaryModality,
  center = TRUE,
  sigma = NA,
  verbose = FALSE ) {
n = length( x )
idim = mask@dimension
p   = sum( mask >= 1 )
# check the size of the inputs matches the mask
for ( i in 1:n )
  if ( !all( dim( x[[n]] ) == dim(mask)  ) )
    stop( paste("dimensionality of input number",i,
    "does not match the mask dimensionality."))
nchannels = x[[ 1 ]]@components
vecmat = matrix( nrow = n, ncol = p*nchannels )
for ( i in 1:n )
  {
  vecmat[ i, ] = multichannelToVector( x[[ i ]] , mask )
  }
  if ( is.numeric( pcaOption ) ) {
    pcak = pcaOption
    pcaOption = "kPCA"
  }
  if ( verbose ) print( paste( "begin decomposition option",pcaOption) )
  if ( is.na( k ) ) k = nrow( vecmat ) - 1
  if ( center ) cx   = sweep( vecmat, 2, colMeans(vecmat), "-") else cx=vecmat
  if ( pcaOption == "randPCA" ) {
    if ( ! usePkg( "rsvd" ) ) stop("please install rsvd")
    vpca = rsvd::rsvd( cx, k )
    } else if ( pcaOption == "rrPCAL" ) {
    vpca = rsvd::rrpca( cx, k=k  )
    vpca = list( d=NA,  u = cx %*% t(vpca$L), v=t( vpca$L ) )
    } else if ( pcaOption == "rrPCAS" ) {
    vpca = rsvd::rrpca( cx, k=k  )
    vpca = list( d=NA,  u = cx %*% t(vpca$S), v=t( vpca$S ) )
    } else if ( pcaOption == "kPCA" ) {
      kpcaopt = 'cov'
      if ( ! is.na( sigma ) ) kpcaopt = 'gaussian'
      if ( ! usePkg( "irlba" ) ) stop("please install irlba")
      if ( missing( "auxiliaryModality") )
        tempdistmat = sparseDistanceMatrix( cx, pcak, kmetric=kpcaopt, sigma=sigma )
      if ( !missing( "auxiliaryModality") )
        {
        if ( center ) cy = sweep( auxiliaryModality, 2, colMeans(auxiliaryModality), "-")
        if ( !center ) cy = auxiliaryModality
        tempdistmat = sparseDistanceMatrixXY( cy, cx, pcak, kmetric=kpcaopt, sigma=sigma )
        }
      vpca = irlba::irlba( tempdistmat, nu=k, nv=k )
      vpca$u = cx %*% vpca$v
    } else if ( pcaOption == "fastICA" ) {
      if ( ! usePkg( "fastICA" ) ) stop("please install fastICA")
      tempica = fastICA::fastICA( t( cx ), k )
      vpca = list( d=NA,  u = cx %*% tempica$S, v=( tempica$S ) )
    } else if ( pcaOption == "eanat" ) {
      # FIXME - implement mask for regularization
      # need to bind mask in proper order to make
      # regularization work
      arr = abind::abind( as.array( mask ), as.array( mask ), along=idim )
      if ( idim == 3 )
        arr = abind::abind( arr,  as.array( mask ), along=idim )
      maska = as.antsImage( arr )
      a = antsCopyImageInfo( mask, maska )
#      eanat = sparseDecom( cx, inmask=maska, sparseness = 1.0/k, nvecs=k,
#        cthresh=10, smooth=0.5, verbose=TRUE )
#      vpca = list( d=eanat$varex,  u=eanat$umatrix,
#       v=t(eanat$eigenanatomyimages ) )
      eanatD = eanatDef( cx, mask = maska, smoother = mean(antsGetSpacing(mask)),
        positivity=TRUE, nvecs = k, cthresh=0, verbose=TRUE )
      vpca = list( d=NA,  u=NA, v=t( eanatD ) )
      k = ncol( vpca$v )
    } else {
      vpca = svd(cx, nv=k, nu=k )
    }
  if ( !verbose ) { rm( vecmat ); vecmat=NA }
  if ( verbose ) {
    print( paste( "convert back to multichannel" ) )
    }
  # now convert the vectors back to warps
  pcaWarps = list( )
  for ( i in 1:k )
    {
    pcaWarps[[ i ]] = vectorToMultichannel( vpca$v[,i], mask )
    }
  datatopcacorrs = NA
  mylms = NA
  if ( verbose ) print( paste( "regression and correlation" ) )
  if ( verbose )
    {
    datatopcacorrs = cor( t( vecmat ) , vpca$v )
    # can also explore regressing the basis against the data
    mydf = data.frame( vpca$v )
    mylms = summary(lm( t( vecmat ) ~ .,data=mydf))
    }
  gc( )
  return( list(
    pca = vpca,
    pcaWarps = pcaWarps,
    vecmat = vecmat,
    datatopcacorrs = datatopcacorrs,
    mylms = mylms ) )
}


#' Convert multichannel \code{antsImage} to a vector.
#'
#' Employs a mask to flatten each channel of a multi-channel image into a
#' contiguous piece of a vector.
#'
#' @param multichannelimage input multichannel image
#' @param mask mask of same dimensionality as multichannelimage
#' @return vector is output
#' @author Avants BB
#' @examples
#'
#' # see vectorToMultichannel
#'
#' @export multichannelToVector
multichannelToVector <- function( multichannelimage, mask )
{
  dd = mask@dimension
  p = sum( mask >= 1 )
  nchannels = multichannelimage@components
  temp = splitChannels( multichannelimage )
  maxn = (nchannels*p)
  myinds = seq( from=1, to=maxn, by=(p) )
  myinds[ nchannels+1 ] = maxn
  v = rep( 0, maxn )
  for ( k in 1:nchannels ) {
    maxind = (myinds[k+1]-1)
    if ( k == nchannels ) maxind = maxn
    locinds = myinds[k]: maxind
    v[ locinds ] = temp[[k]][ mask >= 1 ]
    }
  rm( temp )
  gc()
  return( v )
}


#' Convert vector to multichannel \code{antsImage}.
#'
#' Converts a vector of size n-entries in mask times number of channels to
#' a multi-channel image with the same dimensionality as the mask.
#'
#' @param v input multichannel vector
#' @param mask mask of same dimensionality as multichannelimage
#' @return multichannelimage is output
#' @author Avants BB
#' @examples
#'
#' fi <- antsImageRead(getANTsRData("r16") ) %>%
#'   resampleImage(c(60,60),1,0)
#' mi <- antsImageRead(getANTsRData("r64") ) %>%
#'   resampleImage(c(60,60),1,0)
#' mytx <- antsRegistration(fixed=fi, moving=mi, typeofTransform = c('SyN') )
#' mcimg = antsImageRead( mytx$fwd[1] )
#' msk = getMask( fi )
#' vv=multichannelToVector(mcimg,msk)
#' mcimg2=vectorToMultichannel( vv, msk )
#' vv2=multichannelToVector(mcimg2,msk)
#' cor.test(vv2,vv)
#' stopifnot( all( mcimg2[30,30] == mcimg[30,30] ) )
#'
#' @export vectorToMultichannel
vectorToMultichannel <- function( v, mask ) {
  dd = mask@dimension
  p = sum( mask >= 1 )
  nchannels = round( length( v ) / p )
  if ( length( v ) != (nchannels*p) ) stop("dimensions do not match")
  maxn = (nchannels*p)
  myinds = seq( from=1, to=maxn, by=(p) )
  myinds[ nchannels+1 ] = maxn
  mylist = list( )
  for ( k in 1:nchannels )
    {
    maxind = (myinds[k+1]-1)
    if ( k == nchannels ) maxind = maxn
    locinds = myinds[k]: maxind
    temp = antsImageClone( mask )
    temp[ mask >= 1 ] = v[ locinds ]
    mylist[[ k ]] = temp
    }
  vecimg = mergeChannels( mylist )
  rm( temp, mylist )
  gc()
  return( vecimg )
}