.packageName <- "AnalyzeFMRI"

f.read.analyze.header <- function(file){
  #This function reads in the information from an ANALYZE format .hdr file
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")

    if(file.exists(file.img) == FALSE) return(paste(file.img, "not found"))
    if(file.exists(file.hdr) == FALSE) return(paste(file.hdr, "not found"))

#Detect whether the data is big or little endian. The first part of a .hdr file is the size of the file which is always a C int (i.e. 4 bytes) and always has value 348. Therefore trying to read it in assuming little-endian will tell you if that is the correct mode

    swap <- 0

    if(.C("swaptest_wrap_JM",
          ans = integer(1),
          file.hdr,
          PACKAGE="AnalyzeFMRI")$ans != 348)
        swap <- 1


# A C function is used to read in all the components of the .hdr file
    a<-.C("read_analyze_header_wrap_JM",
          file.hdr,
          as.integer(swap),
          integer(1),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 18), sep = "", collapse = ""),
          integer(1),
          integer(1),
          paste(rep(" ", 1), sep = "", collapse = ""),
          paste(rep(" ", 1), sep = "", collapse = ""),
          integer(8),
          paste(rep(" ", 4), sep = "", collapse = ""),
          paste(rep(" ", 8), sep = "", collapse = ""),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          single(8),
          single(1),
          single(1),
          single(1),
          single(1),
          single(1),
          single(1),
          single(1),
          single(1),
          integer(1),
          integer(1),
          paste(rep(" ", 80), sep = "", collapse = ""),
          paste(rep(" ", 24), sep = "", collapse = ""),
          paste(rep(" ", 1), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ", 10), sep = "", collapse = ""),
          paste(rep(" ",3 ), sep = "",collapse = ""),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          integer(1),
          PACKAGE="AnalyzeFMRI")

#A list (called L) is created containing a selection of the useful components of the .hdr file

    L <- list()
    L$swap <- a[[2]]
    L$file.name <- file.img
    L$dim <- a[[10]]
    L$vox.units <- a[[11]]
    L$cal.units <- a[[12]]
    L$datatype <- a[[14]]
    if(L$datatype == 0 ) L$data.type <- "unknown"
    if(L$datatype == 1) L$data.type <- "binary"
    if(L$datatype == 2) L$data.type <- "unsigned char"
    if(L$datatype == 4) L$data.type <- "signed short"
    if(L$datatype == 8) L$data.type <- "signed int"
    if(L$datatype == 16) L$data.type <- "float"
    if(L$datatype == 32) L$data.type <- "complex"
    if(L$datatype == 64) L$data.type <- "double precision"
    if(L$datatype == 128) L$data.type <- "rgb data"
    if(L$datatype == 255) L$data.type <- "all"
    L$bitpix <- a[[15]]
    L$pixdim <- a[[17]]
    L$glmax <- a[[26]]
    L$glmin <- a[[27]]

    return(L)}


f.analyze.file.summary <- function(file){
#This function prints out a concise summary of the contents of a .img/.hdr image pair

    file.name <- substring(file ,1 , nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
    cat("\n")
    cat("       File name:", file.img, "\n")
    cat("  Data Dimension:", paste(hdr$dim[1], "-D", sep = ""), "\n")
    cat("     X dimension:", hdr$dim[2], "\n")
    cat("     Y dimension:", hdr$dim[3], "\n")
    cat("     Z dimension:", hdr$dim[4], "\n")
    cat("  Time dimension:", hdr$dim[5], "time points", "\n")
    cat("Voxel dimensions:", paste(hdr$pixdim[2], hdr$vox.units, "x",
                                   hdr$pixdim[3], hdr$vox.units, "x",
                                   hdr$pixdim[4], hdr$vox.units), "\n")
    cat("       Data type:", hdr$data.type, paste("(", hdr$bitpix, " bits per voxel)", sep = ""), "\n")
}


f.basic.hdr.list.create <- function(mat, file.hdr){

#creates a basic list that can be used to write a .hdr file

    dim <- dim(mat)
    dim <- c(length(dim), dim, rep(0, 7 - length(dim)))

    l <- list(file = file.hdr,
              size.of.header = 348,
              data.type = paste(rep(" ", 10), sep = "", collapse = ""),
              db.name = paste(rep(" ", 18), sep = "", collapse = ""),
              extents = 0,
              session.error = 0,
              regular = character(1),
              hkey.un0 = character(1),
              dim = as.integer(dim),
              vox.units = "mm",
              cal.units = "voxels",
              unused1 = 0,
              datatype = 0,
              bitpix = 0,
              dim.un0 = 0,
              pixdim = single(8),
              vox.offset = single(1),
              funused1 = single(1),
              funused2 = single(1),
              funused3 = single(1),
              cal.max = single(1),
              cal.min = single(1),
              compressed = single(1),
              verified = single(1),
              glmax = 0,
              glmin = 0,
              descrip = paste(rep(" ", 80), sep = "", collapse = ""),
              aux.file = paste(rep(" ", 24), sep = "", collapse = ""),
              orient = paste(rep(" ", 1), sep = "", collapse = ""),
              originator = paste(rep(" ", 10), sep = "", collapse = ""),
              generated = paste(rep(" ", 10), sep = "", collapse = ""),
              scannum = paste(rep(" ", 10), sep = "", collapse = ""),
              patient.id = paste(rep(" ", 10), sep = "", collapse = ""),
              exp.date = paste(rep(" ", 10), sep = "", collapse = ""),
              exp.time = paste(rep(" ", 10), sep = "", collapse = ""),
              hist.un0 = paste(rep(" ", 3), sep = "", collapse = ""),
              views = integer(1),
              vols.added = integer(1),
              start.field = integer(1),
              field.skip = integer(1),
              omax = integer(1),
              omin = integer(1),
              smax = integer(1),
              smin = integer(1) )
    return(l)
}

f.read.analyze.slice <- function(file, slice, tpt){
  #Reads in a .img file into an array
    file.name <- substring(file ,1 ,nchar(file) - 4)
    file.hdr <- paste(file.name ,".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
    #f.analyze.file.summary(file)

  #num.dim<-hdr$dim[1]
    dim <- hdr$dim[2:3]

    num.data.pts <- dim[1] * dim[2]
    if(tpt < 1 || tpt > hdr$dim[5]) stop("tpt is not in range")
    if(slice < 1 || slice > hdr$dim[4]) stop("slice is not in range")

    offset <- (tpt - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
    if(hdr$datatype == 2){

        vol <- .C("readchar_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 1),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
#this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 4){

        vol <- .C("read2byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 2),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
#this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 8){
        vol <- .C("read4byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 4),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 16){
        vol <- .C("readfloat_v1_JM",
                  mat = single(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 4),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim=dim)
#this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 64){
        vol <- .C("readdouble_v1_JM",
                  mat = numeric(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 8),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }

    if(hdr$datatype == 0 || hdr$datatype == 1 || hdr$datatype == 32 || hdr$datatype == 128 || hdr$datatype == 255) print(paste("The format", hdr$data.type, "is not supported yet. Please contact me if you want me to extend the functions to do this (marchini@stats.ox.ac.uk)"), quote=FALSE)

    return(vol)}


f.read.analyze.tpt <- function(file, tpt){
#Reads in one timepoint of a .img file into an array
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
  #f.analyze.file.summary(file)

  #num.dim <- hdr$dim[1]
    dim <- hdr$dim[2:4]

    num.data.pts <- dim[1] * dim[2] * dim[3]
    if(tpt < 1 || tpt > hdr$dim[5]) stop("tpt is not in range")

    offset <- (tpt - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4]

    if(hdr$datatype == 2){

        vol <- .C("readchar_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 1),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
                                        #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 4){
        vol <- .C("read2byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 2),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 8){
        vol <- .C("read4byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 4),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 16){
        vol <- .C("readfloat_v1_JM",
                  mat = single(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 4),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 64){
        vol <- .C("readdouble_v1_JM",
                  mat = numeric(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(offset * 8),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }

    if(hdr$datatype == 0 || hdr$datatype == 1 || hdr$datatype == 32 || hdr$datatype == 128 || hdr$datatype == 255) print(paste("The format", hdr$data.type, "is not supported yet. Please contact me if you want me to extend the functions to do this (marchini@stats.ox.ac.uk)"), quote = FALSE)

    return(vol)}


f.read.analyze.slice.at.all.timepoints <- function(file, slice){
  #Reads in a slice of a .img file at all time points into an array
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
  #f.analyze.file.summary(file)

  #num.dim <- hdr$dim[1]
    dim <- hdr$dim[2:3]

    num.data.pts <- dim[1] * dim[2]
    if(slice < 1 || slice > hdr$dim[4]) stop("slice is not in range")

    vl <- array(0, dim = hdr$dim[c(2, 3, 5)])

    if(hdr$datatype == 2){
        for(i in 1:hdr$dim[5]){
            offset <- (i - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
            vol <- .C("readchar_v1_JM",
                      mat = integer(num.data.pts),
                      file.img,
                      as.integer(hdr$swap),
                      as.integer(num.data.pts),
                      as.integer(offset * 1),
                      as.integer(1), PACKAGE="AnalyzeFMRI")
            vol <- array(vol$mat, dim = dim)
            vl[, , i] <- vol
        }
    }
    if(hdr$datatype == 4){
        for(i in 1:hdr$dim[5]){
            offset <- (i - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
            vol <- .C("read2byte_v1_JM",
                      mat = integer(num.data.pts),
                      file.img,
                      as.integer(hdr$swap),
                      as.integer(num.data.pts),
                      as.integer(offset * 2),
                      as.integer(1), PACKAGE="AnalyzeFMRI")
            vol <- array(vol$mat, dim = dim)
            vl[, , i] <- vol
        }
    }

    if(hdr$datatype == 8){
        for(i in 1:hdr$dim[5]){
            offset <- (i - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
            vol <- .C("read4byte_v1_JM",
                      mat = integer(num.data.pts),
                      file.img,
                      as.integer(hdr$swap),
                      as.integer(num.data.pts),
                      as.integer(offset * 4),
                      as.integer(1), PACKAGE="AnalyzeFMRI")
            vol <- array(vol$mat, dim = dim)
            vl[, , i] <- vol
        }
    }

    if(hdr$datatype == 16){
        for(i in 1:hdr$dim[5]){
            offset <- (i - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
            vol <- .C("readfloat_v1_JM",
                      mat = single(num.data.pts),
                      file.img,
                      as.integer(hdr$swap),
                      as.integer(num.data.pts),
                      as.integer(offset * 4),
                      as.integer(1), PACKAGE="AnalyzeFMRI")
            vol <- array(vol$mat, dim = dim)
        }
    }

    if(hdr$datatype == 64){
        for(i in 1:hdr$dim[5]){
            offset <- (i - 1) * hdr$dim[2] * hdr$dim[3] * hdr$dim[4] + (slice - 1) * hdr$dim[2] * hdr$dim[3]
            vol <- .C("readdouble_v1_JM",
                      mat = numeric(num.data.pts),
                      file.img,
                      as.integer(hdr$swap),
                      as.integer(num.data.pts),
                      as.integer(offset * 8),
                      as.integer(1), PACKAGE="AnalyzeFMRI")
            vol <- array(vol$mat, dim = dim)
        }}

    if(hdr$datatype == 0 || hdr$datatype == 1 || hdr$datatype == 32 || hdr$datatype == 128 || hdr$datatype == 255) print(paste("The format", hdr$data.type, "is not supported yet. Please contact me if you want me to extend the functions to do this (marchini@stats.ox.ac.uk)"), quote = FALSE)

    return(vl)}

f.read.analyze.ts <- function(file, x, y, z){
  #Reads in a .img file into an array
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
  #f.analyze.file.summary(file)

    if(x < 1 || x > hdr$dim[2]) stop("x is not in range")
    if(y < 1 || y > hdr$dim[3]) stop("y is not in range")
    if(z < 1 || z > hdr$dim[4]) stop("z is not in range")

    offset.start <- (z - 1) * hdr$dim[2] * hdr$dim[3] + (y - 1) * hdr$dim[2] + (x - 1)
    offset.add <- hdr$dim[2] * hdr$dim[3] * hdr$dim[4]

    vol <- 1:hdr$dim[5]

    if(hdr$datatype == 2){
        for(i in 1:hdr$dim[5]){

            v <- .C("readchar_v1_JM",
                    mat = integer(1),
                    file.img,
                    as.integer(hdr$swap),
                    as.integer(1),
                    as.integer(offset.start * 1 + 1 * (i - 1) * offset.add),
                    as.integer(1), PACKAGE="AnalyzeFMRI")
            vol[i] <- v$mat}
    }

    if(hdr$datatype  == 4){
        for(i in 1:hdr$dim[5]){

            v <- .C("read2byte_v1_JM",
                    mat = integer(1),
                    file.img,
                    as.integer(hdr$swap),
                    as.integer(1),
                    as.integer(offset.start * 2 + 2 * (i - 1) * offset.add),
                    as.integer(1), PACKAGE="AnalyzeFMRI")
            vol[i] <- v$mat}
    }

    if(hdr$datatype  == 8){
        for(i in 1:hdr$dim[5]){

            v <- .C("read4byte_v1_JM",
                    mat = integer(1),
                    file.img,
                    as.integer(hdr$swap),
                    as.integer(1),
                    as.integer(offset.start * 4 + 4 * (i - 1) * offset.add),
                    as.integer(1), PACKAGE="AnalyzeFMRI")
            vol[i] <- v$mat}
    }

    if(hdr$datatype == 16){
        for(i in 1:hdr$dim[5]){

            v <- .C("readfloat_v1_JM",
                    mat = integer(1),
                    file.img,
                    as.integer(hdr$swap),
                    as.integer(1),
                    as.integer(offset.start * 4 + 4 * (i - 1) * offset.add),
                    as.integer(1), PACKAGE="AnalyzeFMRI")
            vol[i] <- v$mat}
    }

    if(hdr$datatype == 64){
        for(i in 1:hdr$dim[5]){

            v <- .C("readdouble_v1_JM",
                    mat = integer(1),
                    file.img,
                    as.integer(hdr$swap),
                    as.integer(1),
                    as.integer(offset.start * 8 + 8 * (i - 1) * offset.add),
                    as.integer(1), PACKAGE="AnalyzeFMRI")
            vol[i] <- v$mat}
    }

    if(hdr$datatype == 0 || hdr$datatype == 1 || hdr$datatype == 32 || hdr$datatype == 128 || hdr$datatype == 255) print(paste("The format", hdr$data.type, "is not supported yet. Please contact me if you want me to extend the functions to do this (marchini@stats.ox.ac.uk)"), quote = FALSE)

    return(vol)}



f.read.analyze.volume <- function(file){
  #Reads in a .img file into an array
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
  #f.analyze.file.summary(file)
    num.dim <- hdr$dim[1]
    dim <- hdr$dim[1:num.dim + 1]

    num.data.pts <- 1
    for(i in 1:num.dim){num.data.pts <- num.data.pts * dim[i]}


    if(hdr$datatype == 2){
        vol <- .C("readchar_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(0),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 4){
        vol <- .C("read2byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(0),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
        vol <- array(vol$mat, dim = dim)
#this works because array fills itself with the left most subscript moving fastest
    }
    if(hdr$datatype == 8){
        vol <- .C("read4byte_v1_JM",
                  mat = integer(num.data.pts),
                  file.img,
                  as.integer(hdr$swap),
                  as.integer(num.data.pts),
                  as.integer(0),
                  as.integer(1), PACKAGE="AnalyzeFMRI")
    vol <- array(vol$mat, dim = dim)
#this works because array fills itself with the left most subscript moving fastest
}
if(hdr$datatype == 16){
    vol <- .C("readfloat_v1_JM",
              mat = single(num.data.pts),
              file.img,
              as.integer(hdr$swap),
              as.integer(num.data.pts),
              as.integer(0),
              as.integer(1), PACKAGE="AnalyzeFMRI")
    vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
}
if(hdr$datatype == 64){
    vol <- .C("readdouble_v1_JM",
              mat = numeric(num.data.pts),
              file.img,
              as.integer(hdr$swap),
              as.integer(num.data.pts),
              as.integer(0),
              as.integer(1), PACKAGE="AnalyzeFMRI")
    vol <- array(vol$mat, dim = dim)
    #this works because array fills itself with the left most subscript moving fastest
}

if(hdr$datatype == 0 || hdr$datatype == 1 || hdr$datatype == 32 || hdr$datatype == 128 || hdr$datatype == 255) print(paste("The format", hdr$data.type, "is not supported yet. Please contact me if you want me to extend the functions to do this (marchini@stats.ox.ac.uk)"), quote = FALSE)

return(vol)}



f.spectral.summary <- function(file, mask.file, ret.flag = FALSE)
{
  #for an analyze .img file the periodogram of the time series are divided by a flat spectral estimate using the median periodogram ordinate. The resulting values are then combined within each Fourier frequency and quantiles are plotted against freequency. This provides a fast look at a fMRI dataset to identify any artefacts that reside at single frequencies.

########################
#get info about dataset
########################
    file.name <- substring(file, 1, nchar(file) - 4)
    file.hdr <- paste(file.name, ".hdr", sep = "")
    file.img <- paste(file.name, ".img", sep = "")
    hdr <- f.read.analyze.header(file.hdr)
    nsl <- hdr$dim[4]
    nim <- hdr$dim[5]
    pxdim <- hdr$pixdim[2:4]

#####################
#read in/create mask
#####################

    f.mask.create <- function(dat, pct = .1, slices = c(0)) {
  #function that creats a mask for an fMRI dataset by thresholding the mean of the pixel time series at a percentage point of the maximum intensity of the dataset
        file.name <- substring(dat$file, 1, nchar(dat$file) - 4)
        file.hdr <- paste(file.name, ".hdr", sep = "")
        hdr <- f.read.analyze.header(file.hdr)
        nsl <- hdr$dim[4]
        xc <- hdr$dim[2]
        yc <- hdr$dim[3]
        if(slices[1] == 0){slices <- seq(1, nsl)}
        mask <- array(0, dim = c(xc, yc, length(slices)))

        max.int <- 0
        for(k in 1:length(slices)){
            slice <- f.read.analyze.slice.at.all.timepoints(dat$file, slices[k])
            if(max(slice)>max.int){max.int <- max(slice)}
        }

        for(k in 1:length(slices)){
            slice <- f.read.analyze.slice.at.all.timepoints(dat$file, slices[k])

            for(i in 1:(xc * yc)){
                a <- (i - 1) %/% xc + 1
                b <- i - (a - 1) * xc
                mask[b, a, k] <- mean(slice[b, a, ])

                if(mask[b, a, k] >= (pct * max.int)){mask[b, a, k] <- 1}
                else{mask[b, a, k] <- 0}
            }
        }

        return(mask)

    }

    if(mask.file!=FALSE){mask <- f.read.analyze.volume(mask.file)}
    else{
        dat <- list(file = file, mask.file = mask.file)
        mask <- f.mask.create(dat = dat)}
    dim(mask) <- hdr$dim[2:4]

###############
#set constants
###############

    n <- floor(nim / 2) + 1


##########################
#initialise storage arrays
##########################

    res <- array(NA, dim = c(dim(mask), n))

#####################
#main evaluation loop
#####################
    cat("Processing slices...")
    for(l in 1:nsl){
        cat(" [", l, "]", sep = "")

        slice <- f.read.analyze.slice.at.all.timepoints(file, l)


        for(i in 1:(dim(slice)[1] * dim(slice)[2])){
            a <- (i - 1) %/% dim(slice)[1] + 1
            b <- i - (a - 1) * dim(slice)[1]
            if(mask[b, a, l] == 1){
                t <- Mod(fft(slice[b, a, ]) / sqrt(2 * pi * nim))[1:n]
                s <- median(t)
                res[b, a, l, ] <- t / s
            }
        }
    }
    cat("\n")

    b <- apply(res, 4, FUN = quantile, probs = seq(.5, 1, .05), na.rm = TRUE)
    plot(c(0, n - 1), c(0, 30), type = "n", xlab = "", ylab = "", axes = FALSE)
    axis(1, at = seq(0, n, 5))
    axis(2, at = seq(0, 30, 5))
    for(i in 0:(n - 1)){
        points(rep(i, 11), b[, i + 1])}
    if(ret.flag == TRUE)  return(b)
}


f.write.list.to.hdr <- function(L, file){

# Writes a list to a .hdr file
    a <- .C("write_analyze_header_wrap_JM",
            file,
            as.integer(L$size.of.header),
            as.character(L$data.type),
            as.character(L$db.name),
            as.integer(L$extents),
            as.integer(L$session.error),
            as.character(L$regular),
            as.character(L$hkey.un0),
            as.integer(L$dim),
            as.character(L$vox.units),
            as.character(L$cal.units),
            as.integer(L$unused1),
            as.integer(L$datatype),
            as.integer(L$bitpix),
            as.integer(L$dim.un0),
            as.single(L$pixdim),
            as.single(L$vox.offset),
            as.single(L$funused1),
            as.single(L$funused2),
            as.single(L$funused3),
            as.single(L$cal.max),
            as.single(L$cal.min),
            as.single(L$compressed),
            as.single(L$verified),
            as.integer(L$glmax),
            as.integer(L$glmin),
            as.character(L$descrip),
            as.character(L$aux.file),
            as.character(L$orient),
            as.character(L$originator),
            as.character(L$generated),
            as.character(L$scannum),
            as.character(L$patient.id),
            as.character(L$exp.date),
            as.character(L$exp.time),
            as.character(L$hist.un0),
            as.integer(L$views),
            as.integer(L$vols.added),
            as.integer(L$start.field),
            as.integer(L$field.skip),
            as.integer(L$omax),
            as.integer(L$omin),
            as.integer(L$smax),
            as.integer(L$smin),
            PACKAGE="AnalyzeFMRI")
}


f.write.analyze <- function(mat, file, size = "float", pixdim = c(4, 4, 6), vox.units = "mm", cal.units = "pixels"){

  #Creates a .img and .hdr pair of files froma given array

    if(max(mat) == "NA") return("NA values in array not allowed. Files not written.")

    file.img <- paste(file, ".img", sep = "")
    file.hdr <- paste(file, ".hdr", sep = "")

    L <- f.basic.hdr.list.create(mat, file.hdr)

    if(size == "float"){
        L$datatype <- 16
        L$bitpix <- 32
        L$vox.units <- vox.units
        L$cal.units <- cal.units
        L$pixdim <- c(4, pixdim, 0, 0, 0, 0)
        L$data.type <- "float"
        L$glmax <- as.integer(floor(max(mat)))
        L$glmin <- as.integer(floor(min(mat)))
        f.write.array.to.img.float(mat, file.img)
    }
    if(size == "int"){
        if(max(mat)>32767 || min(mat) < ( -32768)) return("Values are outside integer range. Files not written.")
        L$datatype <- 4
        L$bitpix <- 16
        L$vox.units <- vox.units
        L$cal.units <- cal.units
        L$pixdim <- c(4, pixdim, 0, 0, 0, 0)
        L$data.type <- "signed short"
        L$glmax <- as.integer(floor(max(mat)))
        L$glmin <- as.integer(floor(min(mat)))
        f.write.array.to.img.2bytes(mat, file.img)
    }
    if(size == "char"){
        if(max(mat)>255 || min(mat) < 0) return("Values are outside integer range. Files not written.")
        L$datatype <- 2
        L$bitpix <- 8
        L$vox.units <- vox.units
        L$cal.units <- cal.units
        L$pixdim <- c(4, pixdim, 0, 0, 0, 0)
        L$data.type <- "unsigned char"
        L$glmax <- as.integer(floor(max(mat)))
        L$glmin <- as.integer(floor(min(mat)))
        f.write.array.to.img.8bit(mat, file.img)
    }

    f.write.list.to.hdr(L, file.hdr)
}


f.write.array.to.img.2bytes <- function(mat, file){
  #writes an array into a .img file of 2 byte integers

    dm <- dim(mat)
    dm.ln <- length(dm)
    num.data.pts <- length(mat)

    .C("write2byte_JM",
       as.integer(mat),
       file,
       as.integer(num.data.pts), PACKAGE="AnalyzeFMRI")

}

f.write.array.to.img.8bit <- function(mat, file){
#writes an array into a .img file of 8 bit (1 byte) integers

    dm <- dim(mat)
    dm.ln <- length(dm)
    num.data.pts <- length(mat)

    .C("write8bit_JM",
       as.integer(mat),
       file,
       as.integer(num.data.pts), PACKAGE="AnalyzeFMRI")

}



f.write.array.to.img.float <- function(mat, file){
  #writes an array into a .img file of 4 byte flotas

    dm <- dim(mat)
    dm.ln <- length(dm)
    num.data.pts <- length(mat)

    .C("writefloat_JM",
       as.single(mat),
       file,
       as.integer(num.data.pts), PACKAGE="AnalyzeFMRI")
}

f.analyzeFMRI.gui <- function(){

  #starts GUI that allows user to explore an fMRI dataset stored in an ANALYZEfile

    path <- .path.package(package = "AnalyzeFMRI")
    path.gui <- paste(path, "AnalyzeFMRI.gui.R", sep = .Platform$file.sep)
    source(path.gui)}


f.ica.fmri <- function(file.name, n.comp, norm.col = TRUE, fun = "logcosh", maxit = 1000, alg.type = "parallel", alpha = 1, tol = 0.0001, mask.file.name = NULL, slices = NULL){

  #function for performing Spatial ICA on an fMRI dataset
  #The function avoids reading the dataset into R to minimise the memory used

    hdr <- f.read.analyze.header(file.name)

    if(length(slices) == 0) slices  <-  2:(hdr$dim[4] - 1)
    if(slices == "all") slices  <-  1:(hdr$dim[4])
    if(any(slices < 1 || slices>hdr$dim[4])) {
        return("some of selected slices out of allowable range")}

    ns <- hdr$dim[2] * hdr$dim[3] * hdr$dim[4] * n.comp
    na <- hdr$dim[5] * n.comp

    mask.flag <- 1
    if(length(mask.file.name) == 0){
        mask.flag <- 0
        mask.file.name <- ""}


    col.flag <- 1
    if(norm.col!=TRUE) col.flag <- 0

    fun.flag <- 1
    if(fun == "exp") fun.flag <- 2

    def.flag <- 0
    if(alg.type == "deflation") def.flag <- 1

    W <- matrix(rnorm(n.comp * n.comp), n.comp, n.comp)

    a <- .C("ica_fmri_JM",
            as.character(file.name),
            as.single(t(W)),
            as.integer(n.comp),
            as.integer(1),
            as.integer(col.flag),
            as.integer(fun.flag),
            as.integer(maxit),
            as.integer(def.flag),
            as.single(alpha),
            as.single(tol),
            as.integer(mask.flag),
            as.character(mask.file.name),
            as.integer(slices),
            as.integer(length(slices)),
            S = single(ns),
            A = single(na),
            PACKAGE="AnalyzeFMRI")

    S <- array(a$S, dim = c(hdr$dim[2], hdr$dim[3], hdr$dim[4], n.comp))

    A <- matrix(a$A, hdr$dim[5], n.comp, byrow = TRUE)

    return(list(A = A, S = S, file = file.name, mask = mask.file.name))
}

f.plot.ica.fmri.jpg <- function(ica.obj, file = "./ica", cols = heat.colors(100), width = 700,  height = 700){

    for(i in 1:dim(ica.obj$S)[4]){

        jpeg(file = paste(file, ".comp.", i, ".jpeg", sep = ""), width = width, height = height)
        f.plot.ica.fmri(ica.obj,  i,  cols = cols)
        dev.off()
    }
    return()
}


f.plot.ica.fmri <- function (obj.ica,  comp,  cols = heat.colors(100))
{
    r  <-  range(obj.ica$S[,  ,  ,  comp],  na.rm = TRUE)
    tmp  <-  1000 * (obj.ica$S[,  ,  ,  comp]  == 0) + (obj.ica$S[,  ,  ,  comp]) * (obj.ica$S[,  ,  ,  comp] != 0)
    ncomp  <-  dim(obj.ica$S)[4]
    nsl  <-  dim(obj.ica$S)[3]
    t  <-  nrow(obj.ica$A)
    im  <-  floor(sqrt(nsl + 3)) + 1
    par(mfrow = c(im,  im),  mar = c(.5,  .5,  .5,  .5))
    plot(c(0, 1), c(0, 1), typ = "n", axes = FALSE, xlab = "", ylab = "")

    text(0, .9, "Spatial ICA", pos = 4, cex = 1.5)
    text(0, .75, paste("Component ", comp, sep = ""), pos = 4, cex = 1.5)
    l <- floor(nchar(obj.ica$file) / 20) + 1
    text(0, .6, paste("file: ", substring(obj.ica$file, 1, 20), sep = ""), pos = 4)
    for(i in 2:l){
        text(0, .6 - .07 * (i - 1), paste("       ", substring(obj.ica$file, 20 * (i - 1) + 1, 20 * i), sep = ""), pos = 4)
    }
    text(0, .6 - .07 * (l + 1), paste("Date: ", date(), sep = ""), pos = 4)

    for (i in 1:nsl) {
        image(tmp[,  ,  i],  zlim = r,  axes = FALSE,  col = cols)
        text(.5, .98, paste("slice", i), pos = 1)
        box()
    }
    r <- range(obj.ica$A[,  comp])
    plot(obj.ica$A[,  comp],  typ = "l",  axes = FALSE, ylim = r * 1.5)
    text(length(obj.ica$A[,  comp]) / 2, 1.5 * r[2], "Time course", pos = 1)
    box()
    s  <-  fft(obj.ica$A[,  comp]) / sqrt(2 * pi * t)
    s  <-  Mod(s[2:(floor(t / 2) + 1)])^2
    r <- c(0, max(s))
    plot(s,  axes = FALSE, typ = "l", ylim = r * 1.5)
    text(length(s) / 2, 1.5 * r[2], "Periodogram", pos = 1)
    box()
    par(mfrow = c(1,  1),  mar = c(5,  4,  4,  2))
}

f.ica.fmri.gui <- function(){

  #starts GUI that allows user apply Spatial ICA to an fMRI dataset

    path <- .path.package(package = "AnalyzeFMRI")
    path.gui <- paste(path, "ICA.gui.R", sep = .Platform$file.sep)
    source(path.gui)}




#functions to apply gaussian spatial smoothing to an array

    

GaussSmoothKernel<-function(voxdim = c(1 , 1, 1), ksize = 5, sigma = diag(3, 3))
#calculates a discretized smoothing kernel in up to 3 dimensions given an arbitrary covariance matrix
#sigma is covariance matrix of the gaussian
#doesn't have to be non-singular; zero on the diagonal of sigma indicate no smoothing in that direction
  
{
    if((2 * floor(ksize / 2)) == ksize) stop(paste("ksize must be odd"))
    
    a <- array(0, dim = c(ksize, ksize, ksize))
    centre <- (ksize + 1) / 2
    
    sig.ck <- c(TRUE, TRUE, TRUE)
    for(i in 1:3){
        if(sigma[i, i] == 0){
            sigma[i, i] <- 1
            sig.ck[i] <- FALSE
        }
    }
    sig.inv <- solve(sigma)
    sig.det <- abs(det(sigma))
    
    
    
    for(i in 1:ksize) {
        for(j in 1:ksize) {
            for(k in 1:ksize) {
                x <- (c(i, j, k) - centre) * voxdim
                a[i, j, k] <- ((2 * pi)^(-3 / 2)) * exp(-.5 * (t(x) %*% sig.inv %*% x)) / sqrt(sig.det)
            }
        }
    }
    if(sig.ck[1] == FALSE) a[-centre, , ] <- 0
    if(sig.ck[2] == FALSE) a[, -centre, ] <- 0
    if(sig.ck[3] == FALSE) a[, , -centre] <- 0
    a <- a / sum(a)

    return(a)
}


GaussSmoothArray <- function(x, voxdim = c(1, 1, 1), ksize = 5, sigma = diag(3, 3), mask = NULL, var.norm= FALSE )
{
    filtmat <- GaussSmoothKernel(voxdim, ksize, sigma)
    
    if(!is.array(x)) return("x should be an array")
    if(length(dim(x)) != 3 && length(dim(x)) != 4) return("array x should be 3D or 4D")
    tmp <- FALSE
    if(length(dim(x)) == 3) {
        x <- array(x, dim = c(dim(x), 1))
        tmp <- TRUE
    }
    if(is.null(mask)) mask <- array(1, dim = dim(x)[1:3])
    
    if(var.norm){
        d <- .Fortran("gaussfilter2",
                      as.double(x),
                      as.integer(dim(x)[1]),
                      as.integer(dim(x)[2]),
                      as.integer(dim(x)[3]),
                      as.integer(dim(x)[4]),
                      as.double(filtmat),
                      as.integer(ksize),
                      as.double(mask),
                      double(length(x)),
                      PACKAGE = "AnalyzeFMRI")
        c1 <- array(d[[9]], dim = dim(x))
        if(tmp) c1 <- c1[, , , 1]
    }
        
    else
    {   d <- .Fortran("gaussfilter1",
                      as.double(x),
                      as.integer(dim(x)[1]),
                      as.integer(dim(x)[2]),
                      as.integer(dim(x)[3]),
                      as.integer(dim(x)[4]),
                      as.double(filtmat),
                      as.integer(ksize),
                      as.double(mask),
                      as.double(x),
                      PACKAGE = "AnalyzeFMRI")
        c1 <- array(d[[9]], dim = dim(x))
        if(tmp) c1 <- c1[, , , 1]
    }
    
    return(c1)
}

    
#functions to apply gaussian spatial smoothing to an array


NonLinearSmoothArray<-function(x,voxdim=c(1,1,1),radius=2,sm=3,mask=NULL)
{
    
    if(!is.array(x)) return("x should be an array")
    if(length(dim(x))!=3 && length(dim(x))!=4) return("array x should be 3D or 4D")
    if(is.null(mask)) mask <- array(1,dim=dim(x)[1:3])

    
    if(length(dim(x))==3)
    {
        d<-dim(x)
        a<-.C("non_lin_gauss_smooth",
              as.single(aperm(x,c(3,2,1))),
              as.integer(d),
              as.single(aperm(mask,c(3,2,1))),
              as.single(radius),
              as.single(sm),
              as.single(voxdim),
              single(length(x)),
              PACKAGE = "AnalyzeFMRI")
        
        a<-array(a[[7]],dim=d[3:1])
        a<-aperm(a,c(3,2,1))

        return(a)
    }
    else
    {
        d<-dim(x)
        a<-.C("temporal_non_lin_gauss_smooth",
              as.single(aperm(x,c(4,3,2,1))),
              as.integer(d),
              as.single(aperm(mask,c(3,2,1))),
              as.single(radius),
              as.single(sm),
              as.single(voxdim),
              res=single(length(x)),
              PACKAGE = "AnalyzeFMRI")
        
        a<-array(a$res,dim=d[4:1])
        a<-aperm(a,c(4,3,2,1))

        return(a)
    }
    
}    

## functions for fitting spatial and non-spatial mixture models to fMRI datasets
## the mixture model considered is a mixture of a standard normal distribution and two Gamma functions, denoted N2G
## x ~ p1 * N(0, 1) + p2 * Gamma(a, b) + (1 - p1 -p2) * -Gamma(c, d)
## par = c(a, b, c, d, p1, p2)


N2G <- function(data, par.start = c(4, 2, 4, 2, 0.9, 0.05)) {

    ## fits the N2G model to data using par.start as the starting point of the optimisation
    ## lims is the interval in which observations are assigned to the Normal component
    
    data <- data[data != 0]
    fit <- N2G.Fit(data, par.start, maxit = 500, method = "BFGS")
    lims <- N2G.Region(fit)
    
    return(list(par = fit, lims = lims))
}

N2G.Transform <- function(par) {
    
    ## transform parameters of N2G model so as to lie on the real line
    
    q <- par
    q[1:4] <- log(par[1:4])
    q[5] <- log(par[5] / (1 - par[5] - par[6]))
    q[6] <- log(par[6] / (1 - par[5] - par[6]))
    if(length(par) == 7) q[7] <- log(par[7])
    return(q)
}

N2G.Inverse <- function(par) {
    
    ## transform parameters back to their real domains
    
    q <- par
    q[1:4] <- exp(par[1:4])
    q[5] <- exp(par[5]) / (1 + exp(par[5]) + exp(par[6]))
    q[6] <- exp(par[6]) / (1 + exp(par[5]) + exp(par[6]))
    if(length(par) == 7) q[7] <- exp(par[7])
    
    return(q)
}

N2G.Density <- function (data, par) {
    
    ## density function for the N2G model
    
    pos <- data > 0
    neg <- data < 0
    d <- par[5] * dnorm(data)
    if(par[1] > 0) d[pos] <- d[pos] + par[6] * dgamma(data[pos], par[1], par[2])
    if(par[3] > 0) d[neg] <- d[neg] + (1 - par[5] - par[6]) * dgamma(-data[neg], 
        par[3], par[4])
    return(d)
    
}


N2G.Likelihood <- function(inv.par, data) {

    ## (Negative) Likelihood of the N2G model
    par <- N2G.Inverse(inv.par)
    lik <- -sum(log(N2G.Density(data, par)))

    return(lik)
}

N2G.Class.Probability <- function(data, par){

    ## Posterior Probability of data points being in each class 

    res <- matrix(0, length(data), 3)
    pos <- data > 0
    neg <- data < 0
    res[, 1] <- par[5] * dnorm(data)
    res[pos, 2] <- par[6] * dgamma(data[pos], par[1], par[2])
    res[neg, 3] <- (1 - par[5] - par[6]) * dgamma(-data[neg], par[3], par[4])

    rs <- rowSums(res)
    res <- res / rs
    
    return(res)
}

N2G.Fit <- function(data, par.start, maxit, method){

    ## main fitting function for N2G model
    
    inv.par <- N2G.Transform(par.start)

    ans <- optim(inv.par, fn = N2G.Likelihood, method = method, data = data, control = list(maxit = maxit, trace = 0))
    
    if(ans$convergence!=0) return("No convergence")
    res <- N2G.Inverse(ans$par)

    return(res)
}

N2G.Region <- function(par1) {

    ## calculates the interval within which observations are classified as belonging to the Normal component
    
    f1 <- function(x, p) (p[5] * dnorm(x) - p[6] * dgamma(x, p[1], p[2]))^2
    f2 <- function(x, p) (p[5] * dnorm(x) - (1 - p[5] - p[6]) * dgamma(-x, p[3], p[4]))^2

    ans1 <- optim(par = 1, fn = f1, method = "L-BFGS-B", p = par1, lower = 0, upper = 6, control = list(maxit = 500, trace = 0))$par
    ans2 <- optim(par = -1, fn = f2, method = "L-BFGS-B", p = par1, lower = -6, upper = 0, control = list(maxit = 500, trace = 0))$par
    
    return(c(ans1, ans2))
}



N2G.Likelihood.Ratio <- function(data, par) {

    ## calculates the ratio of the likelihood that data came from the positive Gamma distribution (activation) to the likelihood that data came from the other two distributions (Normal and negative Gamma)
    
    d1 <- array(0, dim = dim(data))
    d0 <- array(0, dim = dim(data))

    p.pos <- par[6]
    p.0   <- par[5]
    p.neg <- 1 - par[5] - par[6]
    
    pos <- data > 0
    neg <- data < 0
    
    d1[pos] <- dgamma(data[pos], par[1], par[2])

    d0 <- (p.0 / (1 - p.pos)) * dnorm(data)
    d0[neg] <- d0[neg] + (p.neg / (1 - p.pos)) * dgamma(-data[neg], par[3], par[4])

    ans <- d1 / d0
    
    return(ans)
}

model.2.cov.func <- function(g, par) {
    
    p.pos <- par[6]
    p0 <- par[5]
    p.neg <- 1 - p0 - p.pos
    
    mu.pos <- par[1] / par[2]
    mu.neg <- - par[3] / par[4]

    a1 <- (mu.neg * p.neg / (p0 + p.neg))^2 * (1 - p.pos * ((1 + g) - 1 / (1 + g)) / g)
    a2 <- 2 * (mu.neg * p.neg / (p0 + p.neg)) * mu.pos * p.pos / (1 + g)
    a3 <- mu.pos * mu.pos * p.pos * g / (1 + g)
    a4 <- p.neg * mu.neg + p.pos * mu.pos

    res <- a1 + a2 + a3 - a4^2
    
    return(res)
}

model.2.est.gamma <- function(cov, par) {
    
    tmp.func <- function(g, cov, par) (cov - model.2.cov.func(g, par))^2
    a <- optimize(f = tmp.func, interval = c(-100, 100), cov = cov, par = par)
    
    return(a)
}

N2G.Spatial.Mixture <- function(data, par.start = c(4, 2, 4, 2, 0.9, 0.05), ksize, ktype = c("2D", "3D"), mask = NULL) {

  ## non-spatial n2g fit ##
  if(is.null(mask)) mask <- data != 0
  data.m <- data[mask == 1]
  fit <- N2G(data.m, par.start)
  ##cat(fit$par)
  p <- fit$par[6]
  
  pos.act.map <- data > fit$lims[1]

  ## lr ##
  lik.ratio <- N2G.Likelihood.Ratio(data, fit$par)
  d <- dim(lik.ratio)
  
  ## neighbourhood (k = 8) ##
  nmat <- matrix(unlist(expand.grid(-1:1, -1:1, 0)), 9, 3)
  nmat <- nmat[-5, ]

  ## number of neighbours ##
  k <- dim(nmat)[1]
  ## possible values of s = number of active voxels in neighbourhood, including central voxel ##
  x <- 0:(k + 1)
  ## number of voxels with s active voxels in the neighbourhood ##
  y <- vector(len = k + 2)

  ## calculate y ##
  for(i in 2:(d[1] - 1))
  {
    for(j in 2:(d[2] - 1))
    {
      for(l in 1:d[3])
      {
        s <- sum(pos.act.map[t(c(i, j, l) + t(nmat))])
        s <- s + pos.act.map[i, j, l]
        y[s + 1] <- y[s + 1] + 1
      }
    }
  }

  ## turn y into a proportion ##
  y <- y / sum(y)

  model.2.func <- function(x, p, gamma, k) {
    ans <- (x == 0) * (1 - (p * ((1 + gamma)^(k + 1) - 1) / (gamma * (1 + gamma)^k)))
    ans <- ans + ((x > 0) * p * gamma^(x - 1) / (1 + gamma)^k)
    
    return(ans)
  }
  
  loglik.func <- function(gamma, x, y, p, k)  {
    -sum(y * log(model.2.func(x, p, gamma, k)))
  }
  
  o <- optim(par = c(gamma = 0.5), fn = loglik.func, method = "L-BFGS-B", lower = c(0), upper = c(Inf), control = list(maxit = 10000, trace = 0), x = x, y = y, k = k, p = p)

  ## parameter estimates ##
  gamma <- o$par[[1]]
  ##cat(gamma)
  
  kt <- switch(ktype[1], "2D" = 2, "3D" = 3)

  ## spatial mixture model ##
  a <- .C("spatial_mixture",
          as.double(aperm(lik.ratio, c(3, 2, 1))),
          as.integer(d),
          as.integer(ksize),
          as.integer(aperm(mask, c(3, 2, 1))),
          as.integer(kt),
          as.double(gamma),
          as.double(p),
          ans = double(prod(d)),
          PACKAGE = "AnalyzeFMRI")
    
  a1 <- array(a$ans, dim = d[3:1])
  a1 <- aperm(a1, c(3, 2, 1))
    
  return(list(p.map = a1, par = fit$par, lims = fit$lims, gamma = gamma, p = p))
}


cov.est <- function(mat, mask, nmat) {

    ## estimate covariance between neighbouring voxels
    
    a <- .C("covariance_est",
            as.double(aperm(mat, c(3, 2, 1))),
            as.integer(dim(mat)),
            as.integer(aperm(mask, c(3, 2, 1))),
            as.integer(t(nmat)),
            as.integer(dim(nmat)),
            ans = double(1),
            PACKAGE = "AnalyzeFMRI")

    return(a$ans)
}




cluster.threshold <- function(x, nmat = NULL, level.thr = 0.5, size.thr) {

  ## thresholds an array at level.thr
  ## calculates the number of contiguous clusters and their sizes
  ## answer is an array in which all voxels that are contained clusters of size greater
  ## than or equal to size.thr are 1, otherwise 0.
  ## nmat is a (Kx3) matrix specifying the neighbourhood system
  ## i.e if a row of nmat is (0, 1, -1) then x[10, 10, 10] and x[10, 11, 9] are neighbours

  if(is.null(nmat)) { ## default is 6 adjacent neighbours
    nmat <- expand.grid(-1:1, -1:1, -1:1)
    nmat <- nmat[c(5, 11, 13, 15, 17, 23), ]
  }   
  
  res <- .C("cluster_mass",
            mat = as.single(aperm(x, c(3, 2, 1))),
            as.integer(dim(x)),
            as.integer(t(nmat)),
            as.integer(dim(nmat)),
            as.integer(level.thr),
            num.c = integer(1),
            res.c = single(1000 * 6),
            PACKAGE = "AnalyzeFMRI")

  res.c <- matrix(res$res.c, 1000, 6, byrow = TRUE)[1:res$num.c, ]

  mat1 <- array(res$mat, dim = dim(x)[3:1])
  mat1  <-  aperm(mat1, c(3, 2, 1))

  m <- (res.c[, 5] < size.thr) * (1:res$num.c)
  m <- m[m != 0]

  for(i in 1:length(m)) 
    mat1[mat1 == m[i]] <- 0
  
  mat1 <- 1 * (mat1 > 0)
  
  return(mat1)
}
## Functions for simulating and thresholding Gaussian Random Fields 

Sim.3D.GRF <- function(d, voxdim, sigma, ksize, mask = NULL, type = c("field", "max")) {

    ## simulates a GRF with covariance matrix sigma of dimension d with voxel dimensions voxdim
    ## if type = "max" then just the maximum of the field is returned
    ## if type = "filed" then just the filed AND the maximum of the field are returned
    
    if(!is.null(mask) && sum(d[1:3] == dim(mask)) < 3)  return("mask is wrong size")
    if(length(d) != 3 && length(d) != 4) return("array x should be 3D or 4D")
    if(length(d) == 3) {
        d <- c(d, 1)
        tmp <- 1
    }
    
    if((2 * floor(ksize / 2)) == ksize) stop(paste("ksize must be odd"))
    
    if(is.null(mask)) mask <- array(1, dim = d[1:3])
    space <- 1 + (type == "field") * prod(d)
    
    filtermat <- GaussSmoothKernel(voxdim, ksize, sigma)
    
    a <- .C("sim_grf",
            as.integer(d),
            as.double(aperm(filtermat, c(3, 2, 1))),
            as.integer(ksize),
            as.integer(aperm(mask, c(3, 2, 1))),
            as.integer((type == "field")),
            mat = double(space),
            max = double(1),
            PACKAGE = "AnalyzeFMRI")
    
    if(type == "field") {
        mat <- array(a$mat, dim = d[4:1])
        mat <- aperm(mat, 4:1)
        if(tmp == 1) mat <- mat[, , , 1]
        return(list(mat = mat, max = a$max))
    }
    
    return(list(mat = NULL, max = a$max))
    
    
}


SmoothEst <- function (mat, mask, voxdim, method = "Forman")
{
  ## Estimate the variance-covariance matrix of a Gaussian random field

  ## set-up ## 
  x <- dim(mat)[1]
  y <- dim(mat)[2]
  z <- dim(mat)[3]
  b1 <- array(0, dim = dim(mat) + 2)
  b2 <- array(0, dim = dim(mat) + 2)
  m1 <- array(0, dim = dim(mat) + 2)
  m2 <- array(0, dim = dim(mat) + 2)
  
  ## X ##
  b1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mat
  b2[1:(x)    , 2:(y + 1), 2:(z + 1)] <- mat
  m1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mask
  m2[1:(x)    , 2:(y + 1), 2:(z + 1)] <- mask
  m3 <- (m1 + m2) == 2
  if(method == "Forman") {
      x.v0 <- mean((b1[m1 == 1])^2)
      x.v1 <- mean(((b1 - b2)[m3 == 1])^2)
      xx <- -(voxdim[1]^2) / (4 * log(1 - x.v1 / (2 * x.v0)))
  } 
  if(method == "Friston") 
      xx <- mean(((b1 - b2)[m3 == 1])^2) / (voxdim[1]^2)
  
  
  ## Y ##
  b1[, , ] <- 0
  b2[, , ] <- 0
  m1[, , ] <- 0
  m2[, , ] <- 0
  b1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mat
  b2[2:(x + 1), 1:(y)    , 2:(z + 1)] <- mat
  m1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mask
  m2[2:(x + 1), 1:(y)    , 2:(z + 1)] <- mask
  m3 <- (m1 + m2) == 2
  if(method == "Forman") {
    y.v0 <- mean((b1[m1 == 1])^2)
    y.v1 <- mean(((b1 - b2)[m3 == 1])^2)
    yy <- -(voxdim[2]^2) / (4 * log(1 - y.v1 / (2 * y.v0)))
  }
  if(method == "Friston")
      yy <- mean(((b1 - b2)[m3 == 1])^2) / (voxdim[2]^2)
  
  ## Z ##
  b1[, , ] <- 0
  b2[, , ] <- 0
  m1[, , ] <- 0
  m2[, , ] <- 0
  b1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mat
  b2[2:(x + 1), 2:(y + 1), 1:(z)    ] <- mat
  m1[2:(x + 1), 2:(y + 1), 2:(z + 1)] <- mask
  m2[2:(x + 1), 2:(y + 1), 1:(z)    ] <- mask
  m3 <- (m1 + m2) == 2
  if(method == "Forman") {
      z.v0 <- mean((b1[m1 == 1])^2)
      z.v1 <- mean(((b1 - b2)[m3 == 1])^2)
      zz <- -(voxdim[3]^2) / (4 * log(1 - z.v1 / (2 * z.v0)))
    }
  if(method == "Friston") 
      zz <- mean(((b1 - b2)[m3 == 1])^2) / (voxdim[3]^2)
  
  ## output ##
  if(method == "Forman") 
    sigma <- diag(c(xx, yy, zz))
  if(method == "Friston")
    sigma <- solve(diag(2 * c(xx, yy, zz)))
  
  return(sigma)
  
}


Threshold.Bonferroni <- function(p.val, n, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL) {

    ## calculate the Bonferroni threshold for n iid tests to give a p-value of p.val
    ## type specifies the univariate distribution of the test statistics under consideration
    
    if(type == "Normal") return(qnorm(1 - p.val / n))
    if(type == "t") return(qt(1 - p.val / n, df = df1))
    if(type == "F") return(qf(1 - p.val / n, df1 = df1, df2 = df2))

}

EC.3D <- function(u, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL) {

    ## The Expectation of the Euler Characteristic for a 3D Random Field above a threshold u
    ## type specifies the marginal distribution of the field
    
    V <- prod(voxdim) * num.vox
    S <- sqrt(det(solve(2 * sigma)))
    EC <- switch(type[1],
                 Normal = V * S * (u^2 - 1) * exp(-u^2 / 2) / (4 * pi^2),
                 t = V * S * (1 + (u^2 / df))^(-(df - 1) / 2) * (u^2 * (df - 1) / df - 1) / (4 * pi^2)
                 )
    return(EC)
}


Threshold.RF <- function(p.val, sigma, voxdim = c(1, 1, 1), num.vox, type = c("Normal", "t"), df = NULL) {

    ## calculates the Random Field theory threshold to give a p-value of p.val
    ## type specifies the marginal distribution of the field
    
    EC.func <- function(u , sigma, voxdim, num.vox, p.val) (EC.3D(u, sigma, voxdim, num.vox, type[1], df) - p.val)^2

    threshold <- optimize(f = EC.func, interval = c(2,10), sigma = sigma, voxdim = voxdim, num.vox = num.vox, p.val = p.val, tol = 1e-9)$minimum
    
    return(threshold)
}


Threshold.FDR <- function(x, q, cV.type = 2, type = c("Normal", "t", "F"), df1 = NULL, df2 = NULL) {

    ## calculates the FDR threshold for a vector of p-values in x
    ## q specifies the desired FDR
    ## cV specfies the type of FDR threshold used (See Genovese et al. (2002))
    
    if(type == "Normal") p <- sort(1 - pnorm(x))
    if(type == "t") p <- sort(1 - pt(x, df = df1))
    if(type == "F") p <- sort(1 - pf(x, df1 = df1, df2 = df2))

    V <- length(p)
    cV <- switch(cV.type, 1, log(V) + 0.5772)
    
    i <- 1
    while (p[i] <= (i * q) / (V * cV)) i <- i + 1
    i <- max(i - 1, 1)

    if(type == "Normal") thr <- qnorm(1 - p[i])
    if(type == "t") thr <- qt(1 - p[i], df = df1)
    if(type == "F") thr <- qf(1 - p[i], df1 = df1, df2 = df2)

    return(thr)

}

Sim.3D.GammaRF <- function(d, voxdim, sigma, ksize, mask, shape, rate) {

  ## simulates a smooth Gamma distributed random field by simulating a GRF and
  ## transforming each statistic value to be a Gammma
  
  field <- Sim.3D.GRF(d = d, voxdim = voxdim, sigma = sigma, ksize = 9,
                      mask = mask, type = "field")$mat

  gamma.field <- qgamma(pnorm(field), shape = shape, rate = rate)
  gamma.field <- gamma.field * mask

  return(gamma.field)

}
.First.lib <- function(lib, pkg){
    library.dynam("AnalyzeFMRI", pkg, lib)
}
