library(ANTsR)
 tasks=c("MR290","MR")
 study="/data/asimmons/data/CDA/ANTs/"
 tem<-antsImageRead('/data/asimmons/data/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c.nii.gz')
 temmask<-antsImageRead('/data/asimmons/data/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c_mask.nii.gz')
 setwd(study)
 subs<-Sys.glob("TC_150_2")
for (sub in subs){
  print(paste("Checking and processing ",sub))
  if(!(file.exists(paste(study,"/",sub,"/",sub,"_t1.nii.gz",sep="")))){
	setwd(paste(study,sub,sep="/"))
	t1<-antsImageRead(paste(sub,".nii.gz",sep=""),3)
	t1c<-abpN4(t1)
	antsImageWrite(t1c,paste(sub,"_c.nii.gz",sep=""))
	print("Skull Strip")
	print(system.time(bm<-abpBrainExtraction(t1c,tem,temmask)))
	bmbrain <- antsApplyTransforms(fixed=tem,moving=bm$brain,transformlist=bm$fwdtransforms )
	antsImageWrite(bm$brain,paste(sub,"_t1orig.nii.gz",sep=""))
	antsImageWrite(bmbrain,paste(sub,"_t1.nii.gz",sep=""))
	file.copy(bm$fwdtransforms[1],paste(sub,"_ft.nii.gz",sep=""))
	file.copy(bm$fwdtransforms[2],paste(sub,"_ft.mat",sep=""))
	file.copy(bm$invtransforms[1],paste(sub,"_invt.nii.gz",sep=""))
	file.copy(bm$invtransforms[2],paste(sub,"_invt.mat",sep=""))
	bmseg <- antsApplyTransforms(fixed=tem,moving=bm$kmeansseg,transformlist=bm$fwdtransforms )
	antsImageWrite(bmseg,paste(sub,"_seg.nii.gz",sep=""))}
 for (task in tasks) {
    taskout<-paste(task,"Prep",sep="")
    if(!(file.exists(paste(study,"/",sub,"/",sub,taskout,".nii.gz",sep="")))){
      if(file.exists(paste(study,"/",sub,"/",sub,task,".nii.gz",sep=""))){
	setwd(paste(study,sub,sep="/"))
	print(paste("running data for ",task," in ",sub,sep=""))
	if(!(file.exists(paste(sub,task,"_ds.nii.gz",sep=""))))
	 {system(paste("3dDespike -ignore 4 -prefix ",sub,task,"_ds.nii.gz ",sub,task,".nii.gz",sep=""))}
	if(!(file.exists(paste(sub,task,"_dst.nii.gz",sep=""))))
	 {system(paste("3dTshift -ignore 2 -tzero 0 -tpattern seqplus -prefix ",sub,task,"_dst.nii.gz ",sub,task,"_ds.nii.gz",sep=""))}
	bold=antsImageRead(paste(sub,task,"_dst.nii.gz",sep=""),4)
	print("n3 of BOLD")
	avg <- getAverageOfTimeSeries(bold)
	print(system.time(avgn4 <- abpN4(avg)))
	mask<-getMask(avgn4,mean(avgn4),Inf,2)
	bold_n3<-timeseriesN3(bold,mask)
	avgm<-as.vector(avgn4[mask==1])	
	boldstm<-timeseries2matrix(bold_n3,mask)
	boldstmp<-100*t(t(boldstm)/avgm)
	boldstmw<-temporalwhiten(boldstmp, myord = 2)
	boldstw<-matrix2timeseries(bold,mask,boldstmw)
	print("preprocess fMRI")
	bmbrainorig=antsImageRead(paste(sub,"_t1orig.nii.gz",sep=""))
	bmfwdtransforms1=paste(sub,"_ft.nii.gz",sep="")
	bmfwdtransforms2=paste(sub,"_ft.mat",sep="")
	print(system.time(boldpre<-preprocessfMRI( boldstw, mask, useMotionCorrectedImage = TRUE)))
	boldtx<-antsRegistration(fixed=avg, moving=bmbrainorig, typeofTransform="SyNBold",
		spatialSmoothingType="perona-malik",spatialSmoothingParameters=c(0.25, 5))
	fulltx<-c(bmfwdtransforms1,bmfwdtransforms2,boldtx$invtransforms[1],boldtx$invtransforms[2])
	temrs<-resampleImage(tem,c(4,4,4),interpType=4)
	antsImageWrite(boldpre$cleanBoldImage,paste(sub,taskout,".nii.gz",sep=""))
	boldpreM<-antsApplyTransforms(fixed=temrs,moving=boldpre$cleanBoldImage,transformlist=fulltx,whichtoinvert=c(F,F,T,F),imagetype=3)
	antsImageWrite(boldpreM,paste(sub,taskout,"MNI.nii.gz",sep=""))
	maskM<-antsApplyTransforms(fixed=temrs, moving=mask, transformlist=fulltx, whichtoinvert=c(F,F,T,F))
	antsImageWrite(maskM,paste(sub,taskout,"mask.nii.gz",sep=""))
	boldpreMa=getAverageOfTimeSeries(boldpreM)
	boldpreMa=boldpreMa*boldpreMa
	write.csv(boldpre$nuisanceVariables,paste("NV",sub,taskout,"1D",sep="."))
	pdf(paste("Prep",sub,taskout,"pdf",sep="."))
	plot(t1,avgn4,window.overlay=c(mean(avgn4),Inf))
	plot(boldpre$DVARS,type='l',ylab='DVARS',xlab='TR')
	plot(bmbrain,boldpreMa,axis=3)
	dev.off()
	}}
  try(rm(list=c("bold","boldpre","boldstm","boldstmw","boldpreM")))
  }}