Index: bin/image2geometry/measure_sepal.py
===================================================================
--- bin/image2geometry/measure_sepal.py	(revision 0)
+++ bin/image2geometry/measure_sepal.py	(revision 1)
@@ -0,0 +1,133 @@
+#!/usr/bin/env python
+
+'''
+Usage:
+measure_sepal.py filename output.csv
+
+- filename = path and name of the stack file containing the object (background=0)
+the object is supposed already oriented as follows:
+- output = output directory
+---------------------------------------------------------------------------------
+- CM in the middle of the image
+- first principal axis along y
+- second principal axis along x
+
+Output:
+
+measured data is written into the file output.csv
+flat projections are written in the output directory
+
+'''
+
+from sys import argv
+import os
+import numpy as np
+from scipy.ndimage.morphology import binary_closing, binary_opening
+
+from titk_tools.io import imread, imsave
+from bib_sepals import *
+
+
+print('============== measure_sepal.py =================')
+
+filename = argv[1]
+outdir = argv[2]
+outputfile=outdir+"/"+outdir+".csv"
+
+
+print(argv)
+
+imname = filename.split('/')[-1]
+imnameroot = imname.split('.')[0]
+
+outputdir = imnameroot+'_sections/'
+projected2Dimage=outdir+"/"+imnameroot+"2D.png"
+projected2Dimage_hight=outdir+"/hight_"+imnameroot+"2D.png"
+
+os.system("mkdir -p "+outputdir)
+
+outname = outputfile.split('/')[-1]
+outroot = outname.split('.')[0]
+
+print("----------------------")
+print("reading the file ", filename)
+print("----------------------")
+sepale=imread(filename)
+voxelsize = sepale.voxelsize
+print("voxelsize=",voxelsize)
+s = sepale.shape
+print("shape=",s)
+
+print("Projecting on 2D...")
+
+def project2D(im3):
+	s2=(s[1],s[0])
+	im2=np.zeros(s2, dtype='uint8')
+	for i in range(s[0]):
+		im2[:,i]=np.array([max(im3[i,j,:]) for j in range(s[1])])
+	return im2
+
+def project2D_hightmap(im3):
+	s2=(s[1],s[0])
+	im2=np.zeros(s2, dtype='uint8')
+	for i in range(s[0]):
+		im2[:,i]=np.array([list(im3[i,j,:]).index(max(im3[i,j,:])) for j in range(s[1])])
+	return im2
+
+
+sepale2D=project2D(sepale)
+imsave2D(projected2Dimage, sepale2D)
+
+'''
+sepale2D=project2D_hightmap(sepale)
+imsave2D(projected2Dimage_hight, sepale2D)
+'''
+
+print('Extracting principal section curves...')
+
+# longitudinal section
+x_section = sepale[int(s[0]/2), :, :]
+xup_filename = outputdir+imnameroot+'_section_xup.txt'
+xdown_filename = outputdir+imnameroot+'_section_xdown.txt'
+x_filename = outputdir+imnameroot+'_section_x.txt'
+x_pixelsize = (voxelsize[1],voxelsize[2])
+xup, xdown = extract_curves(x_section)
+
+write_curve(xup, x_pixelsize, xup_filename)
+
+
+# transversal section
+y_section = sepale[:,int(s[1]/2), :]
+yup_filename = outputdir+imnameroot+'_section_yup.txt'
+ydown_filename = outputdir+imnameroot+'_section_ydown.txt'
+y_filename = outputdir+imnameroot+'_section_y.txt'
+y_pixelsize = (voxelsize[0],voxelsize[2])
+yup, ydown = extract_curves(y_section)
+
+write_curve(yup, y_pixelsize, yup_filename)
+
+
+# 
+# ----------------------------------------
+print('Analyzing principal section curves...')
+
+if not(os.path.isfile(outputfile)):
+	FILE = open(outputfile,"w")
+	FILE.write('Sepal;CurvedLength;FlatLength;LHeight;LR;LR21;LR22;LR41;LR423;LR44;CurvedWidth;FlatWidth;THeight;TR;TR21;TR22;TR41;TR423;TR44\n')
+	FILE.close()
+
+FILE = open(outputfile,"a")
+FILE.write(imnameroot+';')
+
+outfile=outroot+'_up_longitudinal.csv'
+graph_name= outputdir+imnameroot+'_up_longitudinal.png'
+length_x, flatlength_x, height_x, R, R21, R22, R41, R423, R44=analyze_curve1(imnameroot, xup, x_pixelsize, outfile, graph_name)
+
+FILE.write(str(length_x)+';'+str(flatlength_x)+';'+str(height_x)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+';')
+
+outfile=outroot+'_up_transversal.csv'
+graph_name= outputdir+imnameroot+'_up_transversal.png'
+length_y, flatlength_y, height_y, R, R21, R22, R41, R423, R44=analyze_curve1(imnameroot, yup, y_pixelsize, outfile, graph_name)
+
+FILE.write(str(length_y)+';'+str(flatlength_y)+';'+str(height_y)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+'\n')
+

Property changes on: bin/image2geometry/measure_sepal.py
___________________________________________________________________
Added: svn:executable
   + 

Index: bin/image2geometry/lif2tif.ijm
===================================================================
--- bin/image2geometry/lif2tif.ijm	(revision 0)
+++ bin/image2geometry/lif2tif.ijm	(revision 1)
@@ -0,0 +1,66 @@
+//////////////////////////////////////////////////////////////////////////////////
+
+//The script will process all your .lif files in the chosen folder.
+//The stacks will be saved in subfolders named as the treated lif files.
+
+print ("====================================");
+print ("======== Macro lif2tif.ijm =========");
+print ("====================================");
+
+//Choose the directory containing your .lif file//
+dir = getDirectory("Choose a directory")
+setBatchMode(true);
+list = getFileList(dir);
+
+//Find .lif experiment//
+for (FileInd=0; FileInd<list.length; FileInd++){
+	FileName = list[FileInd];
+	if(endsWith (FileName, "lif")){
+		print ("### Processing ",FileName," ###");
+		path = dir+FileName;
+		ShortName=substring(FileName,0,indexOf(FileName,".lif"));
+		print("Creating output directory ",ShortName);
+		File.makeDirectory(dir+ShortName);
+		//Open .lif//
+		run("Bio-Formats Importer", "open=["+path+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack open_all_series ");
+		//Saving as .tif the content of the .lif files, in their destination folders//
+		list = getList("image.titles");
+		//get the name of each image opened//  
+		ids=newArray(nImages);
+		for (j=0; j<list.length; j++){
+			//print(list[j]);
+			selectImage(j+1);
+			ids[j]=getTitle;
+			//Get experiment name (title of the .lif file)//
+			experiment_name=substring(ids[j],0,indexOf(ids[j],".lif"));
+			//Remove the name of the .lif file that is included in each image//
+			short_name1 = substring(ids[j], indexOf(ids[j], "- "), lastIndexOf(ids[j], ""));
+			short_name = replace(short_name1, "- ", "");
+			//If not a stack, save directly as .tif with a short name//
+			if (nSlices==1){				
+				print("Not a stack");
+				run("8-bit");
+				saveAs("Tiff", dir+File.separator+experiment_name+File.separator+short_name);
+				print(short_name, ">>> Image saved as .tif");			
+			//If a stack, reverse and save as .tif with a short name//
+			} else {
+				run("Reverse");
+				run("8-bit");
+				saveAs("Tiff", dir+File.separator+experiment_name+File.separator+short_name);
+				print("----- ",short_name, " >>> Reversed and saved as .tif");
+				};
+			getPixelSize(unit, pw, ph, pd);
+			print("Voxelsize="+pw+"x"+ph+"x"+pd+" "+unit+"^3");
+			print("nSlices="+nSlices);
+			};
+		run("Close All");
+		print("Done with", FileName,"!");
+	} else {
+	print("### ",FileName," Not a .lif file###");
+	};
+};	
+print("Done with this folder!!!");
+
+
+
+
Index: bin/image2geometry/measure_sepals.py
===================================================================
--- bin/image2geometry/measure_sepals.py	(revision 0)
+++ bin/image2geometry/measure_sepals.py	(revision 1)
@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+# Usage:
+# ------
+# measure_sepals.py directory filetype
+# filetype = '.tif' or '.inr.gz'
+
+from sys import path, argv
+import os
+import numpy as np
+from multiprocessing import Pool
+nproc = 11 # number of processors
+
+indir=argv[1]
+filetype=argv[2]
+
+outdir=indir+'measures'
+
+os.system("mkdir "+outdir)
+
+def one_simulation(param):
+	if filetype in param:
+		filename=param
+		os.system("measure_sepal.py "+filename+" "+outdir)
+	return
+
+list_params=[indir+'/'+filename for filename in os.listdir(indir)]
+
+pool = Pool(processes=nproc)
+pool.map(one_simulation, list_params)
+pool.close()
+pool.join()
+

Property changes on: bin/image2geometry/measure_sepals.py
___________________________________________________________________
Added: svn:executable
   + 

Index: bin/image2geometry/orient_sepals.sh
===================================================================
--- bin/image2geometry/orient_sepals.sh	(revision 0)
+++ bin/image2geometry/orient_sepals.sh	(revision 1)
@@ -0,0 +1,37 @@
+#! /bin/bash
+
+# Usage:
+# ------
+# ./orient_sepals.sh directory filetype
+# filetype = '.tif' or '.inr.gz'
+
+
+indir=$1
+filetype=$2
+
+outdir=$indir'_oriented'
+
+
+cd $indir
+
+logfile='orient_sepals-'$indir'.log'
+date > $logfile
+echo '--------------------------' >> $logfile
+
+
+for file in *$filetype
+do
+	echo '--------------------------' >> $logfile
+	echo $file  >> $logfile
+	echo $file $outdir
+	echo '--------------------------' >> $logfile
+	orient_sepal.py $file $outdir
+	echo '--------------------------' >> $logfile
+done
+
+
+oriented='oriented-sepals'
+mkdir $oriented
+scp */*_oriented.inr.gz $oriented
+
+cd ..

Property changes on: bin/image2geometry/orient_sepals.sh
___________________________________________________________________
Added: svn:executable
   + 

Index: bin/image2geometry/bib_sepals.py
===================================================================
--- bin/image2geometry/bib_sepals.py	(revision 0)
+++ bin/image2geometry/bib_sepals.py	(revision 1)
@@ -0,0 +1,447 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import math
+from scipy import optimize
+
+from titk_tools.io import imread, imsave, SpatialImage
+
+def center_on_point(img, cm):
+	'''
+	Centers the image img on the point cm
+	'''
+	x = np.array(np.where(img > 0) )
+	xp=np.zeros_like(x)
+	img1=np.zeros_like(img)
+	s=img.shape
+	center_img=np.array(s)/2
+	for i in [0, 1, 2]:
+		xp[i]=x[i]+center_img[i]-cm[i]
+	for i in range(0, len(x[0])):
+		if ((xp[0][i]>=0) and (xp[0][i]<s[0]) and (xp[1][i]>=0) and (xp[1][i]<s[1]) and (xp[2][i]>=0) and (xp[2][i]<s[2])):
+			img1[xp[0][i],xp[1][i],xp[2][i]]=img[x[0][i],x[1][i],x[2][i]]
+	img1=SpatialImage(img1)
+	img1.resolution=img.resolution
+	return img1
+
+
+def center_on_cm(img):
+	'''
+	Centers the image img on the point cm
+	'''
+	s=img.shape
+	res=img.voxelsize
+	x = np.array(np.where(img > 0) )
+	cm = np.array([np.mean(x[i]) for i in [0,1,2]])
+	xp=np.zeros_like(x)
+	img1=np.zeros_like(img)
+	center_img=np.array(s)/2
+	for i in [0, 1, 2]:
+		xp[i]=x[i]+center_img[i]-cm[i]
+	for i in range(0, len(x[0])):
+		if ((xp[0][i]>=0) and (xp[0][i]<s[0]) and (xp[1][i]>=0) and (xp[1][i]<s[1]) and (xp[2][i]>=0) and (xp[2][i]<s[2])):
+			img1[xp[0][i],xp[1][i],xp[2][i]]=img[x[0][i],x[1][i],x[2][i]]
+	img1=SpatialImage(img1)
+	img1.voxelsize=res
+	return img1
+
+
+def principal_directions(img, threshold=10):	
+	'''
+	Returns the first two principal eigen vectors of an input image
+	Fixed threshold of 10 for a 8 bit image
+	It suppposes an isotropic sampling of the image.
+	'''
+	x,y,z = np.where(img > threshold) ## appropriate thresholding 
+	x = x - np.mean(x)
+	y = y - np.mean(y)
+	z = z - np.mean(z)
+	coords = np.vstack([x,y,z])
+	cov = np.cov(coords)
+	evals,evecs = np.linalg.eig(cov)
+	sort_indices = np.argsort(evals)[::-1]
+	return list(evecs[:,sort_indices][0:2])
+
+def direct_frame(pv):
+	'''
+	Constructs a direct frame with first two unitvectors v0 and v1
+	'''
+	return np.array([pv[0], pv[1], np.cross(pv[0], pv[1])])
+
+
+def write_transformation(T, trsf_file):
+	'''
+	writes the affine transformation T (4x4 matrix) in
+	a textfile read by blockmatching.
+	'''
+	f = open(trsf_file,'w')
+	f.write("( "+"\n")
+	f.write("08 "+"\n")
+	for i in T:
+		for j in i:
+			f.write("  "+str(j)+" ")
+		f.write("\n")
+	f.write(") ")
+	return
+
+
+
+def write_rigid(trsf_file, R, a):
+	'''
+	write a rigid transformation with rotation matrix R and translation a
+	'''
+	T = np.identity(4)
+	T[0:3,0:3] = R
+	T = np.transpose(T)
+	T[0:3,3] = a
+	write_transformation(T, trsf_file)
+	return
+
+
+
+def write_superposing_transfo(pv_flo, pv_ref, img, trsf_file):
+	'''
+	- pv_flo and pv-ref are lists of the first two principal vectors of the floating and the reference images respectively
+	- img is the floating image
+	'''
+	# compute the rotation matrix R
+	pframe_ref = direct_frame(pv_ref)
+	pframe_flo = direct_frame(pv_flo)
+	F = np.transpose(pframe_ref)
+	G = np.transpose(pframe_flo)
+	R = np.matmul(G, np.linalg.inv(F))
+	# blockmatching rotates arround the origin,
+	# so a rotation arround the middle of the image contains an additional translation t
+	s = img.shape
+	res = img.voxelsize
+	com = np.array(res)*np.array(s)/2
+	i = np.identity(3)
+	t = np.dot((i-R),np.array(com))
+	R = np.transpose(R)
+	write_rigid(trsf_file, R, t)
+	return
+
+
+def add_side_slices(Img, x0, x1, y0, y1, z0, z1):
+	"""
+	adds 0 value slices on different sides of the image
+	"""
+	s=Img.shape
+	ss=(s[0]+x0+x1,s[1]+y0+y1,s[2]+z0+z1)
+	Img_new = (SpatialImage(np.zeros(ss))).astype('uint8')
+	Img_new[x0:ss[0]-x1,y0:ss[1]-y1,z0:ss[2]-z1]=Img
+	Img_new.voxelsize = Img.voxelsize
+	return Img_new
+
+
+def remove_side_slices(Img, x0, x1, y0, y1, z0, z1):
+	"""
+	removes slices on different sides of the image
+	"""
+	ss=Img.shape
+	Img_new=Img[x0:ss[0]-x1,y0:ss[1]-y1,z0:ss[2]-z1]
+	#Img_new.resolution = Img.resolution
+	return Img_new
+
+
+def measure_hight(a):
+	b= np.where(a)
+	height = 0
+	if b[0].shape[0]>0 :
+		height = b[0][-1]
+	return height
+
+
+def measure_height(a):
+	b= np.where(a)
+	height = [0,0]
+	if b[0].shape[0]>0 :
+		height = [b[0].min(),b[0].max()]
+	return height
+
+
+def extract_curve(x_section):
+	sx = x_section.shape
+	xx = np.zeros((sx[0],2))
+	for i in range(0,sx[0]):
+		xx[i,0] = i
+		xx[i,1] = measure_hight(x_section[i,:])
+	return xx
+
+
+def extract_curves(x_section):
+	sx = x_section.shape
+	xup = np.zeros((sx[0],2))
+	xdown = np.zeros((sx[0],2))
+	for i in range(0,sx[0]):
+		xup[i,0] = i
+		xdown[i,0] = i
+		xup[i,1] = measure_height(x_section[i,:])[1]
+		xdown[i,1] = measure_height(x_section[i,:])[0]
+	return xup, xdown
+
+
+def write_curve(y, pixelsize, filename) :
+	FILE = open(filename,"w")
+	FILE.write('# '+str(pixelsize[0])+' '+str(pixelsize[1])+'\n')
+	s = y.shape
+	#FILE.write('# '+str(s[0])+'\n')
+	for i in range(0,s[0]) :
+		FILE.write(str(y[i,0])+' '+str(y[i,1])+'\n')
+	FILE.close()
+	return
+
+
+def struct_element_sphere(radius):
+	s=int(2*radius+1)
+	structure=np.zeros((s,s,s))
+	center = np.array([radius, radius, radius])
+	Nl=range(s)
+	for i in Nl:
+		for j in Nl:
+			for k in Nl:
+				p=np.array([i,j,k])
+				d=p-center
+				dist=np.sqrt(sum(d*d))
+				if dist<=radius:
+					structure[i,j,k]=1
+	return structure
+
+
+def create_artificial_sepal(AA, RR, resolution):
+	# A = parameters of an ellipse, as projected on xy plane (2d array)
+	# R = curvature radius in directions x and y (2d array)
+	# (A and R are given in micrometers)
+	# resolution = voxelsize of the output stack
+	A = np.array([AA[i]/resolution[i] for i in [0,1]])
+	R = np.array([RR[i]/resolution[i] for i in [0,1]])
+	h = np.array([math.sqrt(R[i]**2-A[i]**2) for i in range(len(A))])
+	zmax=np.array([R[i]-math.sqrt(R[i]**2-A[i]**2) for i in [0,1]])
+	#zmax= int(math.sqrt(zmax[0]*zmax[1]))
+	zmax=zmax.mean()
+	#zmax = 5
+	print(zmax)
+	marge = 0.2
+	# creating the stack and the surface
+	s = (int(A[0]*2*(1+marge)), int(A[1]*2*(1+marge)), int(zmax*(1+2*marge)))
+	cm = np.array(s)/2
+	x = np.array(range(0,s[0]))
+	y = np.array(range(0,s[1]))
+	xx = x - cm[0]
+	yy = y - cm[1]
+	#xgrid, ygrid = np.meshgrid(xx, yy, sparse=True)  # make sparse output arrays
+	#mask=(((xgrid**2)/float(A[0])**2+(ygrid**2)/float(A[1])**2)<1).astype('uint8')
+	z = np.zeros((s[0],s[1]))
+	stack= np.zeros(s).astype('uint8')
+	for i in x:
+		for j in range(0,s[1]):
+			if xx[i]**2/float(A[0])**2+yy[j]**2/float(A[1])**2<=1 :
+				zx = (math.sqrt(R[0]**2-xx[i]**2)-h[0])*(1-abs(yy[j])/A[1])
+				zy = (math.sqrt(R[1]**2-yy[j]**2)-h[1])*(1-abs(xx[i])/A[0])
+				z[i,j] = math.sqrt(zx*zy)
+				#z[i,j] = (zx+zy)/2
+				#z[i,j] = 5
+				stack[i,j,int(zmax*marge+z[i,j])]=1
+	stack = SpatialImage(stack)
+	stack.resolution = resolution
+	return z, stack
+
+
+
+def create_artificial_sepal_ellipse(AA, H, resolution):
+	'''
+	# A = parameters of an ellipse, as projected on xy plane (2d array)
+	# R = curvature radius in directions x and y (2d array)
+	# (A and R are given in micrometers)
+	# resolution = voxelsize of the output stack
+	'''
+	A = np.array([AA[i]/resolution[i] for i in [0,1]])
+	h = H/resolution[2]
+	zmax=h
+	marge = 0.2
+	# creating the stack and the surface
+	s = (int(A[0]*2*(1+marge)), int(A[1]*2*(1+marge)), int(zmax*(1+5*marge)))
+	cm = np.array(s)/2
+	x = np.array(range(0,s[0]))
+	y = np.array(range(0,s[1]))
+	xx = x - cm[0]
+	yy = y - cm[1]
+	#xgrid, ygrid = np.meshgrid(xx, yy, sparse=True)  # make sparse output arrays
+	#mask=(((xgrid**2)/float(A[0])**2+(ygrid**2)/float(A[1])**2)<1).astype('uint8')
+	z = np.zeros((s[0],s[1]))
+	#z = -h*((xgrid**2)/float(A[0])**2+(ygrid**2)/float(A[1])**2)
+	stack= np.zeros(s).astype('uint8')
+	for i in x:
+		for j in y:
+			if xx[i]**2/float(A[0])**2+yy[j]**2/float(A[1])**2<=1 :
+				z[i,j] = -h*((xx[i]**2)/float(A[0])**2+(yy[j]**2)/float(A[1])**2-1)
+				stack[i,j,int(zmax*4*marge+z[i,j])]=1
+	stack = SpatialImage(stack)
+	stack.voxelsize = resolution
+	return z, stack
+
+# -------------------------------------------------------------------
+
+def read_curve(filename) :
+	"""
+	Reads the coordinates of the point-list defining a 2D curve.
+	**Returns**
+	y : 2-column array
+		defines a curve as an ordered list of points, 
+		each row of the array represents a point on the curve and contains
+		its 2D cartesian coordinates
+	pixelsize : tuple (res0, res1) in micrometers
+	"""	
+	x0 = []
+	pixelsize = (1.0,1.0)
+	for line in file(filename):
+		if line[0] == "#" :
+			line = line.lstrip('# ')
+			line = line.rstrip('\n')
+			line_list = [float(x) for x in line.split(' ')]
+			pixelsize = (line_list[0], line_list[1])
+		else :
+			line = line.rstrip('\n')
+			line_list = [float(x) for x in line.split(' ')]
+			x0.append(line_list)
+	y = np.array(x0)
+	return y, pixelsize
+
+
+
+def distance(p1, p2):
+	return np.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)
+
+def distance_to_line(p1, p2, p0):
+	'''
+	gives the distance of point p0 to the line defined by points p1 and p2 (in 2d)
+	'''
+	return abs((p2[1]-p1[1])*p0[0]-(p2[0]-p1[0])*p0[1]+p2[0]*p1[1]-p2[1]*p1[0])/distance(p1,p2)
+
+def curvature(p1,p2,p3):
+	t1 = (p3[0]**2-p2[0]**2+p3[1]**2-p2[1]**2)/(2.0*(p3[1]-p2[1]))
+	t2 = (p2[0]**2-p1[0]**2+p2[1]**2-p1[1]**2)/(2.0*(p2[1]-p1[1]))
+	n1 = (p3[0]-p2[0])/(p3[1]-p2[1])
+	n2 = (p2[0]-p1[0])/(p2[1]-p1[1])
+	pcx = (t1-t2+p3[1]-p1[1])/(n1-n2)
+	pc = [pcx, -n1*pcx+t1/2+p3[1]/2+p1[1]/2]
+	R = distance(p1, pc)
+	return 1.0/R, pc
+
+
+def curvature2(p1,p2,p3):
+	# http://www.ambrsoft.com/TrigoCalc/Circle3D.htm
+	A = p1[0]*(p2[1]-p3[1])-p1[1]*(p2[0]-p3[0])+p2[0]*p3[1]-p3[0]*p2[1]
+	B = (p1[0]**2+p1[1]**2)*(p3[1]-p2[1])+(p2[0]**2+p2[1]**2)*(p1[1]-p3[1])+(p3[0]**2+p3[1]**2)*(p2[1]-p1[1])
+	C = (p1[0]**2+p1[1]**2)*(p3[0]-p2[0])+(p2[0]**2+p2[1]**2)*(p1[0]-p3[0])+(p3[0]**2+p3[1]**2)*(p2[0]-p1[0])
+	D = (p1[0]**2+p1[1]**2)*(p3[0]*p2[1]-p2[0]*p3[1])+(p2[0]**2+p2[1]**2)*(p1[0]*p3[1]-p3[0]*p1[1])+(p3[0]**2+p3[1]**2)*(p2[0]*p1[1]-p1[0]*p2[1])
+	x = -B/(2*A)
+	y = C/(2*A)
+	pc = [x, y]
+	R = distance(p1, pc)
+	return 1.0/R, pc
+
+
+def curve_perle(y, first, last, step):
+	yy=[]
+	for i in range(first, last, step):
+		if y[i][1]>0:
+			yy.append(y[i])
+	if i<last:
+		yy.append(y[last])
+	return np.array(yy)
+
+
+
+def compute_curvature_radius(x,y):
+	# coordinates of the barycenter
+	x_m = np.mean(x)
+	y_m = np.mean(y)
+	def calc_R(c):
+		""" calculate the distance of each 2D points from the center c=(xc, yc) """
+		return np.sqrt((x-c[0])**2 + (y-c[1])**2)
+	#
+	def calc_ecart(c):
+		""" calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
+		Ri = calc_R(c)
+		return Ri - Ri.mean()
+	#
+	center_estimate = x_m, y_m
+	center_2, ier = optimize.leastsq(calc_ecart, center_estimate)
+	#
+	xc_2, yc_2 = center_2
+	Ri_2       = calc_R(center_2)
+	R_2        = Ri_2.mean()
+	residu_2   = sum((Ri_2 - R_2)**2)
+	pc=[xc_2, yc_2]
+	return R_2, pc
+
+
+def analyze_curve1(name, y, pixelsize, outfilename, graph_name='graph.png'):
+	# isotropic pixels
+	y[:,0] = y[:,0]*pixelsize[0]
+	y[:,1] = y[:,1]*pixelsize[1]
+	maxv = [y[:,i].max() for i in [0,1]]
+	# find first and last points on the curve
+	step = 5
+	sep_indices = np.where(y[:,1])[0]
+	first = sep_indices[0]#+step/2
+	last = sep_indices[-1]#-step/2
+	# compute flat length
+	flatlength = distance(y[first],y[last])
+	print("flat length =", flatlength," um")
+	# compute curved length and height (maximal distance to the base)
+	yy = curve_perle(y, first, last, step)
+	length = 0
+	height = 0
+	height_index = 0
+	p1 = yy[0]
+	p2 = yy[-1]
+	for i in range(0, len(yy)-1):
+		length = length + distance(yy[i+1],yy[i])
+		d = distance_to_line(p1, p2, yy[i])
+		if d>height :
+			height = d
+			height_index = i
+	print("length = ",length," um")
+	print("height = ",height," um")
+	middle_point = yy[int(len(yy)/2)]
+	# compute curvature radius by fitting the sepal section to a circle
+	yyy=y[first:last]
+	R, pc = compute_curvature_radius(yyy[:,0],yyy[:,1])
+	print("R=",R," um")
+	# cutting curve in pieces in order to estimate curvature radius on portions of the curve
+	slength=len(yyy)
+	y21=yyy[0:int(slength/2)]
+	y22=yyy[int(slength/2):slength]
+	y41=yyy[0:int(slength/4)]
+	y423=yyy[int(slength/4):3*int(slength/4)]
+	y44=yyy[3*int(slength/4):slength]
+	# and computing the curvature radius for each portion
+	R21, pc = compute_curvature_radius(y21[:,0],y21[:,1])
+	R22, pc = compute_curvature_radius(y22[:,0],y22[:,1])
+	R41, pc = compute_curvature_radius(y41[:,0],y41[:,1])
+	R423, pc = compute_curvature_radius(y423[:,0],y423[:,1])
+	R44, pc = compute_curvature_radius(y44[:,0],y44[:,1])
+	fig = plt.figure()
+	plt.plot(yy[:,0], yy[:,1], 'bo')
+	plt.plot(y[:,0], y[:,1], 'r')
+	plt.plot(y[first][0], y[first][1], 'ro')
+	plt.plot(middle_point[0], middle_point[1], 'ro')
+	plt.plot(y[last][0], y[last][1], 'ro')
+	plt.plot(yy[height_index][0], yy[height_index][1], 'go')
+	plt.plot(pc[0], pc[1], 'bx')
+	plt.gca().set_aspect('equal', adjustable='box')
+	#plt.show()
+	fig.savefig(graph_name)
+	#plt.close()
+	# write the data into a file
+	FILE = open(outfilename,"a")
+	FILE.write(name+';'+str(length)+';'+str(R)+';'+str(flatlength)+';'+str(height)+'\n')
+	FILE.close()
+	return length, flatlength, height, R, R21, R22, R41, R423, R44
+
+def imsave2D(filename, matrix):
+	fig = plt.figure()
+	plt.gca().imshow(np.transpose(matrix), interpolation='none')
+	fig.savefig(filename)
+	return
Index: bin/image2geometry/inr2tif.ijm
===================================================================
--- bin/image2geometry/inr2tif.ijm	(revision 0)
+++ bin/image2geometry/inr2tif.ijm	(revision 1)
@@ -0,0 +1,41 @@
+//////////////////////////////////////////////////////////////////////////////////
+
+//The script will process all your .inr.gz files in the chosen folder.
+//The stacks will be saved in subfolders named as the treated lif files.
+
+print ("===========================");
+print ("==== Macro inr2tif.ijm ====");
+
+
+//Choose the directory containing your .inr.gz files//
+dir = getDirectory("Choose a directory")
+setBatchMode(true);
+list = getFileList(dir);
+
+print(dir);
+ShortNameDir=substring(dir,0,lastIndexOf(dir,"/"));
+print (ShortNameDir);
+dirout=ShortNameDir+"-tif/";
+print("Creating output directory ",dirout);
+File.makeDirectory(dirout);
+
+for (FileInd=0; FileInd<list.length; FileInd++){
+	FileName = list[FileInd];
+	if(endsWith (FileName, ".inr.gz")){
+		print ("### Processing ",FileName," ###");
+		path = dir+FileName;
+		ShortName=substring(FileName,0,indexOf(FileName,".inr.gz"));
+		run("Bio-Formats Importer", "open=["+path+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack open_all_series ");
+		TiffName=ShortName+".tif";
+		saveAs("Tiff", dirout+TiffName);
+		run("Close All");
+		print("Done with", FileName,"!");
+	} else {
+	print("### ",FileName," Not an .inr.gz file ###");
+	};
+};	
+print("Done with this folder!!!");
+
+
+
+
Index: bin/image2geometry/orient_sepal.py
===================================================================
--- bin/image2geometry/orient_sepal.py	(revision 0)
+++ bin/image2geometry/orient_sepal.py	(revision 1)
@@ -0,0 +1,150 @@
+#!/usr/bin/env python
+
+'''
+Usage:
+orient_object.py filename output.txt
+
+- filename = path and name of the stack file containing the object (background=0)
+
+Output:
+the image with the object reoriented:
+- CM in the middle of the image
+- first principal axis along y
+- second principal axis along x
+
+'''
+
+from sys import argv
+import os
+import numpy as np
+from scipy.ndimage.morphology import binary_closing, binary_opening
+
+from titk_tools.io import imread, imsave, SpatialImage
+from titk_tools.rw.registration import isotropic_resampling_seg
+
+from bib_sepals import *
+
+
+
+print('============== orient_sepal.py =================')
+
+filename = argv[1]
+
+print(argv)
+
+imname = filename.split('/')[-1]
+imnameroot = imname.split('.')[0]
+
+outputdir = imnameroot+'_oriented/'
+os.system("mkdir -p "+outputdir)
+
+isodir = imnameroot+'_iso/'
+os.system("mkdir -p "+isodir)
+
+#outname = outputfile.split('/')[-1]
+#outroot = outname.split('.')[0]
+
+inrname = outputdir+imnameroot+'.inr.gz'
+outfilename = outputdir+imnameroot+'_oriented.inr.gz'
+isoname = isodir+imnameroot+'_iso.inr.gz'
+sname = isodir+imnameroot+'_iso_closed.inr.gz'
+
+# resample isotropically
+print("----------------------")
+print("isotropic resampling of : filename=", filename)
+print("----------------------")
+sepale=imread(filename)
+voxelsize = sepale.voxelsize
+print("voxelsize=",voxelsize)
+s = sepale.shape
+print("shape=",s)
+
+radius=int(voxelsize[2]/voxelsize[0])
+isovs=(voxelsize[0]*voxelsize[1]*voxelsize[2])**(1.0/3)
+#isovs=voxelsize[0]
+print("new voxelsize=",isovs)
+
+sepale=SpatialImage(sepale>0, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])
+imsave(inrname, sepale)
+isotropic_resampling_seg(inrname, isoname, isovs)
+
+sepale=imread(isoname)
+s = sepale.shape
+voxelsize = sepale.voxelsize
+print('new voxelsize read=', voxelsize)
+print('new shape', s)
+
+# alongate the target image so that a sepal originally along the diagonal fits into the target image
+print("----------------------")
+print("elongating field : filename=", filename)
+print("----------------------")
+s_addy=int((np.sqrt(s[0]*s[1]*2)-np.sqrt(s[0]*s[1]))/2)
+#s_addy=int(np.sqrt(s[0]*s[1])(1.-np.sqrt(2.))/2. * 1.05)
+s_addz=int(s[2]*0.20)
+sepale_long=add_side_slices(sepale, 0, 0, s_addy, s_addy, s_addz, s_addz)
+sepale=SpatialImage(sepale_long, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])
+imsave(isoname, sepale)
+print('--------------------')
+
+# smooth it with opening
+sepale=imread(isoname)
+s = sepale.shape
+voxelsize = sepale.voxelsize
+struct_el=struct_element_sphere(radius)
+sepale = binary_closing(sepale, structure=struct_el, iterations=int(radius))
+struct_el=struct_element_sphere(radius/2)
+sepale = binary_closing(sepale, structure=struct_el, iterations=int(radius/2))
+sepale=SpatialImage(sepale*255, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])
+imsave(sname, sepale)
+
+
+# put the center of mass in the center of the image
+sepale_centered=center_on_cm(sepale)
+voxelsize = sepale_centered.resolution
+sepale_centered_filename = outputdir+imnameroot+'_centered.inr.gz'
+imsave(sepale_centered_filename, sepale_centered)
+
+
+# compute principal directions
+pv_flo=principal_directions(sepale_centered)
+print(pv_flo)
+
+# sepals originally are oriented with top in more-or less "up" direction on the images
+# while the eigenvectors' orientation is arbitrary
+# --> do not change the original rough orientation
+
+# check logitudinal direction:
+if pv_flo[0][1]<0 :
+	pv_flo[0] = -pv_flo[0]
+
+# check the transversal direction
+if pv_flo[1][0]<0 :
+	pv_flo[1] = -pv_flo[1]
+
+print("New pv_flo",pv_flo)
+
+'''
+pv_ref = [np.array([0,-1,0]),
+np.array([-1,0,0])]
+
+'''
+pv_ref = [np.array([0,1,0]),
+np.array([1,0,0])]
+
+
+# compute superposing transformation
+trsf_file=outputdir+imnameroot+'_rigid.txt'
+write_superposing_transfo(pv_flo, pv_ref, sepale_centered, trsf_file)
+
+# apply transformation
+sepale_filename = outputdir+imnameroot+'_oriented.inr.gz'
+apply_transfo_on_seg1(trsf_file, sepale_centered_filename, sepale_filename)
+
+# smooth the transformed image
+sepale=imread(sepale_filename)
+sepale = binary_closing(sepale, structure=struct_el, iterations=radius)
+
+# and save the oriented sepal
+sepale=SpatialImage(sepale*255, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])
+imsave(sepale_filename, sepale)
+

Property changes on: bin/image2geometry/orient_sepal.py
___________________________________________________________________
Added: svn:executable
   + 

Index: readme.txt
===================================================================
--- readme.txt	(revision 0)
+++ readme.txt	(revision 1)
@@ -0,0 +1,39 @@
+---------------------------------------------------------------------------
+# "FLORIVAR" project
+(project ID : "florivar")
+---------------------------------------------------------------------------
+Biophysics Team - RDP ENS Lyon
+
+Developers :
+--------------
+Annamaria KISS, Corentin MOLLIER, Frenk Diego HARTASANCHEZ
+
+Download the repository / Rapatrier le dépot
+--------------------------------------------
+cd "directory that will contain the project's directory"
+svn --username $USER checkout http://forge.cbp.ens-lyon.fr/svn/florivar
+
+Update (download changes from the repository)
+--------------------------------------------
+cd florivar
+svn up
+
+Commit (upload your changes into the repository)
+-----------------------------------------------
+cd florivar
+svn commit
+
+Changing the file tree
+----------------------
+if you want to add, delete or move individual files in the file-tree of the repository, you have to do that by svn
+... cd in the right subdirectory
+svn add toto
+svn mv toto titi
+svn rm titi
+... and then make a commit (!)
+
+Useful links
+-------------
+svn useful commands : http://wiki.greenstone.org/wiki/index.php/Useful_SVN_Commands
+formats in a wiki : https://forge.cbp.ens-lyon.fr/redmine/help/wiki_syntax.html
+
