library(ANTsR)
library(robustbase)
task<-("MR")
taskout<-("MRttl")
study<-"/data/asimmons/data/CDA/ANTs"
setwd(study)
subs<-Sys.glob("TC_150_2")
for (sub in subs){
	print(paste("working on ",sub))
      if(!(file.exists(paste(study,"/",sub,"/",sub,taskout,"_rANTS.nii.gz",sep="")))){
      if(file.exists(paste(study,"/",sub,"/",sub,task,"PrepMNI.nii.gz",sep=""))){
	setwd(paste(study,sub,sep="/"))
	if(file.exists(paste(task,"ttlRwav",sep=""))){reg<-read.delim(paste(task,"ttlRwav",sep=""),header=TRUE)}
	if(file.exists(paste(sub,task,"ttlRwav",sep=""))){reg<-read.delim(paste(sub,task,"ttlRwav",sep=""),header=TRUE)}
	system(paste("3dToutcount -automask ",sub,task,".nii.gz > ",sub,task,"_OL.txt",sep=""))
	system(paste("PkCensorSD ",sub,task,"_OL.txt 1.5 ",sub,task,"Censor",sep=""))
	Censor<-read.delim(paste(sub,task,"Censor",sep=""),header=FALSE)
	boldpre<-antsImageRead(paste(sub,task,"PrepMNI.nii.gz",sep=""),4)
	  mask<-antsImageRead(paste(sub,task,"Prepmask.nii.gz",sep=""),pixeltype="unsigned int")
	  avg<-getAverageOfTimeSeries(boldpre)
  	  amask<-avg*avg
	  amask[amask!=0]<-1
  	  tmask=amask*mask
	boldmat<-timeseries2matrix(boldpre,tmask)
	regclip<-reg[1:dim(boldpre)[4],]
	regc<-regclip[Censor==1,]
#	myrlm <- lapply(1:ncol(boldmat), function(x) (summary(lmrob(boldmat[Censor==1,x]~regc$EasyC+regc$MedC+regc$HardC+regc$EasyE+regc$MedE+regc$HardE+regc$Easy+regc$Med+regc$Hard, singular.ok=T))$coefficients))
	myrlm <- lapply(1:ncol(boldmat), function(x) (summary(lmrob(boldmat[Censor==1,x]~regc$EasyMed+regc$Hard, singular.ok=T))$coefficients))
	for (i in 2:(nrow(myrlm[[1]]))) {
	 myrlm_b<-sapply (1:length(myrlm), function(x) myrlm[[x]][i,1])
	 beta<-makeImage(tmask,myrlm_b)
	 antsImageWrite(beta,paste(sub,taskout,"_b",i,".nii.gz",sep=""))
	 myrlm_b<-sapply (1:length(myrlm), function(x) myrlm[[x]][i,3])
	 beta<-makeImage(tmask,myrlm_b)
	 antsImageWrite(beta,paste(sub,taskout,"_b",i,"t.nii.gz",sep=""))}
	system(paste("3dbucket -prefix ",sub,taskout,"_rANTS.nii.gz ",sub,taskout,"_b*.nii.gz",sep=""))
	system(paste("rm ",sub,taskout,"_b*.nii.gz",sep=""))
	refit<-""
	for (i in 2:(nrow(myrlm[[1]]))) {
		refit_b <-paste(" -sublabel",2*(i-1)-2,rownames(data.frame(myrlm[[1]]))[i])
		refit_t <-paste(" -sublabel ",2*(i-1)-1," t_",rownames(data.frame(myrlm[[1]]))[i],sep="")
		refit<-paste(refit,refit_b,refit_t)}
	refit<-gsub("sub","-sub",gsub("regc","",gsub("[[:punct:]]","",refit)))
	system(paste("3drefit -space MNI -view tlrc ",refit," ",sub,taskout,"_rANTS.nii.gz",sep=""))
    }
  }
}