R/sph.odfvmf.R

In gdimap: Generalized Diffusion Magnetic Resonance Imaging

sph.odfvmf <-
function(run=TRUE, fbase=NULL, savedir=tempdir(), rg=NULL, swap=FALSE, btoption=2,
 threshold=0.4, showglyph=FALSE, bview="coronal",  order=4, clusterthr=0.6,  aniso=NULL, ...)
{
  bviews <- c("sagittal", "coronal", "axial")
  kv <- match(bview, bviews)
  stopifnot(is.na(kv) != TRUE)
  ## control parameters for movMF
  E <- list(...)[["E"]]
  if (is.null(E)) E <- "softmax"
  kappa <- list(...)[["kappa"]]
  if (is.null(kappa)) kappa <- "Newton_Fourier"
  minalpha <- list(...)[["minalpha"]]
  if (is.null(minalpha)) minalpha <- 8
  start <- list(...)[["start"]]
  if (is.null(start)) start <- "s"
  startctl=list(E=E, kappa=kappa, minalpha=minalpha, start=start) ## movMF inits
  ##-----------
  ## Read data
  testfilexist(fbase=fbase, btoption=btoption)
  if(btoption == 1){ ## Option 1: S2-shell (DSI 203-point 3mm)
    btable <- as.matrix(readtable(fbase=fbase, filename="btable.txt"))
  } else {
    if(btoption == 2) { ## Option 2: 3D-dsi grid 
      bval <- scantable(fbase=fbase, filename="data.bval")
      # bvec <- readtable(fbase=fbase, filename="data.bvec")
      bvec <- scantable(fbase=fbase, filename="data.bvec")
      bvec <- matrix(bvec, ncol=3)
      btable <- cbind(bval,bvec)
      rm(bval, bvec)
    }
    else stop()
  }
  b0 <- which(btable[,1] == 0)
  odfvertices <- matrix(btable[-b0,2:4], ncol=3)
  tc <-  geometry::delaunayn(odfvertices)
  tcsurf <- t( surf.tri(odfvertices,tc))  
  ##----------------------------
  gc()
  cat("Reading data ...")
  ptm <- proc.time()
  img.nifti  <- readniidata(fbase=fbase, filename="data.nii.gz")
  volimg <- img.nifti@.Data  
  mask.nifti <- readniidata(fbase=fbase, filename="data_brain_mask.nii.gz")
  volmask <- mask.nifti@.Data  
  print(proc.time() - ptm)
  rm(img.nifti, mask.nifti)
  gc()
  ##----------------------------
  d <- dim(volmask)
  volgfa <- array(0, dim=d)   ## gfas map
  V1 <- array(0, dim=c(d, 3)) ## V1 direction
  V2 <- array(0, dim=c(d, 3)) ## V2 direction
  V3 <- array(0, dim=c(d, 3)) ## V3 direction
  if(is.null(rg)) {
    switch(kv,
      { nslices <- d[1]}, # sagittal,
      { nslices <- d[2]}, # coronal
      { nslices <- d[3]})  # axial
    first <- 1; last <- nslices
  }
  else { first <- rg[1]; last <- rg[2] }
  cat("\n")
  ##-----------------------------
  ## SPH process preparation
  gradient <- t(odfvertices)
  z <- design.spheven(order,gradient,lambda=0.006)
  plz <- plzero(order)/2/pi
  ngrad <- dim(gradient)[2]
  ngrad0 <- ngrad
  lord <- rep(seq(0,order,2),2*seq(0,order,2)+1)
  while(length(lord)>=ngrad0){
    order <- order-2
    lord <- rep(seq(0,order,2),2*seq(0,order,2)+1)
     cat("Reduced order of spherical harmonics to",order,"\n")
   }
  cat("Using",length(lord),"spherical harmonics\n")
  L <- -diag(lord*(lord+1)) 
  ##-----------------------------
  ## store 1st vector directions for each non-thresholded voxel 
  ## v1list: vector of lists
  nv1 <- length(first:last)
  v1list <- vector(mode="list", nv1)
  v1count <- 0
  # rglinit()
  npar1 <- 7
  npar2 <- 15
  npar3 <- 23
  for (sl in (first:last)) {
    cat("slice",sl,"\n")
    ptm <- proc.time()
    if(run) {
      slicedata <- read.slice(img=volimg, mask=volmask, slice=sl,
       swap=swap, bview=bview)
      ymaskdata <- premask(slicedata)
      if(ymaskdata$empty) next # empty mask
      maxslicedata <- max(slicedata$niislicets) ##????
      S <- ymaskdata$yn[-b0,]
       S <- S / maxslicedata
      s0 <- 1
      si <- apply(S, 2, datatrans, s0)
      sicoef <- z$matrix%*% si
      sphcoef <- plz%*%L%*%sicoef
      coef0 <- sphcoef[1,]
      sphcoef[1,] <- 1/2/sqrt(pi)
      sphcoef[-1,] <- sphcoef[-1,]/8/pi
      ## odfs
      odfs <- t(z$design) %*% sphcoef
      # odfs <- apply(odfs, 2, norm01) 
    	odfs <- apply(odfs, 2, anisofn, aniso=aniso) 
      ## gfas
      gfas <- apply(odfs, 2, genfa)
      gfas <- norm01(gfas) ## ??
      z2d <- ymaskdata$kin
      ## mask out thresholded values
      zx <- which(gfas <= threshold)
      if(length(zx)) {
        z2d <- z2d[-zx,]
        gfas <- gfas[-zx]
        odfs <- odfs[,-zx]
      }
      if(is.null(dim(z2d))) next
      if(length(gfas) < 8) next # 8 elements as minimum number
      lix <- dim(z2d)[1]
      v1perslice <- matrix(0, nrow=lix,ncol=3) # v1 directions 
      v2perslice <- matrix(0, nrow=lix,ncol=3) # v2 directions 
      v3perslice <- matrix(0, nrow=lix,ncol=3) # v3 directions 
      nullvectors <- NULL
      tperc <- c(20, 40, 60, 80)
      tline <- floor(c(0.2,0.4,0.6,0.8)*lix) 
      cat("processing", lix,"voxels \n")
      for(m in 1:lix) {
        tt <- which(tline == m)
        if(length(tt) != 0) {
          cat(paste(tperc[tt],"% ", sep="")); cflush()
        }
        odf <- odfs[,m]
        ## Find peaks based on clusters
        ith <- which(odf < clusterthr) 
        vx <- odfvertices[-ith,]
        n <- dim(vx)[1]
        lbn <- logb(n)
        ## Fit a vMF mixture  with k=2
        y1 <- movMF::movMF(vx, k=2, control=startctl)
        ## Inspect the fitted parameters:
        par1 <- lbn*npar1
        bic1 <- 2*logLik(y1) - par1
        ## Fit a vMF mixture  with k=4
        y2 <- movMF::movMF(vx, k=4, control=startctl)
        par2 <- lbn*npar2
        bic2 <- 2*logLik(y2) - par2
        if(bic1 >= bic2) { yy <- y1 }
        else {
          ## Fit a vMF mixture  with k=6
          y3 <- movMF::movMF(vx, k=6, control=startctl) 
          par3 <- lbn*npar3
          bic3 <- 2*logLik(y3) - par3
          if(bic2 >= bic3) { yy <- y2 }
          else { yy <- y3 }
        }
        np <- dim(yy$theta)[1]
        if(np < 2) {
          nullvectors <- c(nullvectors, m)
          next
        }
        # reorder by alpha weight
        yys <- sort(yy$alpha, decreasing=TRUE, index.return=TRUE)
        yyn <- yy$theta[yys$ix,]
        yy <- list(theta=yyn, alpha= yys$x)
        pk <- list(np=np , pcoords=t(yy$theta))
        v1perslice[m,] <- pk$pcoords[,1]
        if(np >= 4) 
          v2perslice[m,] <- pk$pcoords[,3]
        if(np == 6) 
          v3perslice[m,] <- pk$pcoords[,5]
        if(showglyph) {
          if(pk$np > 2) { ## show crossing-fiber glyphs only
            # normalize for visualization
            if(length(yy$alpha) < 2)
              pn <- fnorm(yy$theta)
            else
              pn <- apply(yy$theta, 1, fnorm)
            pcoords=yy$theta/pn
            pkv <- list(np=np , pcoords=t(pcoords))
            plotglyph(odfs[,m], odfvertices, pkv, kdir=6)
            pp <- readline(
                "\nmore crossing-fiber glyphs  ? ('n' to exit) ") 
            if(pp == "n" ) { showglyph <- FALSE; }
            else { rgl.clear( type = "shapes" ) }
          }
        }
      }
      cat("100% completed\n")
      ## remove null pk vectors
      nvl <- lix
      nnv <- length(nullvectors)
      if(nnv > 0) {
        nvl <- nvl-nnv
        v1perslice <- v1perslice[-nullvectors,]
        v2perslice <- v2perslice[-nullvectors,]
        v3perslice <- v3perslice[-nullvectors,]
        z2d <- z2d[-nullvectors,]
        gfas <- gfas[-nullvectors]
      }
      if(is.null(dim(z2d))) next
      ## V volumes
      for(k in 1:3) {
        switch(kv,
        { mx <- matrix(0, d[2],d[3])
        mx[z2d] <- v1perslice[,k]
        V1[sl,,,k] <- mx
        mx <- matrix(0, d[2],d[3])
        mx[z2d] <- v2perslice[,k]
        V2[sl,,,k] <- mx 
        mx <- matrix(0, d[2],d[3])
        mx[z2d] <- v3perslice[,k]
        V3[sl,,,k] <- mx }, # sagittal
        { mx <- matrix(0, d[1],d[3])
        mx[z2d] <- v1perslice[,k]
        V1[,sl,,k] <- mx
        mx <- matrix(0, d[1],d[3])
        mx[z2d] <- v2perslice[,k]
        V2[,sl,,k] <- mx
        mx <- matrix(0, d[1],d[3])
        mx[z2d] <- v3perslice[,k]
        V3[,sl,,k] <- mx }, # coronal
        { mx <- matrix(0, d[1],d[2])
        mx[z2d] <- v1perslice[,k]
        V1[,,sl,k] <- mx
        mx <- matrix(0, d[1],d[2])
        mx[z2d] <- v2perslice[,k]
        V2[,,sl,k] <- mx 
        mx <- matrix(0, d[1],d[2])
        mx[z2d] <- v3perslice[,k]
        V3[,,sl,k] <- mx } ) # axial
      }
      ## gfas volume
      fsave <- paste(savedir,"/vol",sl,".RData",sep="")
      switch(kv,
      { mx <- matrix(0, d[2],d[3])
      mx[z2d] <- gfas
      volgfa[sl,,] <- mx 
      res <- list(kv=kv, gfa=volgfa[sl,,],  v1=V1[sl,,,], v2=V2[sl,,,], v3=V3[sl,,,],
        file=fsave) }, # sagittal
      { mx <- matrix(0, d[1],d[3])
      mx[z2d] <- gfas
      volgfa[,sl,] <- mx
      res <- list(kv=kv, gfa=volgfa[,sl,],  v1=V1[,sl,,], v2=V2[,sl,,], v3=V3[,sl,,],
        file=fsave) }, # coronal
      { mx <- matrix(0, d[1],d[2])
      mx[z2d] <- gfas
      volgfa[,,sl] <- mx
      res <- list(kv=kv, gfa=volgfa[,,sl],  v1=V1[,,sl,], v2=V2[,,sl,], v2=V3[,,sl,],
        file=fsave) } ) # axial
      ##
      save(res, file=fsave)
      cat("wrote", fsave,"\n")
    } else {
      fsave <- paste(savedir,"/vol",sl,".RData",sep="")
      load(fsave)
      cat("loaded", fsave, "\n")
      switch(res$kv,
      { V1[sl,,,] <- res$v1
        V2[sl,,,] <- res$v2
        V3[sl,,,] <- res$v3
        volgfa[sl,,] <- res$gfa },
      { V1[,sl,,] <- res$v1
        V2[,sl,,] <- res$v2
        V3[,sl,,] <- res$v3
        volgfa[,sl,] <- res$gfa },
      { V1[,,sl,] <- res$v1
        V2[,,sl,] <- res$v2
        V3[,,sl,] <- res$v3
        volgfa[,,sl] <- res$gfa } )
    }
    print(proc.time() - ptm)
  }
  f <- paste(savedir,"/data_gfa",sep="")
  writeNIfTI(volgfa, filename=f, verbose=TRUE)
  cat("wrote",f,"\n")
  f <- paste(savedir,"/data_V1",sep="")
  writeNIfTI(V1, filename=f, verbose=TRUE)
  cat("wrote",f,"\n")
  f <- paste(savedir,"/data_V2",sep="")
  writeNIfTI(V2, filename=f, verbose=TRUE)
  cat("wrote",f,"\n")
  f <- paste(savedir,"/data_V3",sep="")
  writeNIfTI(V3, filename=f, verbose=TRUE)
  cat("wrote",f,"\n")
  ##
  V123 <- array(0, dim=c(d, 9)) ## V1 direction
	V123[,,,1:3] <- V1
	V123[,,,4:6] <- V2
	V123[,,,7:9] <- V3
  f <- paste(savedir,"/data_V123",sep="")
  writeNIfTI(V123, filename=f, verbose=TRUE)
  cat("wrote",f,"\n")
}