R/spatialnoise.R

In NeuroStat/neuRosim: Functions to Generate fMRI Data Including Activated Data, Noise Data and Resting State Data

spatialnoise <-
function(dim, sigma, nscan, method=c("corr", "gammaRF", "gaussRF"), type=c("gaussian","rician"), rho=0.75, FWHM=4, gamma.shape=6, gamma.rate=1, vee=1, template, verbose=TRUE, voxdim){

        if(length(dim)>3){
                stop("Image space with more than three dimensions is not supported.")
        }
	if(length(dim)==1){
		stop("Spatially correlated noise structures are not defined for vectors.")
	}
	if(missing(method)){
		method <- "corr"
	}
	if(missing(type)&&(method=="corr")){
		type <- "gaussian"
	}

	if(method=="corr"){
		if(length(dim)==2){
			noise <- array(0, dim=c(dim,nscan))
			for(z in 1:nscan){
				if(type=="gaussian"){
				  start <- array(rnorm(prod(dim), 0, sigma), dim=dim)
				}else{
				  start <- array(rrice(prod(dim), vee, sigma), dim=dim)
				}
				noise.scan <- array(0, dim=dim)
				noise.scan[1,1] <- start[1,1]
				for(i in 2:dim[1]){
					noise.scan[i,1] <-rho*noise.scan[(i-1),1] + sqrt(1-rho^2)*start[i,1]
				}
				for(j in 2:dim[2]){
					noise.scan[,j] <- rho*noise.scan[,(j-1)] + sqrt(1-rho^2)*start[,j]
				}
				for(i in 2:dim[1]){
					noise.scan[i,2:dim[2]] <- rho*noise.scan[(i-1),2:dim[2]] + sqrt(1-rho^2)*noise.scan[i,2:dim[2]]
				}
				noise[,,z] <- noise.scan
			}
		} else {
			noise <- array(0, dim=c(dim,nscan))
			for(z in 1:nscan){
				if(type=="gaussian"){
				  start <- array(rnorm(prod(dim), 0, sigma), dim=dim)
				}else{
				  start <- array(rrice(prod(dim), vee, sigma), dim=dim)
				}
				noise.scan <- array(0, dim=dim)
				noise.scan[1,1,1] <- start[1,1,1]
				for(i in 2:dim[1]){
					noise.scan[i,1,1] <- rho*noise.scan[(i-1),1,1] + sqrt(1-rho^2)*start[i,1,1]
				}
				for(j in 2:dim[2]){
					noise.scan[,j,1] <- rho*noise.scan[,(j-1),1] + sqrt(1-rho^2)*start[,j,1]
				}
				for(k in 2:dim[3]){
					noise.scan[,,k] <- rho*noise.scan[,,(k-1)] + sqrt(1-rho^2)*start[,,k]
				}
				for(i in 2:dim[1]){
					noise.scan[i,2:dim[2],2:dim[3]] <- rho*noise.scan[(i-1),2:dim[2],2:dim[3]] + sqrt(1-rho^2)*noise.scan[i,2:dim[2],2:dim[3]]
				}
				noise[,,,z] <- noise.scan
			}
		}
       	        if(!missing(template)){
 			if(length(dim(template))>3){
				stop("Template should be a 2D or 3D array.")
			}
	                template.time <- array(rep(template,nscan), dim=c(dim,nscan))
       		        ix <- which(template.time!=0)
                	noise[-ix] <- 0
        	}
	}
	if(method=="gaussRF" | method=="gammaRF"){
		#require(AnalyzeFMRI, quietly=TRUE)
        ksd = FWHM/(sqrt(8*log(2))) # in mm
		s <- diag(ksd^2, 3)
        if(length(dim)==2){
                dim.RF <- c(dim,1)
        } else {
                dim.RF <- dim
        }
        if (missing(voxdim)){
            cat("spatial noise:  missing voxdim\n")
			voxdim <- rep(1, length(dim.RF))
        }
		if(!missing(template)){
			if(length(dim(template))>3){
				stop("Template should be a 2D or 3D array.")
			}
 			m <- array(ifelse(template!=0, 1, 0), dim=dim.RF)
		} else {
			m <- array(1, dim=dim.RF)
		}

                # make the kernel at least 5sd wide
                ksize = ceiling(5.0*ksd/min(voxdim)) # (mm x vox/mm = vox; kernel is isotropic; width governed by highest resolution / smallest voxel dimension)
                if (ksize%%2 == 0){
                	ksize = ksize + 1
                }

		noise <- array(0, dim=c(dim.RF,nscan))
	        if(method=="gaussRF"){
			for(z in 1:nscan){
				noise[,,,z] <- array(c(Sim.3D.GRF(d=dim.RF,voxdim=voxdim,sigma=s,ksize=ksize,mask=m,type="field")$mat), dim=dim.RF)
		   	}
                	scale = sigma
                } else {
                	for(z in 1:nscan){
                        	n <- Sim.3D.GRF(d=dim.RF,voxdim=voxdim,sigma=s,ksize=ksize,mask=m,type="field")$mat
   				gamma.n <- qgamma(pnorm(c(n)), shape=gamma.shape, rate=gamma.rate)
   				noise[,,,z] <- array(gamma.n, dim=dim.RF)
                   	}
                   	scale = sigma*gamma.rate/sqrt(gamma.shape)
                }
		noise <- array(c(noise), dim=c(dim,nscan))
                noise = noise * scale # scale std. dev. to sigma
	}

  # Noise is already rescaled
	return(noise)
}