diff --git a/anatomicuts/anatomiCutsUtils b/anatomicuts/anatomiCutsUtils
index ba2f84ddd35..8319978865a 100644
--- a/anatomicuts/anatomiCutsUtils
+++ b/anatomicuts/anatomiCutsUtils
@@ -14,6 +14,7 @@ import pandas as pd
 import matplotlib as mpl
 mpl.use('Agg')
 import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
 import freesurfer as fs
 import argparse
 import graph_tools as gt
@@ -24,6 +25,7 @@ from holoviews import opts, dim
 import plotly.express as px#
 import statsmodels.api as sm
 
+colors=['y','b','r','g','coral','purple']
 def getCorrespondingClusters(correspondance,order=True):
 		corr=dict()
 		distances=dict()
@@ -46,8 +48,8 @@ def getCorrespondingClusters(correspondance,order=True):
 					else:
 						header=False
 		return corr, distances, indeces
-def thicknessPerStructure(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, lenght=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
-	ac_folder=f"dmri.ac/{lenght}/{std}"
+def thicknessPerStructure(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, length=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
+	ac_folder=f"dmri.ac/{length}/{std}"
 	clusterNum=200
 	lut = fs.LookupTable.read("/usr/local/freesurfer/dev/FreeSurferColorLUT.txt")
 	print(lut[28].name)
@@ -56,6 +58,7 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 	print(classification_cols, classificationFile, delimiter)
 	groups_cat = {row[classification_cols[0]] :row[classification_cols[1]] for _, row in pd.read_csv(classificationFile).iterrows()}
 	print(groups_cat)
+	motion = {row[0] :row[1] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/outliers_motion.txt", delimiter=' ').iterrows()}
 	significant_childs=set()
 	thicknesses=[[],[],[],[]]
 
@@ -65,15 +68,15 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 
 	ww=0
 	for s , cat2 in groups_cat.items():
-	
-		if cat2 > 10:
+		print(cat2)	
+		if int(cat2) > 10:
 			cat = 0
 		else:
-			cat =  cat2
+			cat =  int(cat2)
 		cat=cat2
 		print(s,cat)
 		gres = glob.glob(f'{subjects_dir}/{s}*')
-		if len(gres)>0 and cat< len(thicknesses):
+		if len(gres)>0 and cat< len(thicknesses) and s in motion and int(motion[s])<350:
 			subject=gres[0].split("/")[-1]
 			try:
 				correspondences= f"{subjects_dir}/{subject}/{ac_folder}/match/{target_subject}_{subject}_c{clusterNum}_hungarian.csv"
@@ -82,9 +85,10 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 					values=pd.read_csv(f"{subjects_dir}/{subject}/{ac_folder}/measures/surf/"+corr[indeces[i]]+".csv")
 					#print(values.shape[0])
 					for j in range(values.shape[0]):
-						val = 1+  abs(float(values['curv.start'][j])) 
+						val = 1 #+  abs(float(values['curv.start'][j])) 
+						val2 =   abs(float(values['thickness.start'][j])) 
 						#val=1
-						if val >0.01:
+						if val2 >0.01:
 							sign= '+'
 							#if abs(float(values['curv.start'][j] )) >0.1:
 							#	sign = '+'
@@ -96,9 +100,10 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 							else:
 								thicknesses[cat][i][sign+str(values['label.start'][j])]= [float(values['thickness.start'][j])/val]
 
-						val =  1+ abs(float(values['curv.end'][j])) 
+						val =  1 #+ abs(float(values['curv.end'][j])) 
+						val2 =  1+ abs(float(values['thickness.end'][j])) 
 						#val=1
-						if val >0.01:
+						if val2 >0.01:
 							sign= '+'
 							#if abs(float(values['curv.end'][j]))  >0.1:
 							#	sign = '+'
@@ -131,30 +136,31 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 		thicknessPolar=[]
 		catsPolar=[]
 		for ks,v in thicknesses[0][i].items():
-			ku = ks[1:]
-			if   ks in thicknesses[1][i]  and  ks in thicknesses[2][i]  and str(ku).isnumeric() :
-				print(	 len(v)*100 /totals[0],  len(thicknesses[1][i][ks])*100/totals[1] , len(thicknesses[2][i][ks])*100/totals[2] ) 
-			if int(ku) > 0 and  len(v)*100 /totals[0] > 2 and ks in thicknesses[1][i]  and len(thicknesses[1][i][ks])*100/totals[1] > 2 and ks in thicknesses[2][i]  and len(thicknesses[2][i][ks])*100/totals[2] > 2  and str(ku).isnumeric() :
-				toplot.append(v)				
-				toplot.append(thicknesses[1][i][ks])				
-				toplot.append(thicknesses[2][i][ks])				
-				labels.append(lut[int(ku)].name+str(ks[0]))
-				labels.append("")
-				labels.append("")
-				print(lut[int(ku)].name , int(ku), ks)
-				colors.append('r')
-				colors.append('b')
-				colors.append('g')
-				labelsPolar.append(lut[int(ku)].name+str(ks[0]))
-				labelsPolar.append(lut[int(ku)].name+str(ks[0]))
-				labelsPolar.append(lut[int(ku)].name+str(ks[0]))
-				thicknessPolar.append(np.mean(v))				
-				thicknessPolar.append(np.mean(thicknesses[1][i][ks]))				
-				thicknessPolar.append(np.mean(thicknesses[2][i][ks]))			
-				catsPolar.append("H")
-				catsPolar.append("MD")
-				catsPolar.append("D")
-		
+			if str(ks[1:]).isnumeric():
+				ku = int(float(ks[1:]))
+			
+				if int(ku) > 0 and  len(v)*100 /totals[0] > 5 and ks in thicknesses[1][i]  and len(thicknesses[1][i][ks])*100/totals[1] > 5 and ks in thicknesses[2][i]  and len(thicknesses[2][i][ks])*100/totals[2] > 5  :
+					print(i,totals,	 len(v)*100 /totals[0],  len(thicknesses[1][i][ks])*100/totals[1] , len(thicknesses[2][i][ks])*100/totals[2] ) 
+					toplot.append(v)				
+					toplot.append(thicknesses[1][i][ks])				
+					toplot.append(thicknesses[2][i][ks])				
+					labels.append(lut[int(ku)].name+str(ks[0]))
+					labels.append("")
+					labels.append("")
+					print(lut[int(ku)].name , int(ku), ks)
+					colors.append('r')
+					colors.append('b')
+					colors.append('g')
+					labelsPolar.append(lut[int(ku)].name+str(ks[0]))
+					labelsPolar.append(lut[int(ku)].name+str(ks[0]))
+					labelsPolar.append(lut[int(ku)].name+str(ks[0]))
+					thicknessPolar.append(np.mean(v))				
+					thicknessPolar.append(np.mean(thicknesses[1][i][ks]))				
+					thicknessPolar.append(np.mean(thicknesses[2][i][ks]))			
+					catsPolar.append("H")
+					catsPolar.append("MD")
+					catsPolar.append("D")
+			
 				
 		if len(toplot) >0:
 			parts =	plt.violinplot(toplot, showextrema=False)
@@ -164,29 +170,109 @@ def thicknessPerStructure(classificationFile, classification_cols, target_subjec
 				pc.set_alpha(1)
 			plt.gcf().set_size_inches(len(labels),5)	
 			plt.xticks(range(1,len(labels)+1), labels)
-			plt.savefig(f"{subjects_dir}/average/{i}.png")
+			plt.savefig(f"{subjects_dir}/average/dmri.ac/{length}/{std}/{i}.png")
 		
 			thick = pd.DataFrame({'Structure':labelsPolar, 'Thickness':thicknessPolar, 'Cat':catsPolar})
 			fig = px.line_polar(thick, r="Thickness", theta="Structure", color="Cat", line_close=True)
-			fig.write_html(f"{subjects_dir}/average/{i}.html")
+			fig.write_html(f"{subjects_dir}/average/dmri.ac/{length}/{std}/{i}.html")
+
+
+def diffAlongTheLine(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, length=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
+	ac_folder=f"dmri.ac/{length}/{std}"
+	clusterNum=200
+	lut = fs.LookupTable.read("/usr/local/freesurfer/dev/FreeSurferColorLUT.txt")
+	print(lut[28].name)
+
+	indeces=[]
+	print(classification_cols, classificationFile, delimiter)
+	groups_cat = {row[classification_cols[0]] :row[classification_cols[1]] for _, row in pd.read_csv(classificationFile).iterrows()}
+	print(groups_cat)
+	motion = {row[0] :row[2] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/motion_hcp.txt", delimiter=',').iterrows()}
+	significant_childs=set()
+	alongTheLine=[[],[],[],[]]
+	diagnosis=["HC","HD","MCID","DD"]
+	measures=["FA","MD","AD","RD"]
+	for j in range(len(alongTheLine)):
+		for k in range(len( measures)):
+			alongTheLine[j].append([])
+			for i in range(200):
+				alongTheLine[j][k].append([])
 
-	#except Exception as e :
-	#	print(e)
-	#bash /space/erebus/2/users/vsiless/code/freesurfer/anatomicuts/dmri_ac.sh GA BAKP64e_nrecon-trio3_vb17 40 5 /space/snoke/1/public/vivros/hd_tracula/labels_hd.csv  0:2 999999999 2 [0,1,2] True)
-def connectivityGraph(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, lenght=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
+	ww=0
+	for s , cat2 in groups_cat.items():
+		if int(cat2) == 7:
+			cat = 1000
+		elif int(cat2)==0:
+			cat = 0
+		elif int(cat2)==1:
+			cat = 1
+		elif int(cat2)==2 or int(cat2)==3:
+			cat =2
+		else:
+			cat = 1000
+		gres = glob.glob(f'{subjects_dir}/{s}*')
+		if len(gres)>0 and cat< len(alongTheLine)  and (s in motion and float(motion[s])<1.5):
+			subject=gres[0].split("/")[-1]
+			try:
+				correspondences= f"{subjects_dir}/{subject}/{ac_folder}/match/{target_subject}_{subject}_c{clusterNum}_hungarian.csv"
+				corr, distances, indeces =  getCorrespondingClusters(correspondences, False)
+				print(s,cat, cat2)
+				for i in range(200):
+					values=pd.read_csv(f"{subjects_dir}/{subject}/{ac_folder}/measures/surf/"+corr[indeces[i]]+".csv")
+					
+					if distances[indeces[i]] < 1e-5:
+						for mi, m  in enumerate(measures):
+							allVals=[]
+							for j in range(values.shape[0]):
+								vals=[]
+								for p in range(10):
+									val =   float(values[str(p)+"_"+m][j]) 
+									vals.append(val)
+								allVals.append(vals)
+							allVals= np.array(allVals).mean(axis=0)
+							alongTheLine[cat][mi][i].append(allVals)
+			except Exception as e:
+				print("ERROR",e)
+		ww+=1
+		#if ww >8:
+		#	break
+
+	for i in range(200):
+		total =0
 
+		fg , ax = plt.subplots(1,4, figsize=(20,5))
+		for mi, m in enumerate(measures):
+			for j in range(len(alongTheLine)):		
+				print(m, j, len(alongTheLine[j][mi][i]))
+				lines= np.array(alongTheLine[j][mi][i])
+				#for line in alongTheLine[j][mi][i]:		
+				#	plt.plot(line, color=colors[j])
+					
+				if len(lines)>0:
+					meanvals= lines.mean(axis=0)
+					#stdvals= lines.std(axis=0)
+					stdvals= stats.sem(lines,axis=0)
+					print(meanvals, stdvals)
+					ax[mi].plot(meanvals, color=colors[j], label=diagnosis[j])
+					ax[mi].fill_between(range(10),meanvals-stdvals, meanvals+stdvals, facecolor=colors[j], alpha=.5)
+		plt.legend()
+		plt.savefig(f"{subjects_dir}/average/dmri.ac/{length}/{std}/{i}_DTI.png")
+		
 
+def connectivityGraph(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, length=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
 
-	ac_folder=f"dmri.ac/{lenght}/{std}"
+	ac_folder=f"dmri.ac/{length}/{std}"
 	clusterNum=200
 	lut = fs.LookupTable.read("/usr/local/freesurfer/dev/FreeSurferColorLUT.txt")
 	print(lut[28].name)
 
 	indeces=[]
 	groups_cat = {row[classification_cols[0]] :row[classification_cols[1]] for _, row in pd.read_csv(classificationFile, delimiter=delimiter).iterrows()}
+	motion = {row[0] :row[1] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/outliers_motion.txt", delimiter=' ').iterrows()}
+
 	print(groups_cat)
 	significant_childs=set()
-	thicknesses=[[],[],[],[]]
+	thicknesses=[[],[],[],[], [],[]]
 
 	for j in range(len(thicknesses)):
 		thicknesses[j]= dict()
@@ -194,14 +280,11 @@ def connectivityGraph(classificationFile, classification_cols, target_subject, s
 	ww=0
 	for s , cat2 in groups_cat.items():
 	
-		if cat2 > 10:
-			cat =3 
-		else:
-			cat =  cat2
+		cat =  int(cat2)
 		#cat=cat2
 		print(s,cat)
 		gres = glob.glob(f'{subjects_dir}/{s}*')
-		if len(gres)>0 and cat< len(thicknesses):
+		if len(gres)>0 and cat< len(thicknesses) and s in motion and int(motion[s])<350:
 			subject=gres[0].split("/")[-1]
 			try:
 				correspondences= f"{subjects_dir}/{subject}/{ac_folder}/match/{target_subject}_{subject}_c{clusterNum}_hungarian.csv"
@@ -210,25 +293,28 @@ def connectivityGraph(classificationFile, classification_cols, target_subject, s
 					values=pd.read_csv(f"{subjects_dir}/{subject}/{ac_folder}/measures/surf/"+corr[indeces[i]]+".csv")
 					#print(values.shape[0])
 					for j in range(values.shape[0]):
-
-	
-						start = values['label.start'][j]
-						end = values['label.end'][j]
-						if start>0 and end>0:
-							ts = values['thickness.start'][j]
-							te = values['thickness.end'][j]
-							if int(start) < int(end):
-								end = start
-								start= values['label.end'][j]
-
-							if not start in thicknesses[cat]:
-								thicknesses[cat][start] = dict()
-							if not end in thicknesses[cat][start]:
-								thicknesses[cat][start][end]=[]
-							if len(thicknesses[cat][start][end])==0:
-								thicknesses[cat][start][end]=[]
-							thicknesses[cat][start][end].append((ts+te)/2)
+						try:		
+							start = int(values['label.start'][j])
+							end = int(values['label.end'][j])
+							if start>0 and end>0:
+								ts = float(values['thickness.start'][j])
+								te = float(values['thickness.end'][j])
+								if ts>0.1 and te>0.1:
+									if int(start) < int(end):
+										end = start
+										start= values['label.end'][j]
+
+									if not start in thicknesses[cat]:
+										thicknesses[cat][start] = dict()
+									if not end in thicknesses[cat][start]:
+										thicknesses[cat][start][end]=[]
+									if len(thicknesses[cat][start][end])==0:
+										thicknesses[cat][start][end]=[]
+									thicknesses[cat][start][end].append((ts+te)/2)
+						except Exception as e:
+							print("ERROR",e)
 			except Exception as e:
+				print(correspondences, values.shape[0], values.shape[1])
 				print("ERROR",e)
 		ww+=1
 	for cat in range(len(thicknesses)):
@@ -240,7 +326,7 @@ def connectivityGraph(classificationFile, classification_cols, target_subject, s
 		nodesId = []
 		nodesNames = []
 		nodesSet =set()
-		if len(thicknesses[cat]) >5:
+		if len(thicknesses[cat]) >1:
 
 			for start, items in thicknesses[cat].items():	
 				for end, values in items.items():
@@ -263,40 +349,47 @@ def connectivityGraph(classificationFile, classification_cols, target_subject, s
 			nodes = hv.Dataset(thenodes,'NodeID','NodeName')
 			chord = hv.Chord((links,nodes),['SourceID','DestinID'], ['values'])
 			chord.opts(opts.Chord(cmap='fire', edge_line_width=dim('values'), edge_color=dim('values'),height=200, labels='NodeName', width=200))
-			hv.save(chord,f"{subjects_dir}/average/cat_{cat}.html")
+			hv.save(chord,f"{subjects_dir}/average/dmri.ac/{length}/{std}/cat_{cat}.html")
 		
 
-def averageCorrespondingClusters(correspondences, imagesFolder, outputFolder,  clusterIndeces):
-		averages=dict()
-		norm=dict()
-
-		for s_i, correspondance in enumerate(correspondences):
-				try:
-					for clusterIndex in clusterIndeces:
-						corr, distances, indeces =  getCorrespondingClusters(correspondance, True)
-						#print(correspondance)
-						image=imagesFolder[s_i]+""+indeces[clusterIndex]+".nii.gz"
-						im =nib.load(image)
-						b= im.get_data()
-						b = b/b.max()
-						b = np.ceil(b) 
-						if clusterIndex in averages:
-							b += averages[clusterIndex].get_data() 
-							averages[clusterIndex]=nib.Nifti1Image(b, averages[clusterIndex].get_affine())
-							norm[clusterIndex]+=1
-						else:
-							averages[clusterIndex] = nib.Nifti1Image(b, im.get_affine())    
-							norm[clusterIndex]=1
-				except Exception as e:
-					print(str(e))
-
-		for clusterIndex in clusterIndeces:
-				directory=outputFolder
-				if not os.path.exists(directory):
-					os.makedirs(directory)
-				data=averages[clusterIndex].get_data()/norm[clusterIndex]
-				nib.Nifti1Image(data, averages[clusterIndex].get_affine()).to_filename(directory+"/"+str(clusterIndex)+'.nii.gz')
-				print ("saving",directory+"/"+str(clusterIndex)+'.nii.gz')
+def averageCorrespondingClusters(correspondences, imagesFolder, outputFolder,  clusterIndeces, subjectsNames , classificationFile, classification_cols,group, delimiter):
+		groups_cat = {row[classification_cols[0]] :row[classification_cols[1]] for _, row in pd.read_csv(classificationFile, delimiter=delimiter).iterrows()}
+
+		print(groups_cat)
+
+		for g in group:
+			averages=dict()
+			norm=dict()
+			print("group", g)	
+			for s_i, correspondance in enumerate(correspondences):
+				if subjectsNames[s_i] in groups_cat and groups_cat[subjectsNames[s_i]] in g:
+					try:
+						for clusterIndex in clusterIndeces:
+							corr, distances, indeces =  getCorrespondingClusters(correspondance, True)
+							#print(correspondance)
+							image=imagesFolder[s_i]+""+indeces[clusterIndex]+".nii.gz"
+							im =nib.load(image)
+							b= im.get_data()
+							b = b/b.max()
+							b = np.ceil(b) 
+							if clusterIndex in averages:
+								b += averages[clusterIndex].get_data() 
+								averages[clusterIndex]=nib.Nifti1Image(b, averages[clusterIndex].get_affine())
+								norm[clusterIndex]+=1
+							else:
+								averages[clusterIndex] = nib.Nifti1Image(b, im.get_affine())    
+								norm[clusterIndex]=1
+					except Exception as e:
+						print(str(e))
+
+			for clusterIndex in clusterIndeces:
+					directory=outputFolder
+					if not os.path.exists(directory):
+						os.makedirs(directory)
+					if clusterIndex in averages:
+						data=averages[clusterIndex].get_data()/norm[clusterIndex]
+						nib.Nifti1Image(data, averages[clusterIndex].get_affine()).to_filename(directory+"/"+str(clusterIndex)+'_'+str(g)+'.nii.gz')
+						print ("saving",directory+"/"+str(clusterIndex)+"_"+str(g)+'.nii.gz')
 
 def readTree(numNodes, histogramFile,header=True):
     almostFoundAllClusters=False
@@ -338,14 +431,15 @@ def readTree(numNodes, histogramFile,header=True):
             
     return nodes_childs, whos_dad
 
-def groupAnalysis( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, lenght, std, delimiter=",", groupA=[0], groupB=[1], thickness=False):
-	ac_folder=f"dmri.ac/{lenght}/{std}"
+def groupAnalysis( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, length, std, delimiter=",", groupA=[0], groupB=[1], thickness=False):
+	ac_folder=f"dmri.ac/{length}/{std}"
 	
 
 	indeces=[]
 	groups_cat = {row[classification_columns[0]] :row[classification_columns[1]] for _, row in pd.read_csv(groups_classification, delimiter=delimiter).iterrows()}
 	print(groups_cat)
 
+	motion = {row[0] :row[1] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/outliers_motion.txt", delimiter=' ').iterrows()}
 	significant_childs=set()
 	for c_i, clusterNum in enumerate(clustersToAnalyze):
 		
@@ -363,30 +457,33 @@ def groupAnalysis( headers, cols , groups_classification, classification_columns
 		X=[]
 		
 		for s in groups_cat.keys():
-			#print(f'{subjects_dir}/{s}*')
-			gres = glob.glob(f'{subjects_dir}/{s}*')
-			if len(gres)>0:
-				subject=gres[0].split("/")[-1]
-				if thickness:
-					measures=f"{subjects_dir}/{subject}/{ac_folder}/measures/surfaceMeasures.csv"
-				else:
-					measures=f"{subjects_dir}/{subject}/{ac_folder}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
-				if os.path.exists(measures):
-					data = pd.read_csv(measures,delimiter=",", header=0, names=headers, usecols=cols)
-					#print(measures, headers, cols)
-					if len(data[headers[0]])>= clusterNum:
-						#print(int(groups_cat[s]	), groupA, groupB)
-						if int(groups_cat[s]) in groupA:
-							X.append([1, 0])
-			
-						elif int(groups_cat[s]) in groupB:
-							X.append([0, 1])
+			if s in motion and int(motion[s]) <200:
+	
+				gres = glob.glob(f'{subjects_dir}/{s}*')
+				similarities = {row[0] :row[2] for _, row in pd.read_csv(f"{subjects_dir}/{subject}/{ac_folder}/match/{target_subject}_{subject}_c{clusterNum}.csv", delimiter=delimiter).iterrows()}
+
+				if len(gres)>0:
+					subject=gres[0].split("/")[-1]
+					if thickness:
+						measures=f"{subjects_dir}/{subject}/{ac_folder}/measures/surfaceMeasures.csv"
+					else:
+						measures=f"{subjects_dir}/{subject}/{ac_folder}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
+					if os.path.exists(measures):
+						data = pd.read_csv(measures,delimiter=",", header=0, names=headers, usecols=cols)
+						#print(measures, headers, cols)
+						if len(data[headers[0]])>= clusterNum:
+							#print(int(groups_cat[s]	), groupA, groupB)
+							if int(groups_cat[s]) in groupA:
+								X.append([1, 0])
 				
-						if int(groups_cat[s]) in np.concatenate( (groupA,groupB), axis=None):
-							#print("hola")
-							for i,h in enumerate(headers):
-								for j in range(clusterNum):
-									ys[i][j].append(data[h][j])
+							elif int(groups_cat[s]) in groupB:
+								X.append([0, 1])
+					
+							if int(groups_cat[s]) in np.concatenate( (groupA,groupB), axis=None):
+								#print("hola")
+								for i,h in enumerate(headers):
+									for j in range(clusterNum):
+										ys[i][j].append(data[h][j+1])
 
 		#print(X, ys)
 		X= np.array(X).reshape(len(ys[0][0]),2)
@@ -411,13 +508,101 @@ def groupAnalysis( headers, cols , groups_classification, classification_columns
 					indeces.append(index)
 					
 				 
+	return range(200) #indeces
+"""
+classification colums contains column index for subjects name, folder name, diagnosis label and timepoint
+"""
+def groupAnalysisL( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, length, std, delimiter=",", groupA=[0], groupB=[1], thickness=False):
+	ac_folder=f"dmri.ac/{length}/{std}"
+	indeces=[]
+	groups_cat = {row[classification_columns[0]] :[row[classification_columns[1]],row[classification_columns[2]],row[classification_columns[3]]] for _, row in pd.read_csv(groups_classification, delimiter=delimiter).iterrows()}
+	print(groups_cat)
+
+	motion = {row[0] :row[1] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/outliers_motion.txt", delimiter=' ').iterrows()}
+	subjectsVals=dict()
+	significant_childs=set()
+	for c_i, clusterNum in enumerate(clustersToAnalyze):
+		
+		childs, dads = readTree(clusterNum, f"{subjects_dir}/{target_subject}/{ac_folder}/HierarchicalHistory.csv")
+		#print(childs, dads)
+		order_nodes= pd.read_csv(f"{subjects_dir}/{target_subject}/{ac_folder}/measures/{target_subject}_{target_subject}_c{clusterNum}.csv",delimiter=",", header=0,usecols=[0])
+		#print(order_nodes["Cluster"][0])
+		X=[]
+		
+		for s, elems in groups_cat.items():
+			if s in motion and int(motion[s]) <350:
+				#print(elems)
+				if elems[0] in subjectsVals:
+					ys=subjectsVals[elems[0]]
+				else:
+					ys=[]
+					for j in range(len(headers)):
+						ys.append([])
+						for i in range(clusterNum):
+							ys[j].append([[],[],[],[]])
+					subjectsVals[elems[0]]=ys
+		
+				#print(f'{subjects_dir}/{s}*')
+				gres = glob.glob(f'{subjects_dir}/{s}*')
+				if len(gres) >0:	
+					subject=gres[0].split("/")[-1]
+					measures=f"{subjects_dir}/{subject}/{ac_folder}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
+					#similarities = {row[0] :row[2] for _, row in pd.read_csv(f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/match/{target_subject}_{subject}_c{clusterNum}_hungarian.csv", delimiter=delimiter).iterrows()}
+					if os.path.exists(measures):
+						data = pd.read_csv(measures,delimiter=",", header=0, names=headers, usecols=cols)
+						if len(data[headers[0]])>= clusterNum:
+							for i,h in enumerate(headers):
+								for j in range(clusterNum):
+									ys[i][j][0].append(elems[2])
+									ys[i][j][1].append(data[h][j])
+									ys[i][j][2].append(elems[1])
+									#ys[i][j][3].append(similatities[similarities.keys()[j+1]])
+									
+	labels=["HD", "HD-MCI","HD-D", "MCI-D", "MCI-MCI", "D-D"]
+	labels_dic=dict()
+	labels_dic[(0,0)]=0
+	labels_dic[(0,1)]=1
+	labels_dic[(0,2)]=2
+	labels_dic[(0,3)]=2
+	labels_dic[(1,1)]=4
+	labels_dic[(1,3)]=3
+	labels_dic[(2,2)]=5
+	labels_dic[(2,3)]=5
+	labels_dic[(3,3)]=5
+	patches=[]
+	for i, c in enumerate(colors):
+		patches.append(mpatches.Patch(color=c, label=labels[i]))
+
+
+	for i, m  in enumerate(headers):
+		#print( p_vals)
+		for index in range(clusterNum):
+		
+			plt.figure(index)
+			for s, ys in subjectsVals.items():
+				if(len(ys[i][index][0])>0) :
+					if len(ys[i][index][2] ) >1 and ys[i][index][2][0] <5  and ys[i][index][2][1]<5:
+						col = labels_dic[(ys[i][index][2][0], ys[i][index][2][1])]
+					
+						plt.plot([0,1], ys[i][index][1][0:2], color=colors[col])
+		
+					if len(ys[i][index][2] ) >2 and ys[i][index][2][1] <5  and ys[i][index][2][2]<5:
+						col = labels_dic[(ys[i][index][2][1] , ys[i][index][2][2])]
+					
+						plt.plot([0,1], ys[i][index][1][1:3], color=colors[col])
+		
+			
+			plt.legend(handles=patches, loc="upper right")
+			plt.savefig(f"{subjects_dir}/average/dmri.ac/{length}/{std}/"+str(index)+"_"+headers[i]+"_lines.png")		
+			plt.close()
+	
 	return indeces
 					
 
-def plotAverageMeasures( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, lenght, std, delimiter=",", groups=[[0],[1],[2],[3]], thickness=False):
+def plotAverageMeasures( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, length, std, delimiter=",", groups=[[0],[1],[2],[3]], thickness=False):
 	groups_cat = {row[classification_columns[0]] : row[classification_columns[1]] for _, row in pd.read_csv(groups_classification, delimiter=delimiter).iterrows()}
 	clusterNum=200
-	order_nodes= pd.read_csv(f"{subjects_dir}/{target_subject}/dmri.ac/{lenght}/{std}/measures/{target_subject}_{target_subject}_c{clusterNum}.csv",delimiter=",", header=0,usecols=[0])
+	order_nodes= pd.read_csv(f"{subjects_dir}/{target_subject}/dmri.ac/{length}/{std}/measures/{target_subject}_{target_subject}_c{clusterNum}.csv",delimiter=",", header=0,usecols=[0])
 	measures=[]
 	for g in groups:
 		measures.append([])
@@ -435,9 +620,9 @@ def plotAverageMeasures( headers, cols , groups_classification, classification_c
 		if len(gres)>0:
 			subject=gres[0].split("/")[-1]
 			if thickness:
-				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{lenght}/{std}/measures/surfaceMeasures.csv"
+				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/measures/surfaceMeasures.csv"
 			else:
-				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{lenght}/{std}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
+				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
 			try:
 				data = pd.read_csv(measuresFile,delimiter=",", header=0, names=headers, usecols=cols)
 				#print(measures)
@@ -460,61 +645,80 @@ def plotAverageMeasures( headers, cols , groups_classification, classification_c
 				#print(measures[gi][j][i] ,[gi])
 				plt.violinplot(measures[gi][j][i], [gi],  showmeans=True )
 				#print("dwhola")
-			plt.savefig(f"{subjects_dir}/average/dmri.ac/{lenght}_{std}/"+str(i)+"_"+headers[j]+".png")		
+			plt.savefig(f"{subjects_dir}/average/dmri.ac/{length}/{std}/"+str(i)+"_"+headers[j]+".png")		
 	#plt.show()
-def plotScatter( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, lenght, std, delimiter=",", groups=[[0],[1],[2],[3]], thickness=False):
+def plotScatter( headers, cols , groups_classification, classification_columns, clustersToAnalyze, subjects_dir, target_subject, length, std, delimiter=",", groups=[[0],[1],[2],[3]], thickness=False):
 	groups_cat = {row[classification_columns[0]] : row[classification_columns[1]] for _, row in pd.read_csv(groups_classification, delimiter=delimiter).iterrows()}
+	motion = {row[0] :row[1] for _, row in pd.read_csv("/space/snoke/2/public/vivros/ADDS/FS_STRUCTURE_HCP/outliers_motion.txt", delimiter=' ').iterrows()}
 	clusterNum=200
-	order_nodes= pd.read_csv(f"{subjects_dir}/{target_subject}/dmri.ac/{lenght}/{std}/measures/{target_subject}_{target_subject}_c{clusterNum}.csv",delimiter=",", header=0,usecols=[0])
+	order_nodes= pd.read_csv(f"{subjects_dir}/{target_subject}/dmri.ac/{length}/{std}/measures/{target_subject}_{target_subject}_c{clusterNum}.csv",delimiter=",", header=0,usecols=[0])
 	measures=[]
+	sim=[]
 	for g in groups:
 		measures.append([])
+		sim.append([])
 	for gi,g in enumerate(groups):
 		for a in headers:
 			measures[gi].append([])
+			sim[gi].append([])
 		
 	for gi, g in enumerate(groups):
 		for i in range(clusterNum):
 			for j in range(len(headers)):
 				measures[gi][j].append([])
-	
+				sim[gi][j].append([])
+
+	cv=[]
+	for a in range(clusterNum):
+		cv.append([])	
 	for s in groups_cat.keys():
+	
 		gres = glob.glob(f'{subjects_dir}/{s}*')
-		if len(gres)>0:
+		if len(gres)>0 and (s not in motion or int(motion[s])<250):
 			subject=gres[0].split("/")[-1]
 			if thickness:
-				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{lenght}/{std}/measures/surfaceMeasures.csv"
+				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/measures/surfaceMeasures.csv"
 			else:
-				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{lenght}/{std}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
+				measuresFile=f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/measures/{target_subject}_{subject}_c{clusterNum}.csv"
 			try:
+				similarities = {row[0] :row[2] for _, row in pd.read_csv(f"{subjects_dir}/{subject}/dmri.ac/{length}/{std}/match/{target_subject}_{subject}_c{clusterNum}_hungarian.csv", delimiter=delimiter).iterrows()}
 				data = pd.read_csv(measuresFile,delimiter=",", header=0, names=headers, usecols=cols)
-				#print(measures)
+				print(groups_cat[s])
 				if len(data[headers[0]])>= clusterNum:
 					val  = [ i for i in range(len(groups))  if int(groups_cat[s]) in groups[i]]
-					if len(val)>0:
+					if len(val)>0 :
 						group=val[0]
 						for i,h in enumerate(headers):
-							for j in range(clusterNum):
-								if not math.isnan(data[h][j]):
+							for j, k in enumerate(list(similarities.keys())[1:]):
+								if not math.isnan(data[h][j]) :
 									measures[group][i][j].append(data[h][j])
+									sim[group][i][j].append(similarities[k])
+									cv[j].append(similarities[k])	
+											
 			except Exception as e:
 				print (e, measuresFile)
+	print(groups)
 	for j, h in enumerate(headers):
 		for i in clustersToAnalyze:
 			plt.figure()
 			plt.title(headers[j])
 			#print(headers[j])
-			todo=[]
+			todos=[]
 			for gi, g in enumerate( groups):
+				todo=[]
 				#print(measures[gi][j][i] ,[gi])
-				todo.append	(measures[gi][j][i])
-			plt.boxplot( todo )
+				for a in range(len(sim[gi][j][i])):
+					if np.abs(np.mean(cv[i]) - sim[gi][j][i][a]) < np.std(cv[i]):
+						#print(np.abs(np.mean(cv[i]) - sim[gi][j][i][a]) , np.std(cv[i]))
+						todo.append	(measures[gi][j][i][a])
+				#todos.append(todo)
+				plt.boxplot( todo, positions=[gi] )
 			#print("dwhola")
-			plt.savefig(f"{subjects_dir}/average/dmri.ac/{lenght}_{std}/"+str(i)+"_"+headers[j]+"_scatter.png")		
+			plt.savefig(f"{subjects_dir}/average/dmri.ac/{length}/{std}/"+str(i)+"_"+headers[j]+"_scatter.png")		
 			plt.close()
 	#plt.show()
 
-def GA(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, lenght=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
+def GA(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, length=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
 	measures=['FA','MD','RD','AD']
 	if model == 'DKI':
 		measures=measures+['MK','RK','AK']
@@ -524,40 +728,72 @@ def GA(classificationFile, classification_cols, target_subject, subjects_dir, gr
 	columns = []
 
 
-	#indeces += groupAnalysis(headers=["N"],cols=[1], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, lenght=lenght, std=std)
-	#headers+= ["N"]
-	#columns+=[1]
+	#indeces += groupAnalysis(headers=["N"],cols=[1], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+	headers+= ["N"]
+	columns+=[1]
 
 	for i, m in enumerate(measures):
 		print("mean"+m, " " , 2+4*i, classification_cols)
 		headers += [" mean"+m+" "]
 		columns+=[2+i*4]
-		indeces += groupAnalysis(headers=[" mean"+m+" "],cols=[2+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, lenght=lenght, std=std)
-		#plotAverageMeasures(headers=["mean"+m],cols=[2+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std) #groups=[groupA,groupB], 
-
-	"""
-	for i, m in enumerate(measures):
-		headers += ["   meanCentroid"+m+" "]
-		columns+=[4+i*4]
-		print("meanCentroid"+m, " " , 2+4*i, classification_cols)
-		indeces += groupAnalysis(headers=["   meanCentroid"+m+" "],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, lenght=lenght, std=std)
-		#plotAverageMeasures(headers=["meanCentroid"+m],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std) #groups=[groupA,groupB], 
-		plotScatter(headers=["meanCentroid"+m],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std) #groups=[groupA,groupB], 
+		#indeces += groupAnalysis(headers=[" mean"+m+" "],cols=[2+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+
+	#for i, m in enumerate(measures):
+	#	headers += ["   meanCentroid"+m+" "]
+	#	columns+=[4+i*4]
+	#	print("meanCentroid"+m, " " , 2+4*i, classification_cols)
+		#indeces += groupAnalysis(headers=["   meanCentroid"+m+" "],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+		#plotScatter(headers=["meanCentroid"+m],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std) #groups=[groupA,groupB], 
 	
+	"""
 	if thickness:
 		measures=["curv.star","curv.end","thickness.start", "thickness.end"]
 
 		for i, m in enumerate(measures):
-			indeces += groupAnalysis(headers=[m],cols=[1+i], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, lenght=lenght, std=std, thickness=thickness)
+			indeces += groupAnalysis(headers=[m],cols=[1+i], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std, thickness=thickness)
 
-			#plotAverageMeasures(headers=measures,cols=[0,1,2,3], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std, groups=groups, thickness=True) #groups=[groupA,groupB], 
-			plotScatter(headers=measures,cols=[0,1,2,3], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std, groups=groups, thickness=True) #groups=[groupA,groupB], 
+			#plotAverageMeasures(headers=measures,cols=[0,1,2,3], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=True) #groups=[groupA,groupB], 
+			plotScatter(headers=measures,cols=[0,1,2,3], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=True) #groups=[groupA,groupB], 
 	"""
 	print(columns, headers,indeces)
-	#plotAverageMeasures(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
-	plotScatter(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, lenght=lenght, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
+	#plotAverageMeasures(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
+	indeces=range(200)
+	plotScatter(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
 	
 
+def GAL(classificationFile, classification_cols, target_subject, subjects_dir, groupA, groupB, length=45, std=5,clustersToAnalyze=[200], model="DTI", delimiter="," , groups=[0,1], thickness=False):
+	measures=['FA','MD','RD','AD']
+	if model == 'DKI':
+		measures=measures+['MK','RK','AK']
+
+	indeces = []	
+	headers = []
+	columns = []
+
+
+	#indeces += groupAnalysisL(headers=["N"],cols=[1], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+	#headers+= ["N"]
+	#columns+=[1]
+
+	for i, m in enumerate(measures):
+		print("mean"+m, " " , 2+4*i, classification_cols)
+		headers += [" mean"+m+" "]
+		columns+=[2+i*4]
+		indeces += groupAnalysisL(headers=[" mean"+m+" "],cols=[2+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+
+	"""for i, m in enumerate(measures):
+		headers += ["   meanCentroid"+m+" "]
+		columns+=[4+i*4]
+		print("meanCentroid"+m, " " , 2+4*i, classification_cols)
+		indeces += groupAnalysisL(headers=["   meanCentroid"+m+" "],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=clustersToAnalyze,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, groupA=groupA, groupB=groupB, length=length, std=std)
+		plotScatter(headers=["meanCentroid"+m],cols=[4+i*4], groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std) #groups=[groupA,groupB], 
+	"""
+	#print(columns, headers,indeces)
+	#plotAverageMeasures(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
+	#plotScatter(headers=headers,cols=columns, groups_classification=classificationFile, classification_columns=classification_cols,clustersToAnalyze=indeces,target_subject=target_subject,subjects_dir=subjects_dir, delimiter=delimiter, length=length, std=std, groups=groups, thickness=False) #groups=[groupA,groupB], 
+	
+
+
 
 
 cta=[200]
@@ -574,10 +810,11 @@ parser.add_argument('-cc','--classification_columns',required=False,  help='clas
 parser.add_argument('-cta','--clusters_to_analyze',required=False,  help='Clusters to analyze')
 parser.add_argument('-ts','--target_subject',required=False,  help='target subject')
 parser.add_argument('-s','--subjects_folder', required=False, help='subject folder')
+parser.add_argument('-sn','--subjectsNames', required=False, help='subject names')
 parser.add_argument('-ga','--group_a',required=False,  help='group a')
 parser.add_argument('-gb','--group_b',required=False,  help='group b')
 parser.add_argument('-d','--delimiter',required=False,  help='delimiter for classification file')
-parser.add_argument('-l','--lenght',required=False,  help='minimum lenght')
+parser.add_argument('-l','--length',required=False,  help='minimum length')
 parser.add_argument('-std','--std',required=False,  help='standard deviation of clusters')
 parser.add_argument('-pt','--plot_groups',required=False,  help='groups for plot')
 parser.add_argument('-t','--thickness',required=False,  help='analyze thickness')
@@ -585,37 +822,32 @@ parser.add_argument('-t','--thickness',required=False,  help='analyze thickness'
 
 args = parser.parse_args()
 
-if args.function == "GA":
+groups=[[7],[0],[1],[2,3]]
+print(groups)
+if args.function == "GA" or args.function == "GAL":
 	columns=  np.asarray(args.classification_columns.split(":"), dtype=np.int)
 	print(columns, columns[0])
 	group_a = np.asarray(args.group_a.split(":"), dtype=np.int)
 	group_b =np.asarray(args.group_b.split(":"), dtype=np.int)
 	clusters_to_analyze= np.asarray(args.clusters_to_analyze.split(":"), dtype=np.int)
-		
-	#groups=np.asarray(args.plot_groups.split(":"), dtype=np.int)
-	groups=[[999999999,0],[1],[2]]
-	print(groups)
 
-	eval(args.function)(args.classification_file, columns, args.target_subject, args.subjects_folder,group_a, group_b,args.lenght, args.std, clusters_to_analyze, args.model, args.delimiter,groups, args.thickness )
-elif args.function == "thicknessPerStructure" or args.function == "connectivityGraph":
+	eval(args.function)(args.classification_file, columns, args.target_subject, args.subjects_folder,group_a, group_b,args.length, args.std, clusters_to_analyze, args.model, args.delimiter,groups, args.thickness )
+elif args.function == "thicknessPerStructure" or args.function == "connectivityGraph" or args.function == "diffAlongTheLine":
 	columns=  np.asarray(args.classification_columns.split(":"), dtype=np.int)
 	print(columns, columns[0])
 	group_a = np.asarray(args.group_a.split(":"), dtype=np.int)
 	group_b =np.asarray(args.group_b.split(":"), dtype=np.int)
 	clusters_to_analyze= np.asarray(args.clusters_to_analyze.split(":"), dtype=np.int)
-		
-	#groups=np.asarray(args.plot_groups.split(":"), dtype=np.int)
-	groups=[[999999999,0],[1],[2]]
-	print(groups)
-
-	eval(args.function)(args.classification_file, columns, args.target_subject, args.subjects_folder,group_a, group_b,args.lenght, args.std, clusters_to_analyze, args.model, args.delimiter,groups )
 
+	eval(args.function)(args.classification_file, columns, args.target_subject, args.subjects_folder,group_a, group_b,args.length, args.std, clusters_to_analyze, args.model, args.delimiter,groups )
 
 elif args.function == "averageCorrespondingClusters":
 	corresponding_folders =  args.corresponding_file_list.replace("]","").replace("[","").split(",")
 	images_folders =  args.imagesFolder.replace("]","").replace("[","").split(",")
+	subjectsNames=  args.subjectsNames.replace("]","").replace("[","").split(",")
 	cluster_indeces =  np.asarray(args.clusterIndeces.replace("]","").replace("[","").split(","),dtype=np.int)
-	eval(args.function)(corresponding_folders, images_folders, args.outputFolder,cluster_indeces)
+	columns=  np.asarray(args.classification_columns.split(":"), dtype=np.int)
+	eval(args.function)(corresponding_folders, images_folders, args.outputFolder,cluster_indeces,subjectsNames , args.classification_file, columns,groups, args.delimiter)
 
 
 """ts="BAKP64e_nrecon-trio3_vb17"
diff --git a/anatomicuts/dmri_ac.sh b/anatomicuts/dmri_ac.sh
index 09094a3b607..a79474bd776 100644
--- a/anatomicuts/dmri_ac.sh
+++ b/anatomicuts/dmri_ac.sh
@@ -40,7 +40,9 @@ function forAll()
 	echo ${cluster_call}
 	
 	cd ${SUBJECTS_DIR}                                                
+	
 	for s in */;       
+	#for s in `ls -dir [a-zA-Z]*`;       
 	do
 		subject=${s//[\/]/}                                                                                                    
 		echo ${subject}
@@ -259,12 +261,12 @@ function preGA()
 	std=$4
 	bb=$5	
 
-	Clean ${subject} ${targetSubject} ${lenght} ${std} ${bb}
-	Hungarian ${subject} ${targetSubject} ${lenght} ${std} ${bb}
-	Measures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
-	ToAnat  ${subject} ${targetSubject} ${lenght} ${std} ${bb}
+	#Clean ${subject} ${targetSubject} ${lenght} ${std} ${bb}
+	#Hungarian ${subject} ${targetSubject} ${lenght} ${std} ${bb}
+	#Measures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
+	#ToAnat  ${subject} ${targetSubject} ${lenght} ${std} ${bb}
 	SurfaceMeasures ${subject} ${targetSubject} ${lenght} ${std} ${bb}
-	ToTarget ${subject} ${targetSubject} ${lenght} ${std} ${bb}
+	#ToTarget ${subject} ${targetSubject} ${lenght} ${std} ${bb}
 	
 }
 function GA()
@@ -278,8 +280,12 @@ function GA()
 	groupB=$7
 	groups=$8
 	thickness=$9
+#	mkdir -p ${ODMRI_DIR}/average/dmri.ac/${2}/${3}/
+#	rm 	 ${ODMRI_DIR}/average/dmri.ac/${2}/${3}/*
+#${FREESURFER_HOME}/bin/anatomiCutsUtils -f GAL -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} 
 #${FREESURFER_HOME}/bin/anatomiCutsUtils -f GA -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} 
-${FREESURFER_HOME}/bin/anatomiCutsUtils -f thicknessPerStructure -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness} 
+${FREESURFER_HOME}/bin/anatomiCutsUtils -f diffAlongTheLine -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} 
+#${FREESURFER_HOME}/bin/anatomiCutsUtils -f thicknessPerStructure -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness} 
 #${FREESURFER_HOME}/bin/anatomiCutsUtils -f connectivityGraph -m "DTI" -cf "${labels_file}" -cc ${labels_cols} -cta 50:100:150:200 -ts ${targetSubject} -s ${ODMRI_DIR} -d "," -ga $groupA -gb $groupB -l ${lenght} -std ${std} -pt ${groups} -t ${thickness} 
 
 #${FREESURFER_HOME}/bin/anatomiCutsUtils -f GA -m "DKI" -cf "/space/snoke/1/public/vivros/data/demos_fullID.csv" -cc 0:6 -cta 200 -ts ${targetSubject} -s ${ODMRI_DIR} -d " " -ga 3 -gb 1 -l ${lenght} -std ${std} 
@@ -323,6 +329,7 @@ function Hungarian()
 		ci=${ODMRI_DIR}/${targetSubject}/dmri.ac/${lenght}/${std}
 		cj=${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}
 		mkdir -p ${cj}/match/	
+		echo ${HungarianBin} -s1 ${si} -s2 ${sj} -h1 ${ci}/ -h2 ${cj}/ -o ${cj}/match/${targetSubject}_${subject}_c${c}_hungarian.csv  -labels -hungarian -c ${c} ${bb}
 		${HungarianBin} -s1 ${si} -s2 ${sj} -h1 ${ci}/ -h2 ${cj}/ -o ${cj}/match/${targetSubject}_${subject}_c${c}_hungarian.csv  -labels -hungarian -c ${c} ${bb}
 	done
 }
@@ -341,11 +348,12 @@ function Measures()
 
 		if [[ -e   ${SUBJECTS_DIR}/${subject}/dmri/DKI/dki_RD.nii.gz ]] ;
 		then 
-			string="${stats_ac_bin} -i ${anatomicuts}/  -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 7 FA ${diff}/DKI/dki_FA.nii   MD ${diff}/DKI/dki_MD.nii   RD ${diff}/DKI/dki_RD.nii   AD ${diff}/DKI/dki_AD.nii   MK ${diff}/DKI/dki_MK.nii   RK ${diff}/DKI/dki_RK.nii   AK ${diff}/DKI/dki_AK.nii   -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
+			string="${stats_ac_bin} -i ${anatomicuts}/  -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 7 FA ${diff}/DKI/dki_FA.nii.gz   MD ${diff}/DKI/dki_MD.nii.gz   RD ${diff}/DKI/dki_RD.nii.gz   AD ${diff}/DKI/dki_AD.nii.gz   MK ${diff}/DKI/dki_MK.nii.gz   RK ${diff}/DKI/dki_RK.nii.gz   AK ${diff}/DKI/dki_AK.nii.gz   -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
 		else
-			string="${stats_ac_bin} -i ${anatomicuts}/  -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 4 FA ${diff}/DTI/dti_FA.nii.gz   MD ${diff}/DTI/dti_MD.nii   RD ${diff}/DTI/dti_RD.nii   AD ${diff}/DTI/dti_AD.nii  -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
+			string="${stats_ac_bin} -i ${anatomicuts}/  -n ${c} -c ${anatomicuts}/match/${targetSubject}_${subject}_c${c}_hungarian.csv -m 4 FA ${diff}/DTI/dti_FA.nii.gz.gz   MD ${diff}/DTI/dti_MD.nii.gz   RD ${diff}/DTI/dti_RD.nii.gz   AD ${diff}/DTI/dti_AD.nii.gz  -o ${anatomicuts}/measures/${targetSubject}_${subject}_c${c}.csv"
 		fi
 
+		echo ${string}
 		${string}
 	done
 }
@@ -363,8 +371,9 @@ function SurfaceMeasures()
 	annot=${SUBJECTS_DIR}/${subject}/label/
 	mri=${SUBJECTS_DIR}/${subject}/mri/
 	
+	diff=${ODMRI_DIR}/${subject}/dmri/
 
- 	${FREESURFER_HOME}/bin/dmri_extractSurfaceMeasurements -i ${anatomicuts}/*trk -sl ${surf}/lh.pial -tl ${surf}/lh.thickness -cl ${surf}/lh.curv.pial -sr ${surf}/rh.pial -tr ${surf}/rh.thickness -cr ${surf}/rh.curv.pial -rid ${mri}/brain.nii.gz -ria ${mri}/brain.nii.gz  -al ${annot}/lh.aparc.annot -ar ${annot}/rh.aparc.annot -o ${anatomicutsdiff}/measures/	-p ${anatomicutsdiff}/match/${targetSubject}_${subject}_c200_hungarian.csv 
+ 	${FREESURFER_HOME}/bin/dmri_extractSurfaceMeasurements -i ${anatomicuts}/*.trk -sl ${surf}/lh.pial -tl ${surf}/lh.thickness -cl ${surf}/lh.curv.pial -sr ${surf}/rh.pial -tr ${surf}/rh.thickness -cr ${surf}/rh.curv.pial -rid ${mri}/brain.nii.gz -ria ${mri}/brain.nii.gz  -al ${annot}/lh.aparc.annot -ar ${annot}/rh.aparc.annot -o ${anatomicutsdiff}/measures/	-p ${anatomicutsdiff}/match/${targetSubject}_${subject}_c200_hungarian.csv -fa 7 FA ${diff}/DKI/dki_FA.nii.gz   MD ${diff}/DKI/dki_MD.nii.gz   RD ${diff}/DKI/dki_RD.nii.gz   AD ${diff}/DKI/dki_AD.nii.gz   MK ${diff}/DKI/dki_MK.nii.gz   RK ${diff}/DKI/dki_RK.nii.gz   AK ${diff}/DKI/dki_AK.nii.gz   
 
 }
 function ToAnat()
@@ -392,6 +401,7 @@ function ToAnat()
 	mkdir -p ${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}/${common}/
         
 	common_clustering=${ODMRI_DIR}/${subject}/dmri.ac/${lenght}/${std}/${common}/
+	mkdir ${diff}/xfms/
 
 	flirt -in ${diff}/DTI/dti_FA.nii.gz -ref ${mri}/brain.nii.gz  -omat ${diff}/xfms/fa2brain.mat
 
@@ -464,7 +474,13 @@ function average()
 	targetSubject=$1	
 	lenght=$2
 	std=$3
-
+	labels_file=$4
+	labels_cols=$5
+	groupA=$6
+	groupB=$7
+	groups=$8
+	thickness=$9
+#	
 	mkdir -p ${ODMRI_DIR}/average/dmri.ac/${lenght}/${std}/images
         #correspondences="["
         #imagesFolder="["
@@ -480,13 +496,17 @@ function average()
                         if [ ${#correspondences} -ge 3 ]; then 
                                 correspondences=${correspondences},
                                 imagesFolder=${imagesFolder},
+                                subjectsNames=${subjectsNames},
                         fi
                         correspondences=${correspondences}${ODMRI_DIR}/${s}/dmri.ac/${lenght}/${std}/match/${s2}_${s}_c200_hungarian.csv   
                         imagesFolder=${imagesFolder}${ODMRI_DIR}/${s}/dmri.ac/${lenght}/${std}/to${targetSubject}/images/
+			subjectsNames=${subjectsNames}${s}
+			
                 fi
         done    
         correspondences=${correspondences}
         imagesFolder=${imagesFolder}
+	subjectsNames=${subjectsNames}	
 
         echo $correspondences
         echo $imagesFolder
@@ -494,9 +514,9 @@ function average()
         echo $clusterIndeces
         mkdir -p ${outputFolder}
         #correspondences, imagesFolder, outputFolder,  clusterIndeces   
-        cd /space/erebus/2/users/vsiless/code/freesurfer/anatomicuts/ 
+        cd ${SUBJECTS_DIR} 
         #pbsubmit -n 1 -c "python3 -c \"import anatomiCutsUtils;  anatomiCutsUtils.averageCorrespondingClusters($correspondences, $imagesFolder, $outputFolder,$clusterIndeces) \" "
-	anatomiCutsUtils -f averageCorrespondingClusters -co ${correspondences} -if ${imagesFolder} -of ${outputFolder} -in ${clusterIndeces} 
+	${FREESURFER_HOME}/bin/anatomiCutsUtils -f averageCorrespondingClusters -co ${correspondences} -if ${imagesFolder} -of ${outputFolder} -in ${clusterIndeces} -sn ${subjectsNames} -cf "${labels_file}" -cc ${labels_cols} -d "," -pt ${groups} 
         #python3 -c "import anatomiCutsUtils;  anatomiCutsUtils.averageCorrespondingClusters(${correspondences}, ${imagesFolder}, ${outputFolder},${clusterIndeces}) "
         
 }
diff --git a/anatomicuts/dmri_extractSurfaceMeasurements.cxx b/anatomicuts/dmri_extractSurfaceMeasurements.cxx
index 5799f7b04c8..ddf873162e5 100644
--- a/anatomicuts/dmri_extractSurfaceMeasurements.cxx
+++ b/anatomicuts/dmri_extractSurfaceMeasurements.cxx
@@ -226,7 +226,7 @@ int main(int narg, char* arg[])
 	
 	ClusterToolsType::Pointer clusterTools = ClusterToolsType::New();
 	clusterTools->GetPolyDatas(TRKFiles, &polydatas, ref_Image.at(0));
-	meshes = clusterTools->PolydataToMesh(polydatas);
+	meshes = clusterTools->FixSampleClusters(polydatas,10);
 	
 	//Loading the surface for each hemisphere and metric
 	//Left Curvature
@@ -322,8 +322,8 @@ int main(int narg, char* arg[])
 		getline ( file, value, ',' ); 
 		long long v2 = atoll(value.c_str());
 		correspondences.push_back(v2);		
-		std::cout << " v2 " <<  v2 << std::endl;
 	} 	
+	std::cout << " Files " <<  meshes->size() << std::endl;
 	// Cycling through the TRK files
 	for(int i = 0; i < meshes->size(); i++)
 	{ 
@@ -343,9 +343,12 @@ int main(int narg, char* arg[])
 		if (FA_FOUND)
 		{
 			for (int a = 0; a < image_fileNames.size(); a++)
-				oFile << ", mean" << image_fileNames.at(a) << ", stde" << image_fileNames.at(a);
+				oFile << ",mean" << image_fileNames.at(a) << ",stde" << image_fileNames.at(a);
+			for (int a = 0; a < image_fileNames.size(); a++)
+				for (int b =0; b<10; b++)
+					oFile << ","<< b << "_"<< image_fileNames.at(a);
 		} 
-		oFile << endl;
+		oFile << ","<<endl;
 
 		// Initialization of a new stream for every TRK files
 		
@@ -357,22 +360,25 @@ int main(int narg, char* arg[])
 		for (; inputCellIt != input->GetCells()->End(); ++inputCellIt, ++counter)
 		{
 			vector<float> meanFA;
+			vector<vector<float>> allFA;
 			vector<float> stdeFA;
 		
 			// If there are image files, then find the mean and stde of FA	
 			if (FA_FOUND)
-			{
+			{	
+				
 				for (int p = 0; p < volumes.size(); p++)
 				{
+					allFA.push_back(vector<float>());
 					// Creating streamline variable and finding first point
 					CellType::PointIdIterator it = inputCellIt.Value()->PointIdsBegin();
 					input->GetPoint(*it, &firstPt);
-
+					input->GetPoint(*inputCellIt.Value()->PointIdsEnd(), &lastPt);
+					float val = (firstPt[0] -lastPt[0])*(firstPt[1]-lastPt[1])*(firstPt[2]-lastPt[2]);
+					
 					// Getting the FA value at all points
 					vector<float> FA_values;
 					ImageType::IndexType index;
-					if (volumes.at(p)->TransformPhysicalPointToIndex(firstPt, index))
-						FA_values.push_back(volumes.at(p)->GetPixel(index));
 				
 					// Cycling through the points in the stream
 					for (; it != inputCellIt.Value()->PointIdsEnd(); it++)
@@ -381,7 +387,13 @@ int main(int narg, char* arg[])
 
 						// If FA options is used, then add the FA values to the vector	
 						if (volumes.at(p)->TransformPhysicalPointToIndex(lastPt, index))
+						{
 							FA_values.push_back(volumes.at(p)->GetPixel(index));	
+							if( val>0)
+								allFA[p].push_back(volumes.at(p)->GetPixel(index));	
+							else
+								allFA[p].insert(allFA[p].begin(),volumes.at(p)->GetPixel(index));	
+						}	
 					}
 			
 					// Calculating the MeanFA and stdeFA
@@ -437,81 +449,123 @@ int main(int narg, char* arg[])
 
 
 			// Outputting values to the file
+			std::cout << ID1 << " " <<ID2<< std::endl;
 			oFile << "StreamLine" << counter << ", ";
 			int structure;
-			if (ID1 == Left_ID1)
+			if (ID1 >=0)
 			{
-		 		CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID1].annotation, &structure);
-				std::cout << ID1 << " " <<   surfCL->vertices[ID1].annotation << " " <<structure<< std::endl;
+				if (ID1 == Left_ID1)
+				{
+					CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID1].annotation, &structure);
+					std::cout <<"L" << ID1 << " " <<   surfCL->vertices[ID1].annotation << " " <<structure<< std::endl;
+				}
+				else
+				{
+					CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID1].annotation, &structure);
+					std::cout << "R" << ID1 << " " <<   surfCR->vertices[ID1].annotation << " " <<structure<< std::endl;
+				}
+				oFile << structure << ",";
 			}
 			else
 			{
-		 		CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID1].annotation, &structure);
-				std::cout << ID1 << " " <<   surfCR->vertices[ID1].annotation << " " <<structure<< std::endl;
+				oFile << ",";
 			}
-			oFile << structure << ",";
-			if (ID2 == Left_ID2)
+			if(ID2>=0)
 			{
-		 		CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID2].annotation, &structure);
-				std::cout << ID2 << " " <<   surfCL->vertices[ID2].annotation << " " <<structure<< std::endl;
+				if (ID2 == Left_ID2)
+				{
+					CTABfindAnnotation(surfCL->ct , surfCL->vertices[ID2].annotation, &structure);
+					std::cout << "L"<< ID2 << " " <<   surfCL->vertices[ID2].annotation << " " <<structure<< std::endl;
+				}
+				else
+				{
+					CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID2].annotation, &structure);
+					std::cout <<"R"<< ID2 << " " <<   surfCR->vertices[ID2].annotation << " " <<structure<< std::endl;
+				}
+
+				oFile << structure << ",";
 			}
 			else
 			{
-		 		CTABfindAnnotation(surfCR->ct , surfCR->vertices[ID2].annotation, &structure);
-				std::cout << ID2 << " " <<   surfCR->vertices[ID2].annotation << " " <<structure<< std::endl;
+				oFile << ",";
 			}
-
-			oFile << structure << ",";
-
-			if (ID1 == Left_ID1)
+			if (ID1 >=0)
 			{
-				oFile << surfCL->vertices[ID1].curv << ",";
-				values[0]+= surfCL->vertices[ID1].curv ;
+				if (ID1 == Left_ID1)
+				{
+					oFile << surfCL->vertices[ID1].curv << ",";
+					values[0]+= surfCL->vertices[ID1].curv ;
+				}
+				else
+				{
+					oFile << surfCR->vertices[ID1].curv << ",";
+					values[0]+= surfCR->vertices[ID1].curv ;
+				}
 			}
 			else
 			{
-				oFile << surfCR->vertices[ID1].curv << ",";
-				values[0]+= surfCR->vertices[ID1].curv ;
+				oFile << ",";
 			}
-	
-			if (ID2 == Left_ID2)
+			if (ID2 >=0)
 			{
-				oFile << surfCL->vertices[ID2].curv << ",";
-				values[1]+= surfCL->vertices[ID2].curv ;
-	
+				if (ID2 == Left_ID2)
+				{
+					oFile << surfCL->vertices[ID2].curv << ",";
+					values[1]+= surfCL->vertices[ID2].curv ;
+
+				}
+				else
+				{
+					oFile << surfCR->vertices[ID2].curv << ",";
+					values[1]+= surfCR->vertices[ID2].curv ;
+				}
 			}
 			else
 			{
-				oFile << surfCR->vertices[ID2].curv << ",";
-				values[1]+= surfCR->vertices[ID2].curv ;
-			}
-			      
-			if (ID1 == Left_ID1)
+				oFile << ",";
+			}	
+
+			if (ID1 >=0)
 			{
-				oFile << surfTL->vertices[ID1].curv << ","; 
-				values[2]+= surfTL->vertices[ID1].curv ;
+				if (ID1 == Left_ID1)
+				{
+					oFile << surfTL->vertices[ID1].curv << ","; 
+					values[2]+= surfTL->vertices[ID1].curv ;
+				}
+				else
+				{
+					oFile << surfTR->vertices[ID1].curv << ",";
+					values[2]+= surfTR->vertices[ID1].curv ;
+				}
 			}
 			else
 			{
-				oFile << surfTR->vertices[ID1].curv << ",";
-				values[2]+= surfTR->vertices[ID1].curv ;
+				oFile << ",";
 			}
-
-
-			if (ID2 == Left_ID2)
+			if (ID2 >=0)
 			{
-				oFile << surfTL->vertices[ID2].curv;
-				values[3]+= surfTL->vertices[ID2].curv ;
+				if (ID2 == Left_ID2)
+				{
+					oFile << surfTL->vertices[ID2].curv <<",";
+					values[3]+= surfTL->vertices[ID2].curv ;
+				}		
+				else
+				{
+					oFile << surfTR->vertices[ID2].curv <<",";
+					values[3]+= surfTR->vertices[ID2].curv ;
+				}
 			}		
 			else
 			{
-				oFile << surfTR->vertices[ID2].curv;
-				values[3]+= surfTR->vertices[ID2].curv ;
-			}		
+				oFile << ",";
+			}
 			if (FA_FOUND)
 			{
-				for (int m = 0; m < stdeFA.size(); m++)
-					oFile << "," << meanFA.at(m) << "," << stdeFA.at(m);
+				for (int m = 0; m < volumes.size(); m++)
+					oFile << meanFA.at(m) << "," << stdeFA.at(m)<<",";
+				for (int m = 0; m < volumes.size(); m++)
+					for (int b = 0; b < allFA[m].size(); b++)
+						oFile <<  allFA[m][b]<<",";
 			}
 
 			oFile << endl;
diff --git a/anatomicuts/dmri_neighboringRegions.cxx b/anatomicuts/dmri_neighboringRegions.cxx
index b42e2698bd6..d99e1f409b6 100644
--- a/anatomicuts/dmri_neighboringRegions.cxx
+++ b/anatomicuts/dmri_neighboringRegions.cxx
@@ -14,6 +14,6 @@ using namespace std;
 int main(int argc, char *argv[])
 {
 	// TODO
-
+	std::cout << " chau " << std::endl;
 	return 0;
 }
diff --git a/freeview/LayerSurface.cpp b/freeview/LayerSurface.cpp
index 02d83207222..a6a4fc45326 100644
--- a/freeview/LayerSurface.cpp
+++ b/freeview/LayerSurface.cpp
@@ -809,7 +809,12 @@ void LayerSurface::InitializeActors()
 
     vtkSmartPointer<vtkPolyDataMapper> mapper3 = vtkSmartPointer<vtkPolyDataMapper>::New();
     m_vertexPoly2D[i] = vtkSmartPointer<vtkPolyData>::New();
+#if VTK_MAJOR_VERSION > 5
     mapper3->SetInputData(m_vertexPoly2D[i]);
+#else
+mapper3->SetInput(m_vertexPoly2D[i]);
+#endif 
+
     mapper3->ScalarVisibilityOff();
     m_vertexActor2D[i]->SetMapper(mapper3);
     m_vertexActor2D[i]->SetProperty( m_vertexActor2D[i]->MakeProperty() );
diff --git a/mris_make_surfaces/mris_make_surfaces.cpp b/mris_make_surfaces/mris_make_surfaces.cpp
index 06ae6d82411..7685328c005 100644
--- a/mris_make_surfaces/mris_make_surfaces.cpp
+++ b/mris_make_surfaces/mris_make_surfaces.cpp
@@ -1650,8 +1650,8 @@ int main(int argc, char *argv[])
 		MRIS_MultimodalRefinement* refine = new MRIS_MultimodalRefinement();
 		MRI* whiteMR= MRIcopy(mri_T1_pial,NULL);
 		MRI* vesselMR= MRIcopy(mri_T1_pial,NULL);
-		refine->SegmentWM(mri_T1_pial,mri_flair, whiteMR);
-		refine->SegmentVessel(mri_T1_pial,mri_flair, vesselMR);
+		refine->SegmentWM(mri_T1_pial,mri_flair, whiteMR, contrast_type);
+		refine->SegmentVessel(mri_T1_pial,mri_flair, vesselMR, contrast_type);
 
 		refine->SetStep(.4);
 		refine->SetNumberOfSteps(8);
diff --git a/mris_make_surfaces/mris_place_surface.cpp b/mris_make_surfaces/mris_place_surface.cpp
index 63f6fad66e4..57077703705 100644
--- a/mris_make_surfaces/mris_place_surface.cpp
+++ b/mris_make_surfaces/mris_place_surface.cpp
@@ -222,6 +222,7 @@ double shrinkThresh = -1;
 int mm_contrast_type = -1;
 char *mmvolpath = NULL;
 MRI *mmvol = NULL;
+std::map<int, MRI*> mmvols;
 float T2_min_inside = 110 ;
 float T2_max_inside = 300 ;
 double T2_min_outside = 130;
@@ -234,6 +235,8 @@ double Ghisto_left_outside_peak_pct = 0.5;
 double Ghisto_right_outside_peak_pct = 0.5;
 int n_averages=0;
 int UseMMRefine = 0;
+float MMRefineMinPGrey = 20;
+int UseMMRefineWeights = 0;
 AutoDetGWStats adgws;
 char *coversegpath = NULL;
 MRI *mri_cover_seg = NULL;
@@ -485,8 +488,11 @@ int main(int argc, char **argv)
     printf("Reading in multimodal volume %s\n",mmvolpath);
     mmvol = MRIread(mmvolpath);
     if(mmvol==NULL) exit(1);
-    // should check that the dimensions, resolution are the same
-    parms.l_intensity = 0.000;
+  }
+ // should check that the dimensions, resolution are the same
+ if (mmvolpath || UseMMRefine || UseMMRefineWeights ||mmvols.size()>0){
+  	std::cout << " using multi modal weights "<< std::endl;
+	  parms.l_intensity = 0.000;
     parms.l_location  = 0.500;
     parms.l_repulse   = 0.025;
     parms.l_nspring   = 1.000;
@@ -586,7 +592,7 @@ int main(int argc, char **argv)
     parms.n_averages = n_averages ;
 
     INTEGRATION_PARMS_copy(&old_parms, &parms) ;
-    if(mmvol == NULL){
+    if(mmvol == NULL &&  mmvols.size()==0 && !UseMMRefine){
       // Compute the target intensity value (l_intensity)
       printf("Computing target border values \n");
       // The outputs are set in each vertex structure:
@@ -635,23 +641,43 @@ int main(int argc, char **argv)
 					       T2_min_outside, T2_max_outside, 
 					       max_outward_dist,
 					       Ghisto_left_inside_peak_pct, Ghisto_right_inside_peak_pct, 
-					       Ghisto_left_outside_peak_pct, Ghisto_right_outside_peak_pct, 
+					     Ghisto_left_outside_peak_pct, Ghisto_right_outside_peak_pct, 
 					       wm_weight, pial_sigma, invol) ;
       }
       else{
-	printf("UseMMRefine\n");
+	std::cout<<"UseMMRefine, num vols"  << mmvols.size() << std::endl;
 	MRIS_MultimodalRefinement* refine = new MRIS_MultimodalRefinement();
 	MRI* whiteMR  = MRIcopy(invol,NULL);
 	MRI* vesselMR = MRIcopy(invol,NULL);
-	refine->SegmentWM(invol,mmvol, whiteMR);
-	refine->SegmentVessel(invol,mmvol, vesselMR);
+	if(mmvols.count(CONTRAST_T2)==0 && mmvols.count(CONTRAST_FLAIR)==0)
+	{
+		refine->SegmentWM(mmvols[CONTRAST_T1],seg, whiteMR, -1);
+		refine->SegmentVessel(mmvols[CONTRAST_T1],seg, vesselMR, -1);
+	}
+	else if (mmvols.count(CONTRAST_T1) ==0)
+	{
+		int contrast = mmvols.count(CONTRAST_T2)>0?CONTRAST_T2:CONTRAST_FLAIR;
+		refine->SegmentWM(invol,mmvols[contrast], whiteMR, contrast);
+		refine->SegmentVessel(invol,mmvols[contrast], vesselMR, contrast);
+	}
+	else
+	{
+		int contrast = mmvols.count(CONTRAST_T2)>0?CONTRAST_T2:CONTRAST_FLAIR;
+		refine->SegmentWM(mmvols[CONTRAST_T1],mmvols[contrast], whiteMR, contrast);
+		refine->SegmentVessel(mmvols[CONTRAST_T1],mmvols[contrast], vesselMR, contrast);
+	}
 	refine->SetStep(.4);
 	refine->SetNumberOfSteps(12);
 	refine->SetGradientSigma(.3);
 	refine->SetSegmentation(seg);
-	refine->FindMaximumGradient(mm_contrast_type == CONTRAST_T2);
-	refine->addImage(invol);
-	refine->addImage(mmvol);
+	refine->SetMinPGrey(MMRefineMinPGrey);
+	//refine->FindMaximumGradient(mm_contrast_type == CONTRAST_T2);
+	//refine->addImage(invol);
+	for( auto& x:mmvols)
+	{
+		std::cout << x.first << std::endl;
+		refine->addImage(x.second);
+	}
 	refine->SetWhiteMR(whiteMR);
 	refine->SetVesselMR(vesselMR);
 	refine->getTarget(surf); //, debugVertex);
@@ -844,8 +870,43 @@ static int parse_commandline(int argc, char **argv) {
       coversegpath = pargv[0];
       nargsused = 1;
     } 
-    else if(!strcasecmp(option, "--mm-refine"))   UseMMRefine = 1;
-    else if(!strcmp(option, "--i")){
+    else if(!strcasecmp(option, "--mm-min-p-grey")){
+	MMRefineMinPGrey=atof(pargv[0]);
+	nargsused=1;
+    }else if(!strcasecmp(option, "--mm-weights")){
+		
+	UseMMRefineWeights=1;
+    } else if(!strcasecmp(option, "--mm-refine")){
+	    UseMMRefine = 1;
+	    int numImages  =atoi(pargv[0]);
+	    
+	for(int i=0;i<numImages;i++)
+	    {		
+		if(!stricmp(pargv[i*2+1],"t2"))
+		{
+			std::cout << CONTRAST_T2 << " T2 " << pargv[i*2+2]<< std::endl;	
+		   	mmvols[CONTRAST_T2]=MRIread(pargv[i*2+2]);
+		}
+		else if(!stricmp(pargv[i*2+1],"flair"))
+		{
+			std::cout << CONTRAST_FLAIR << " FLAIR " << pargv[i*2+2]<< std::endl;	
+			 mmvols[ CONTRAST_FLAIR]=MRIread(pargv[i*2+2]);
+		}
+		else if(!stricmp(pargv[i*2+1],"t1"))
+		{
+			std::cout << CONTRAST_T1 << " T1 " << pargv[i*2+2]<< std::endl;	
+			 mmvols[ CONTRAST_T1]=MRIread(pargv[i*2+2]);
+		}
+		else 
+		{
+			printf("ERROR: mmvol must be either t2 or flair weighted\n");
+			exit(1);
+		}
+	    }
+      surftype = GRAY_CSF;
+      nargsused = numImages*2+1;
+
+    }else if(!strcmp(option, "--i")){
       if(nargc < 1) CMDargNErr(option,1);
       insurfpath = pargv[0];
       parms.l_tspring = 0.3;
diff --git a/resurf/Code/mris_multimodal_refinement.h b/resurf/Code/mris_multimodal_refinement.h
index 16fcf18e24c..5aabfa86449 100644
--- a/resurf/Code/mris_multimodal_refinement.h
+++ b/resurf/Code/mris_multimodal_refinement.h
@@ -30,10 +30,12 @@ class MRIS_MultimodalRefinement {
 		void SetSphere(MRIS* sphere){ this->sphere=sphere;}
 		void getNormal(MRIS* surf,MRI* image, int vtxno, double* x, double* y, double* z);
 		std::vector<std::pair<float, float>> GetMeanAndVariance(MRIS* surf, MRI* image);
-		void  SegmentWM(MRI* t1, MRI* t2, MRI* output);
-		void  SegmentVessel(MRI* t1, MRI* t2, MRI* output);
+		void  SegmentWM(MRI* t1, MRI* t2, MRI* output, int contrast_type);
+		void  SegmentVessel(MRI* t1, MRI* t2, MRI* output, int contrast_type);
 		void GeneratePlanes();
 		float DistanceToMidline(double x, double y, double z);
+		void SetMinPGrey(float p){ this->minPGrey = p;}
+
 	private:
 		MRIS* white;
 		MRIS* sphere;
@@ -55,6 +57,7 @@ class MRIS_MultimodalRefinement {
 		double center[3];	
 		double normal[3];
 		float GetVentricleIntensity();	
+		float minPGrey=20;
 		
 };
 void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
@@ -63,6 +66,8 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 	std::cout << " number of steps "<< this->numberOfSteps <<std::endl;
 	std::cout << " step "<< step << std::endl;
 	std::cout <<" gradient sigma " << this->gradientSigma << std::endl;
+	std::cout <<" Min P Grey  " << this->minPGrey << std::endl;
+	std::cout << " min p " << 1.0 / pow(10,this->minPGrey) << " " << 1.0 / pow(10,this->minPGrey/3)<< std::endl;
 	//std::map<std::tuple<int,int,int>, int> priors =getPs(surf);
 			
   	MRIS_HASH_TABLE *mht ;
@@ -72,7 +77,7 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 	//this->GeneratePlanes();
 	std::vector<float> intensities(this->numberOfSteps*5+1,0);
 	std::vector<float> magnitudes(this->numberOfSteps*5+1,0);
-	std::vector<float> ps(this->numberOfSteps*5+1,0);
+	std::vector<float> ps(this->numberOfSteps*5+1,1);
 	std::vector<float> masks(this->numberOfSteps*5+1,0);
 	std::vector<float> masksW(this->numberOfSteps*5+1,0);
 	std::vector<float> segs(this->numberOfSteps*5+1,0);
@@ -178,7 +183,7 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 		double mag, val;
 		bool count =false;
 		int counting = 0;
-		for(int t=-this->numberOfSteps; (t<=this->numberOfSteps || (ps[t+numberOfSteps-1] >1.0e-20  && t +numberOfSteps< ps.size()));t++)
+		for(int t=-this->numberOfSteps; (t<=this->numberOfSteps || (ps[t+numberOfSteps-1] > (1.0/pow(10,this->minPGrey))  && t +numberOfSteps< ps.size()));t++)
 		{
 			//ps[t+numberOfSteps] =1;
 
@@ -205,7 +210,7 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 
 				magnitudes[t+numberOfSteps]+= fabs(mag) /images.size();
 				intensities[t+numberOfSteps]+=val/images.size();
-				ps[t+numberOfSteps]+=(1-whiteIntensity) * p/images.size();
+				ps[t+numberOfSteps]+=(1-whiteIntensity) * p; ///images.size();
 			}
 			if(vertexDebug==j)
 			{
@@ -238,7 +243,7 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 			float changeHemisAseg=0;
 			bool good=true; 
 			int zeroLabel=0;
-			for(int i=1;(i< numberOfSteps+1 || ( ps[i-1] > 1e-15 && i < ps.size() -1));i++)
+			for(int i=1;(i< numberOfSteps+1 || ( ps[i-1] > (1.0/pow(10,this->minPGrey)) && i < ps.size() -1));i++)
 			{
 				if (vertexDebug ==j)
 					std::cout <<" opt mag " << opt_mag << " mag " << magnitudes[i] << " " << ps[i] << " " << ps[i-1] << " " << ps[i-2] << std::endl;
@@ -275,18 +280,18 @@ void MRIS_MultimodalRefinement::getTarget(MRIS* surf )
 				if (vertexDebug ==j)
 					std::cout << lastCortexLabel << " " << segs[i+1] << " " << touchedStructure <<   std::endl;
 			//label = (label>.5)?1:0;
-				if((fabs(magnitudes[i])>opt_mag && ps[i] > 1e-15)  || touchedStructure>0 || (changeHemis > 0  && fabs(magnitudes[i]) *1.5 > opt_mag) ||( changeHemisAseg >0&& zeroLabel <3 &&opt_mag < 5))
+				if((fabs(magnitudes[i])>opt_mag && ps[i] > (1.0/pow(10,this->minPGrey)))  || touchedStructure>0 || (changeHemis > 0  && fabs(magnitudes[i]) *1.5 > opt_mag) ||( changeHemisAseg >0&& zeroLabel <3 &&opt_mag < 5))
 				{
 
 					int w=i;
 					good = true;
-					for( ; w>0 && good && (ps[w-1] <1e-05  || fabs(i-w) <1);w--)
+					for( ; w>0 && good && (ps[w-1] < (1.0/pow(10, this->minPGrey /3.0)) || fabs(i-w) <1);w--)
 					{
 						good = ps[w-1] >= ps [w] && masks[w] <.1 ;
 						if(vertexDebug==j)
 							std::cout << ps[w-1 ] << " " << ps[w] <<  " " << w << " " << masks[w]<< std::endl;
 					}
-					if((((ps[w-1] > 1e-05 )|| good ) || (touchedStructure >0  && ps[i-1] > 1e-1 )) && label != 41 && label !=2 )  
+					if((((ps[w-1] > (1.0/pow(10, this->minPGrey/3.0)) )|| good ) || (touchedStructure >0  && ps[i-1] > 0.0001 )) && label != 41 && label !=2 )  
 					{
 						opt_mag = fabs( magnitudes[i]);
 						opt_val= intensities[i];
@@ -647,7 +652,7 @@ std::vector<std::pair<float, float>> MRIS_MultimodalRefinement::GetMeanAndVarian
 	std::cout << " hi " << std::endl;
 	return mean_variance;
 }
-void  MRIS_MultimodalRefinement::SegmentWM(MRI* t1, MRI* t2, MRI* output)
+void  MRIS_MultimodalRefinement::SegmentWM(MRI* t1, MRI* t2, MRI* output, int contrast_type)
 {
 	for (int x = 0 ; x < t1->width ; x++)
 	{
@@ -658,8 +663,17 @@ void  MRIS_MultimodalRefinement::SegmentWM(MRI* t1, MRI* t2, MRI* output)
 				//int label =(imageAseg)? MRIgetVoxVal(imageAseg, x, y, z, 0) :0;
 				float T1val = MRIgetVoxVal(t1, x, y, z, 0) ;
 				float T2val = MRIgetVoxVal(t2, x, y, z, 0) ;
+
 				int val = 0;
-				if(1.3*T1val > T2val &&  (T1val +  T2val) >= 180 && T1val > 80 && T2val > 80)
+				if(contrast_type <0)
+				{
+					double label;
+					MRIsampleVolumeFrameType(t2, x, y, z, 0, SAMPLE_NEAREST, &label);
+					if(label==2 || label==41)
+					{
+						val=1;
+					}
+				} else 	if(1.3*T1val > T2val &&  (T1val +  T2val) >= 180 && T1val > 80 && T2val > 80)
 				{
 					val =1;
 				}
@@ -671,8 +685,19 @@ void  MRIS_MultimodalRefinement::SegmentWM(MRI* t1, MRI* t2, MRI* output)
 
 
 }
-void  MRIS_MultimodalRefinement::SegmentVessel(MRI* t1, MRI* t2, MRI* output)
+void  MRIS_MultimodalRefinement::SegmentVessel(MRI* t1, MRI* t2, MRI* output, int contrast_type)
 {
+
+	int sum=180;
+	int max=80;
+	int min=110;
+	float mult=1.4;
+	std::cout << contrast_type << " " << CONTRAST_FLAIR<<" flair" <<std::endl;
+	if (contrast_type == CONTRAST_FLAIR)
+	{
+		mult=1;
+		max=50;
+	}
 	for (int x = 0 ; x < t1->width ; x++)
 	{
 		for (int y = 0 ; y < t1->height ; y++)
@@ -683,7 +708,16 @@ void  MRIS_MultimodalRefinement::SegmentVessel(MRI* t1, MRI* t2, MRI* output)
 				float T1val = MRIgetVoxVal(t1, x, y, z, 0) ;
 				float T2val = MRIgetVoxVal(t2, x, y, z, 0) ;
 				int val = 0;
-				if (  (1.4*T1val > T2val && (T1val +T2val <180)) || (T1val <80 && T2val <80) || (T1val > 110 && T2val>110))
+				if( contrast_type <0)
+				{  
+					double label;
+					MRIsampleVolumeFrameType(t2, x, y, z, 0, SAMPLE_NEAREST, &label);
+					/*if(label==0)
+					{
+						val=1;
+					}*/
+				
+				} else if (  (mult*T1val > T2val && (T1val +T2val <sum)) || (T1val <max && T2val <max) || (T1val > min && T2val>min))
 				{
 					val=1;
 				}
diff --git a/resurf/mri_vessel_segment.cxx b/resurf/mri_vessel_segment.cxx
index 07d62014c87..69b6615b2de 100644
--- a/resurf/mri_vessel_segment.cxx
+++ b/resurf/mri_vessel_segment.cxx
@@ -57,7 +57,7 @@ int main(int narg, char * arg[])
 		MRI *vesselMR =  MRIcopy(imageT1, NULL) ;
 		//MRI *whiteMR=  MRIcopy(imageT1, NULL) ;
 		MRIS_MultimodalRefinement refinement;
-		refinement.SegmentVessel(imageT1, imageT2, vesselMR);
+		refinement.SegmentVessel(imageT1, imageT2, vesselMR, 0);
 
 		MRIwrite(vesselMR,outputName) ;
 		MRIfree(&imageT1);	
diff --git a/resurf/mris_extract_values.cxx b/resurf/mris_extract_values.cxx
index ae7056eec6c..3d7c9019e9f 100644
--- a/resurf/mris_extract_values.cxx
+++ b/resurf/mris_extract_values.cxx
@@ -41,6 +41,8 @@
 #include <vtkPolyLine.h>
 #include <vtkPolyData.h>
 
+#include "annotation.h"
+#include "surfcluster.h"
 int main(int narg, char*  arg[])
 {
 //	try{
@@ -58,78 +60,123 @@ int main(int narg, char*  arg[])
 		if(cl.size()==1 || cl.search(2,"--help","-h"))
 		{
 			std::cout<<"Usage: " << std::endl;
-			std::cout<< arg[0] << " -i surface -v overlay -o csvfile -m numberOfImages image1 image2 imaga3"  << std::endl;   
+			std::cout<< arg[0] << " -i surface -v overlay -a annotation -o csvfile -m numberOfImages image1 image2 imaga3"  << std::endl;   
 			return -1;
 		}
 		const char *inSurf= cl.follow ("", "-i");
 		const char *outCsv = cl.follow ("", "-o");
 		const char *overlayFilename= cl.follow ("", "-v");
+		const char *annotFilename= cl.follow ("", "-a");
+		
 		MRI_SURFACE *surf;
 		surf = MRISread(inSurf);
 		MRIScomputeMetricProperties(surf);
 		MRISstoreMetricProperties(surf);
 
 		int imageNumber = cl.follow(0,"-m");
-		
-		std::vector<MRI*> images; 
-
-		std::vector<std::string> fileNames;
-		for(;imageNumber>0;imageNumber--)
+		if( cl.search("-a") )
 		{
-			fileNames.push_back(cl.next(""));
-			MRI *imageFS =  MRIread(fileNames[fileNames.size()-1].c_str()) ;
-			images.push_back(imageFS);
+			MRISreadCurvature(surf,overlayFilename) ;
+			MRISreadAnnotation(surf,annotFilename) ;
+			sclustCountClusters(surf);
+			//read_named_annotation_table(annotFilename);
+			std::vector<int> vector = readAnnotationIntoVector(annotFilename);
+			std::map<int,int> myannot; 
+			for(int i=0;i<vector.size();i++)
+			{
+				myannot[vector[i]]=i;
+			}
+		
+			std::map<int,std::vector<float>> clusterVals;
+			for (unsigned j=0;j<surf->nvertices;j++)
+			{
+				//std::cout <<  surf->vertices[j].undefval ;
+				if( clusterVals.count(surf->vertices[j].annotation) ==0)
+					clusterVals[surf->vertices[j].annotation]= std::vector<float>();
+				clusterVals[surf->vertices[j].annotation].push_back(surf->vertices[j].curv);
+
+			}
+
+			MRISfree(&surf);	
+			std::ofstream csv;
+			csv.open(outCsv);
+			csv << " cluster, mean,std  "<< std::endl;
+			int i =0;
+			for( auto& x:clusterVals)
+			{	
+				float mean = std::accumulate(x.second.begin(), x.second.end(), 0.0)/x.second.size();
+				double accum = 0.0;
+				std::for_each (std::begin(x.second), std::end(x.second), [&](const double d) {
+					accum += (d - mean) * (d - mean);
+				});
+
+				double stdev = sqrt(accum / (x.second.size()-1));
+				csv << i << "," <<   x.first <<  ","<< mean << ","<<stdev << std::endl;
+				i++;
+			}
+			csv.close();
 		}
+		else
+		{
+			std::vector<MRI*> images; 
 
-		double x,y,z;
-		MRISreadCurvature(surf,overlayFilename) ;
-		
-		std::ofstream csv;
-		csv.open(outCsv);
-		csv << " Intensity, PrevIntensity, NextIntensity, Overlay "<< std::endl;
+			std::vector<std::string> fileNames;
+			for(;imageNumber>0;imageNumber--)
+			{
+				fileNames.push_back(cl.next(""));
+				MRI *imageFS =  MRIread(fileNames[fileNames.size()-1].c_str()) ;
+				images.push_back(imageFS);
+			}
 
-		for (unsigned j=0;j<surf->nvertices;j++)
-		{
-			
-			double xv, yv,zv, x,y,z, nx,ny,nz,val;
-		        float point[3]={0,0,0};
-			float normal[3]={0,0,0};
-	
-			nx= surf->vertices[j].nx;
-			ny=surf->vertices[j].ny;
-			nz=surf->vertices[j].nz;
-			float dist = sqrt(nx*nx + ny*ny + nz*nz);
-			if( dist>0)
+			double x,y,z;
+			MRISreadCurvature(surf,overlayFilename) ;
+
+			std::ofstream csv;
+			csv.open(outCsv);
+			csv << " Intensity, PrevIntensity, NextIntensity, Overlay "<< std::endl;
+
+			for (unsigned j=0;j<surf->nvertices;j++)
 			{
-				nx /= dist;
-				ny /= dist;
-				nz /= dist;
-			}	
-
-			MRISsurfaceRASToVoxel(surf, images[0], surf->vertices[j].x,surf->vertices[j].y,surf->vertices[j].z, &xv,&yv,&zv);
-			MRIsampleVolume(images[0], xv, yv, zv, &val);
-			csv << val <<", ";
-			x=surf->vertices[j].x -nx;
-			y=surf->vertices[j].y -ny;
-			z=surf->vertices[j].z -nz;
-
-			MRISsurfaceRASToVoxel(surf, images[0], x,y,z, &xv,&yv,&zv);
-
-			MRIsampleVolume(images[0], xv, yv, zv, &val);
-			csv << val <<", ";
-			x=surf->vertices[j].x +nx;
-			y=surf->vertices[j].y +ny;
-			z=surf->vertices[j].z +nz;
-
-			MRISsurfaceRASToVoxel(surf, images[0], x,y,z, &xv,&yv,&zv);
-
-			MRIsampleVolume(images[0], xv, yv, zv, &val);
-			csv << val <<", ";
-			csv<< surf->vertices[j].curv << std::endl;
-		}
-		MRISfree(&surf);	
-		csv.close();
 
+				double xv, yv,zv, x,y,z, nx,ny,nz,val;
+				float point[3]={0,0,0};
+				float normal[3]={0,0,0};
+
+				nx= surf->vertices[j].nx;
+				ny=surf->vertices[j].ny;
+				nz=surf->vertices[j].nz;
+				float dist = sqrt(nx*nx + ny*ny + nz*nz);
+				if( dist>0)
+				{
+					nx /= dist;
+					ny /= dist;
+					nz /= dist;
+				}	
+
+				MRISsurfaceRASToVoxel(surf, images[0], surf->vertices[j].x,surf->vertices[j].y,surf->vertices[j].z, &xv,&yv,&zv);
+				MRIsampleVolume(images[0], xv, yv, zv, &val);
+				csv << val <<", ";
+				x=surf->vertices[j].x -nx;
+				y=surf->vertices[j].y -ny;
+				z=surf->vertices[j].z -nz;
+
+				MRISsurfaceRASToVoxel(surf, images[0], x,y,z, &xv,&yv,&zv);
+
+				MRIsampleVolume(images[0], xv, yv, zv, &val);
+				csv << val <<", ";
+				x=surf->vertices[j].x +nx;
+				y=surf->vertices[j].y +ny;
+				z=surf->vertices[j].z +nz;
+
+				MRISsurfaceRASToVoxel(surf, images[0], x,y,z, &xv,&yv,&zv);
+
+				MRIsampleVolume(images[0], xv, yv, zv, &val);
+				csv << val <<", ";
+				csv<< surf->vertices[j].curv << std::endl;
+			}
+			MRISfree(&surf);	
+			csv.close();
+		}
 	/*
 	}catch(...)
 	{
diff --git a/resurf/mris_multimodal_surface_placement.cxx b/resurf/mris_multimodal_surface_placement.cxx
index f9c7737922d..0239806af41 100644
--- a/resurf/mris_multimodal_surface_placement.cxx
+++ b/resurf/mris_multimodal_surface_placement.cxx
@@ -99,6 +99,7 @@ int main(int narg, char*  arg[])
 
 		MRIS_MultimodalRefinement* t2refinement = new MRIS_MultimodalRefinement();	
 		std::vector<MRI*> images; 
+		int modality=0;
 		if( cl.search("-t1"))
 		{
 			MRI *t1 = MRIread(cl.follow("","-t1"));
@@ -110,18 +111,20 @@ int main(int narg, char*  arg[])
 			MRI *t2 = MRIread(cl.follow("","-t2"));
 			images.push_back(t2);
 			t2refinement->addImage(t2);
+			modality=1;
 		}
 		if( cl.search("-flair"))
 		{
 			MRI *flair = MRIread(cl.follow("","-flair"));
 			images.push_back(flair);
 			t2refinement->addImage(flair);
+			modality=2;
 		}
 
 		MRI* whiteMR= MRIcopy(images[0], NULL);
 		MRI* vesselMR= MRIcopy(images[0], NULL);
-		t2refinement->SegmentVessel(images[0], images[1],vesselMR);
-		t2refinement->SegmentWM(images[0], images[1],whiteMR);
+		t2refinement->SegmentVessel(images[0], images[1],vesselMR, modality);
+		t2refinement->SegmentWM(images[0], images[1],whiteMR,modality);
 		t2refinement->SetWhiteMR(whiteMR);
 		t2refinement->SetVesselMR(vesselMR);