Laboratoire RDP » Biophysics Team » FLORIVAR

#!/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()

#! /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 ..

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

//////////////////////////////////////////////////////////////////////////////////

//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!!!");

#!/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)

#!/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')

//////////////////////////////////////////////////////////////////////////////////

//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!!!");

---------------------------------------------------------------------------

# "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