
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)
 print("",quote=F)
 print(paste("       File name:",file.img),quote=F)
 print(paste("  Data Dimension:",paste(hdr$dim[1],"-D",sep="")),quote=F)
 print(paste("     X dimension:",hdr$dim[2]),quote=F)
 print(paste("     Y dimension:",hdr$dim[3]),quote=F)
 print(paste("     Z dimension:",hdr$dim[4]),quote=F)
 print(paste("  Time dimension:",hdr$dim[5],"time points"),quote=F)
 print(paste("Voxel dimensions:",paste(hdr$pixdim[2],hdr$vox.units,"x",
                                   hdr$pixdim[3],hdr$vox.units,"x",
                                   hdr$pixdim[4],hdr$vox.units)),quote=F)
 print(paste("       Data type:",hdr$data.type,paste("(",hdr$bitpix," bits per voxel)",sep="")),quote=F)
}


f.analyzeIO.gui<-function(){

  #starts GUI that allows user to explore an fMRI dataset stored in an ANALYZEfile
  
  path<-.path.package(package="AnalyzeIO")
  path.gui<-paste(path,"/AnalyzeIOgui.R",sep="")
  source(path.gui)}


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.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="")


 
#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("read4byte1",
	ans = integer(1),
	file.hdr,as.integer(1),
	as.integer(1),
	as.integer(0))$ans==348) swap<-1 


# A C function is used to read in all the components of the .hdr file
a<-.C("read_analyze_header",
      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))

#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.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==4){
     vol<-.C("read2byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*2))
     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("read4byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*4))
     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("readfloat1",
             mat=single(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*4))
     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("readdouble1",
             mat=numeric(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*8))
     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==2||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=F)
       
   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==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("read2byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*2))
     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("read4byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*4))
     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("readfloat1",
             mat=single(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*4))
     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("readdouble1",
             mat=numeric(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(offset*8))
     vol<-array(vol$mat,dim=dim)
   }}
  
     if(hdr$datatype==0||hdr$datatype==1||hdr$datatype==2||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=F)
       
   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==4){
  for(i in 1:hdr$dim[5]){

     v<-.C("read2byte1",
             mat=integer(1),
             file.img,
             as.integer(hdr$swap),
             as.integer(1),
	     as.integer(offset.start*2+2*(i-1)*offset.add))
     vol[i]<-v$mat}
   }
  
  if(hdr$datatype==8){
  for(i in 1:hdr$dim[5]){

     v<-.C("read4byte1",
             mat=integer(1),
             file.img,
             as.integer(hdr$swap),
             as.integer(1),
	     as.integer(offset.start*4+4*(i-1)*offset.add))
     vol[i]<-v$mat}
   }
  
  if(hdr$datatype==16){
  for(i in 1:hdr$dim[5]){

     v<-.C("readfloat1",
             mat=integer(1),
             file.img,
             as.integer(hdr$swap),
             as.integer(1),
	     as.integer(offset.start*4+4*(i-1)*offset.add))
     vol[i]<-v$mat}
   }
  
  if(hdr$datatype==64){
  for(i in 1:hdr$dim[5]){

     v<-.C("readdouble1",
             mat=integer(1),
             file.img,
             as.integer(hdr$swap),
             as.integer(1),
	     as.integer(offset.start*8+8*(i-1)*offset.add))
     vol[i]<-v$mat}
   }
  
     if(hdr$datatype==0||hdr$datatype==1||hdr$datatype==2||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=F)
       
   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==4){
     vol<-.C("read2byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(0))
     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("read4byte1",
             mat=integer(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(0))
     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("readfloat1",
             mat=single(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(0))
     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("readdouble1",
             mat=numeric(num.data.pts),
             file.img,
             as.integer(hdr$swap),
             as.integer(num.data.pts),
	     as.integer(0))
     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==2||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=F)
       
   return(vol)}


f.spectral.summary<-function(file,mask.file,ret.flag=F)
{
  #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(length(slices),xc,yc))

  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[k,b,a]<-mean(slice[b,a,])}

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

  return(mask)

}

  if(mask.file!=F){mask<-f.read.analyze.volume(mask.file)}
  else{
dat<-list(file=file,mask.file=mask.file)
    mask<-f.mask.create(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
#####################
  
  for(l in 1:nsl){
    print(c("Processing slice",l))
 
    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
      }
    }
  }
    
b<-apply(res,4,FUN=quantile,probs=seq(.5,1,.05),na.rm=T)
  plot(c(0,n-1),c(0,30),type="n",xlab="",ylab="",axes=F)
  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==T)  return(b)
}





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)
  }
  
  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",
     as.integer(mat),
     file,
     as.integer(num.data.pts))

  }



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",
     as.single(mat),
     file,
     as.integer(num.data.pts))
}

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

# Writes a list to a .hdr file
a<-.C("write_analyze_header",
      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))
}

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