Laboratoire RDP » Biophysics Team » FLORIVAR

import numpy as np

import matplotlib.pyplot as plt

import math

from scipy import optimize

from scipy.ndimage import binary_dilation

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, volumic=True):

	'''

	Centers the image img on the point cm

	'''

	s=img.shape

	res=img.voxelsize

	if volumic:

		obj=img

	else :

		obj=binary_dilation(img)-img

	x = np.array(np.where(obj > 0) )

	cm = np.array([np.mean(x[i]) for i in [0,1,2]])

	x = np.array(np.where(img > 0) )

	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, volumic=True):	

	'''

	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.

	'''

	if volumic:

		obj=img

	else :

		obj=binary_dilation(img)-img

	x,y,z = np.where(obj*255 > 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)

	np.savetxt(trsf_file, T, delimiter=";")

	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 circle(R, c):

	phi=np.linspace(np.pi/6., 5*np.pi/6., 51)

	x=c[0]+R*np.cos(phi)

	y=c[1]+R*np.sin(phi)

	return x,y

def analyze_curve1(name, y, pixelsize, outfilename,  title, 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, pcR = 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, pc21 = compute_curvature_radius(y21[:,0],y21[:,1])

	R22, pc22 = compute_curvature_radius(y22[:,0],y22[:,1])

	R41, pc41 = compute_curvature_radius(y41[:,0],y41[:,1])

	R423, pc423 = compute_curvature_radius(y423[:,0],y423[:,1])

	R44, pc44 = compute_curvature_radius(y44[:,0],y44[:,1])

	fig = plt.figure(figsize=(12, 5))

	plt.plot(yy[:,0], yy[:,1], color='lightgreen', linewidth=10, label='length='+"{:.0f}".format(length)+" um")

	#plt.plot(y[:,0], y[:,1], '--')

	#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(pcR[0], pcR[1], 'ro', color='royalblue')

	plt.plot(pc423[0], pc423[1], 'ro', color='tomato')

	xR,yR = circle(R, pcR)

	plt.plot(xR, yR, "--",color='royalblue', linewidth=3, label='R='+"{:.0f}".format(R)+" um")

	xR,yR = circle(R423, pc423)

	plt.plot(xR, yR, "--", color='tomato', linewidth=3, label='R423='+"{:.0f}".format(R423)+" um")

	plt.gca().set_aspect('equal', adjustable='box')

	plt.xlabel('x (um)')

	plt.ylabel('z (um)')

	plt.ylim([0,1000])

	plt.xlim([0,3000])

	plt.legend()

	plt.grid(True)

	plt.title(title)

	#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

def register_objects(img_flo, img_ref, preregisterd=True, ph=6, pl=2, es=5, fs=3, save_files=True, iterations=2):

	# -- "Specify paths to image files" --

	# ------------------------------------------------

	sequence_name = "register_"

	file_times = [1, 2]

	filenames = [sequence_name + str(t).zfill(2)  for t in file_times]

	list_images = [img_flo, img_ref]

	images = {}

	for i in [0,1]:

		images[filenames[i]] = list_images[i]

	rigid_images = {}

	rigid_transformations = {}

	affine_images = {}

	affine_transformations = {}

	registered_images = {}

	non_linear_transformations = {}

	[reference_filename, floating_filename, reference_time, floating_time] = [ filenames[0], filenames[1], file_times[0], file_times[1]]

	#for reference_filename, floating_filename, reference_time, floating_time, in zip(filenames[1:],filenames[:-1],file_times[1:],file_times[:-1]):

	print("================================")

	print("Registering ", floating_time, " on ", reference_time)

	print("================================")

	###################################################

	# Output structure

	###################################################

	namestring = "_ph"+str(ph)+"_pl"+str(pl)+"_es"+str(es)+"_fs"+str(fs)

	output_dirname = dirname + '/' + sequence_name + "_" + str(floating_time).zfill(2) + "_on_" + str(reference_time).zfill(2) + "" + namestring

	if not os.path.exists(output_dirname):

		os.makedirs(output_dirname)

	reference_img = images[reference_filename]

	floating_img = images[floating_filename]

	####################################################################

	# Registration of the floating images on ref

	# -----------------------------------------------------

	if (not preregistered) :

		# A. Find optimised rigid transformations

		rigid_img, rigid_trsf = find_rigid_transfo2(reference_img, floating_img)

		rigid_images[floating_filename] = rigid_img

		rigid_transformations[floating_filename] = rigid_trsf

		if save_files:

			# rigid-registered images

			rigid_filename = output_dirname + '/' + floating_filename + "_rigid_on_" + str(reference_time).zfill(2) + ".inr.gz"

			imsave(rigid_filename, rigid_img)

			# optimised (with block-matching) rigid transformations

			rigid_trsf_filename = output_dirname + '/tr_' + floating_filename + '_rigid.txt'

			tfsave(rigid_trsf_filename, rigid_trsf)

		# B. Find optimised affine transformations (initialised by the rigid tfs computed above)

		affine_img, affine_trsf = find_affine_transfo(reference_img, floating_img, init_trsf=rigid_trsf)

		affine_images[floating_filename] = affine_img

		affine_transformations[floating_filename] = affine_trsf

		if save_files:

			# affine-registered images

			affine_filename = output_dirname + '/' + floating_filename + "_affine_on_" + str(reference_time).zfill(2) + ".inr.gz"

			imsave(affine_filename, affine_img)

			# optimised (with block-matching) affine transformations

			affine_trsf_filename = output_dirname + '/tr_' + floating_filename + '_affine.txt'

			tfsave(affine_filename, affine_trsf)

		floating_img0=affine_img

	else :

		floating_img0=floating_img

	# C. Find nonlinear transformations which register the affine-transformed images

	#registered_img, nl_trsf = find_nl_transfo(reference_img, affine_img, ph=ph, pl=pl, es=es, fs=fs)

	for it in range(iterations):

		registered_img, nl_trsf = find_nl_transfo(reference_img, floating_img0, ph=ph, pl=pl, es=es, fs=fs)

		registered_images[floating_filename] = registered_img

		non_linear_transformations[floating_filename] = nl_trsf

		# paths to nonlinearly registered images

		registered_filename = output_dirname + '/' + floating_filename + "_registered_on_" + str(reference_time).zfill(2) + "_it"+str(it)+".inr.gz"

		imsave(registered_filename, registered_img)

		floating_img0=registered_img

		if save_files:

			# non-linear transformations

			nl_trsf_filename = output_dirname + '/' + floating_filename + "_it"+str(it)+"_vectorfield.inr.gz"

			save_trsf(nl_trsf, nl_trsf_filename, compress=True)

	return registered_img

#!/usr/bin/env python

# before lounching the executable activate the timagetk conda environement

# conda activate florivar

# general imports :

from sys import argv

import pandas as pd

import os

def read_areas(dirname):

    listfiles=os.listdir(dirname)

    listfiles=[i for i in listfiles if '-area.csv' in i]

    sepalname=[]

    area=[]

    for i in listfiles:

        name=i.split('-area.csv')[0]

        sepalname.append(name)

        df=pd.read_csv(dirname+i, delimiter=',')[:-2]

        A=sum(df['Value'])/2.

        area.append(A)

    data = {'Sepal': sepalname, 'Area': area}

    newdf=pd.DataFrame.from_dict(data)    

    newdf=newdf.sort_values("Sepal", ascending=True, inplace=False)

    return newdf

def read_measurements(fname):

    df=pd.read_csv(fname, delimiter=';')

    print("Reading file : ",fname)

    print("Number of entries at reading: ",df.shape)

    df['Sepal']=[name.split('_oriented')[0] for name in df['Sepal']]

    df=df.sort_values("Sepal", ascending=True, inplace=False)

    return df

measures_filename = argv[1]

areadir = argv[2]+'/'

root=measures_filename.split('.csv')[0]

outputfile=root+"_full.csv"

df=read_measurements(measures_filename)

newdf=read_areas(areadir)

#

df = pd.merge(df, newdf, left_on="Sepal", right_on="Sepal")

#

df['CurvedSF']=df['CurvedWidth']/df['CurvedLength']

df['FlatSF']=df['FlatWidth']/df['FlatLength']

#

df['Thickness']=df['volume']/df['Area']

#

df['GaussK']=1./df['LR']/df['TR']

df['CentralGaussK']=1./df['LR423']/df['TR423']

df.to_csv(outputfile, index=False, sep=";")

print("Number of entries : ",df.shape)

labels=list(df.columns[1:])

print('Columns written:',list(df.columns))

print('Sepals:',list(df['Sepal']))

print('Measurements written in file:',outputfile)

#!/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, tfread, SpatialImage

from titk_tools.registration import isotropic_resampling_seg, isotropic_resampling

from titk_tools.registration import apply_transfo_on_seg1, apply_transfo_on_image

from bib_sepals import *

print('============== orient_sepal.py =================')

filename = argv[1]

outputdir=argv[2]+'/'

segmented=True

if len(argv)>3:

	imagetype=argv[3]

	if imagetype=='grey':

		segmented=False

print(argv)

print('segmented=',segmented)

imname = filename.split('/')[-1]

imnameroot = imname.split('.')[0]

outfilename = outputdir+imnameroot+'_oriented.inr.gz'

# read the image

img=imread(filename)

# resample isotropically

print("----------------------")

print("isotropic resampling of : filename=", filename)

print("----------------------")

voxelsize = img.voxelsize

print("voxelsize=",voxelsize)

s = img.shape

imgtype=img.dtype

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(img>0, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])

if segmented:

	sepale=isotropic_resampling_seg(sepale,isovs)

else:

	img=isotropic_resampling(img,isovs)

	sepale=isotropic_resampling_seg(sepale,isovs)

s = sepale.shape

voxelsize = sepale.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_addz=int(s[2]*0.20)

sepale=add_side_slices(sepale, 0, 0, s_addy, s_addy, s_addz, s_addz)

if not segmented :

	img=add_side_slices(img, 0, 0, s_addy, s_addy, s_addz, s_addz)

'''

# 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, volumic=True)

if not segmented:

	img=center_on_cm(img, volumic=True)

# compute principal directions

pv_flo=principal_directions(sepale_centered, volumic=True)

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)

trsf=tfread(trsf_file)

os.system("rm "+trsf_file)

# apply transformation

if segmented :

	oriented=apply_transfo_on_seg1(trsf, sepale_centered)

else :

	oriented=apply_transfo_on_image(trsf, img)

# smooth the transformed image

#sepale = binary_closing(oriented, structure=struct_el, iterations=radius)

# and save the oriented sepal

sepale_filename = outputdir+imnameroot+'_oriented.inr.gz'

if segmented :

	sepale=SpatialImage(oriented*255, dtype='uint8', voxelsize=voxelsize, origin=[0, 0, 0])

else:

	sepale=SpatialImage(oriented, dtype=imgtype, voxelsize=voxelsize, origin=[0, 0, 0])

imsave(sepale_filename, sepale)

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

//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+"-newtif/";

print("Creating output directory ",dirout);

File.makeDirectory(dirout);

for (FileInd=0; FileInd<list.length; FileInd++){

	FileName = list[FileInd];

	if(endsWith (FileName, ".tif")){

		print ("### Processing ",FileName," ###");

		path = dir+FileName;

		ShortName=substring(FileName,0,indexOf(FileName,".tif"));

		//run("Bio-Formats Importer", "open=["+path+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack open_all_series ");

		TiffName=ShortName+".tif";

		open(path);

		saveAs("Tiff", dirout+TiffName);

		run("Close All");

		print("Done with", FileName,"!");

	} else {

	print("### ",FileName," Not a .tif file ###");

	};

};	

print("Done with this folder!!!");

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

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

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

//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 flip-vertically-mgx.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+"-flipped/";

print("Creating output directory ",dirout);

File.makeDirectory(dirout);

for (FileInd=0; FileInd<list.length; FileInd++){

	FileName = list[FileInd];

	if(endsWith (FileName, ".tif")){

		print ("### Processing ",FileName," ###");

		path = dir+FileName;

		pathout = dirout+FileName;

		print("avant open", FileName,"!");

		//run("Bio-Formats Importer", "open=["+path+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack open_all_series ");

		open(path);

		print("apres open", FileName,"!");

		run("Flip Vertically", "stack");

		saveAs("Tiff", pathout);

		print("avant close", FileName,"!");

		run("Close All");

		print("Done with", FileName,"!");

	} else {

	print("### ",FileName," Not a .tif file ###");

	};

};	

print("Done with this folder!!!");

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

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

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

//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, ".lsm")){

		print ("### Processing ",FileName," ###");

		path = dir+FileName;

		ShortName=substring(FileName,0,indexOf(FileName,".lsm"));

		//run("Bio-Formats Importer", "open=["+path+"] color_mode=Default view=Hyperstack stack_order=XYCZT use_virtual_stack open_all_series ");

		TiffName=ShortName+".tif";

		open(path);

		saveAs("Tiff", dirout+TiffName);

		run("Close All");

		print("Done with", FileName,"!");

	} else {

	print("### ",FileName," Not an .lsm file ###");

	};

};	

print("Done with this folder!!!");

#!/usr/bin/env python

'''

Usage:

measure_sepal.py filename outputdir

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

outputdir=outdir+"/sections/"

projectiondir=outdir+"/projections"

print(argv)

imname = filename.split('/')[-1]

imnameroot = imname.split('.')[0]

projected2Dimage=projectiondir+"/"+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

imtype=sepale.dtype

print("shape=",s)

print("Computing the volume and volumic asymmetry...")

vol=len(np.where(sepale)[0])

# left-right symmetry

sep_left=sepale[:int(s[0]/2), :, :]

sep_right=sepale[int(s[0]/2):, :, :]

vol_left=len(np.where(sep_left)[0])

vol_right=len(np.where(sep_right)[0])

Avolx=(vol_left-vol_right+0.0)/(vol+0.0)

# top-base symmetry

sep_left=sepale[:, :int(s[1]/2), :]

sep_right=sepale[:, int(s[1]/2):, :]

vol_left=len(np.where(sep_left)[0])

vol_right=len(np.where(sep_right)[0])

Avoly=(vol_left-vol_right+0.0)/(vol+0.0)

vol=vol*voxelsize[0]*voxelsize[1]*voxelsize[2]

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;volume;AsymL;AsymT\n')

	FILE.close()

FILE = open(outputfile,"a")

FILE.write(imnameroot+';')

print("Longitudinal section:")

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, "Longitudinal", graph_name)

os.system("rm "+outfile)

FILE.write(str(length_x)+';'+str(flatlength_x)+';'+str(height_x)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+';')

print("Transversal section:")

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, "Transversal", graph_name)

os.system("rm "+outfile)

FILE.write(str(length_y)+';'+str(flatlength_y)+';'+str(height_y)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+';')

FILE.write(str(vol)+';'+str(Avoly)+';'+str(Avolx)+'\n')

import os

'''

Input:

 - directory with detected sepal contour stacks .tif

Output:

 - segmented meshes (up and down surfaces treated as cells)

 - csv file / sepale containing : up and down surface areas

'''

# directory with detected sepal contours

filedata1="/home/akiss/RDP/Projects/Sepals/RNA-Seq_Final/06.ContourDetection/JL_22C/mgx-contours"

################### Parameters ###########

# cubesize for mesh creation

cubesize=str(15)

# size of neighbourhood for the segmented curvature map

curv_neighb=str(100)

# segmentation

signal_ratio=str(3) # minimal border / cell signal ratio (if less, cells are fused)

# - cell

cell_blur_radius=str(100)

seed_radius=str(150)

# - border

border_distance=str(100)

border_blur_radius=str(100)

#########################################

filedata=filedata1+'/'

meshresult=filedata1+'-MGX-mesh'+cubesize+'/'

curvaturecsvresult=filedata1+'-MGX-curvature-csv/'

arearesult=filedata1+'-MGX-area/'

segresult=filedata1+'-MGX-seg/'

os.system("mkdir "+meshresult)

os.system("mkdir "+curvaturecsvresult)

os.system("mkdir "+arearesult)

os.system("mkdir "+segresult)

listfiles=os.listdir(filedata)

listfiles=[ i for i in listfiles if '.tif' in i]

Process.Stack__System__Clear_Work_Stack('0')

Process.Stack__System__Clear_Main_Stack('0')

Process.Stack__System__Clear_Main_Stack('1')

Process.Stack__System__Clear_Work_Stack('1')

for i in listfiles:	

	Process.Stack__System__Set_Current_Stack('Main', '0')

	Process.Stack__System__Open(filedata+i, 'Main', '0')

	# --- creating surface mesh

	Process.Mesh__Creation__Marching_Cubes_Surface(cubesize, '5000')

	Process.Mesh__Structure__Smooth_Mesh('1', 'No')

	Process.Mesh__Structure__Smooth_Mesh('1', 'No')

	Process.Mesh__Structure__Smooth_Mesh('1', 'No')

	Process.Mesh__Structure__Smooth_Mesh('1', 'No')

	# --- computing maximal curvature

	csvname=i[:-4]+'-carte-courbure-Maximal'+curv_neighb+'.csv'

	Process.Mesh__Signal__Project_Mesh_Curvature(curvaturecsvresult+csvname, 'Maximal', curv_neighb, 'Yes', '-50.0', '50.0', '85')

	meshname=i[:-4]+'-mesh.mgxm'

	Process.Mesh__System__Save(meshresult+meshname, 'no', '0')

	Process.Stack__System__Clear_Work_Stack('0')

	Process.Stack__System__Clear_Main_Stack('0')

	Process.Stack__System__Clear_Main_Stack('1')

	Process.Stack__System__Clear_Work_Stack('1')

	# --- watershed segmentation of the maximal curvature map

	Process.Mesh__Segmentation__Auto_Segmentation('Yes', 'No', cell_blur_radius, seed_radius, border_blur_radius, '100.0', border_distance, signal_ratio)

	segname=i[:-4]+'-seg.mgxm'

	Process.Mesh__System__Save(segresult+segname, 'no', '0')

	# --- computing the area of the faces

	csvname=i[:-4]+'-area.csv'

	Process.Mesh__Heat_Map__Heat_Map_Classic('Area', 'Geometry', arearesult+csvname, 'Geometry', 'No', '0', '65535', 'Yes', 'No', 'None', 'No', 'Decreasing', 'Ratio', '0.001', '1')

	# --- computing teh perimeter

	#csvname=i[:-4]+'-perimeter.csv'

	#Process.Mesh__Heat_Map__Heat_Map('/Geometry/Perimeter', 'No', 'No', '', 'No', '', '', 'Yes', 'Yes')

	# --- computing gaussian curvature

	#Process.Mesh__Signal__Project_Mesh_Curvature('', 'Gaussian', '100', 'Yes', '-50.0', '50.0', '85')

	"""print("Please put two seeds on the two sides of the sepal :-) ... ")

	a=input("")

	Process.Mesh__Segmentation__Watershed_Segmentation('50000')

	csvname=i[:-4]+'-area.csv'

	Process.Mesh__Heat_Map__Heat_Map_Classic('Area', 'Geometry', csvresult+csvname, 'Geometry', 'No', '0', '65535', 'Yes', 'No', 'None', 'No', 'Decreasing', 'Ratio', '0.001', '1')

	"""

import os

'''

Input:

 - directory with grey images of sepals

Output:

 - detected sepal contour stacks

'''

# directory of grey images

filedata="/home/biophysics/Desktop/Fanfan/Col-0_A-tif"

################### Parameters ###########

# smoothing

gaussianblurvalue=str(1)

# threshold value for the Edgedetect procedure 

edgedetectvalue=str(8000)

# xy dilation value

dilatevalue=str(0)

#########################################

path=filedata+'/'

outputdir=filedata+'-MGX-EdgeDetect_Blur'+gaussianblurvalue+'_ED'+edgedetectvalue+'_d'+dilatevalue+'/'

os.system("mkdir "+outputdir)

listfiles=os.listdir(path)

listfiles=[ i for i in listfiles if '.tif' in i]

Process.Stack__System__Clear_Work_Stack('0')

Process.Stack__System__Clear_Main_Stack('0')

Process.Stack__System__Clear_Main_Stack('1')

Process.Stack__System__Clear_Work_Stack('1')

for i in listfiles:	

	Process.Stack__System__Set_Current_Stack('Main', '0')

	Process.Stack__System__Open(path+i, 'Main', '0')

	Process.Stack__System__Open(path+i, 'Main', '1')

	Process.Stack__Filters__Gaussian_Blur_Stack(gaussianblurvalue, gaussianblurvalue, gaussianblurvalue)

	Process.Stack__MultiStack__Copy_Work_to_Main_Stack()

	Process.Stack__MultiStack__Swap_or_Copy_Stack_1_and_2('Main', '1 -> 2')

	Process.Stack__Morphology__Edge_Detect(edgedetectvalue, '2.0', '0.3', '30000')

	Process.Stack__MultiStack__Copy_Work_to_Main_Stack()

	Process.Stack__MultiStack__Swap_or_Copy_Stack_1_and_2('Main', '1 <-> 2')

	Process.Stack__System__Set_Current_Stack('Main', '0')

	Process.Stack__Canvas__Reverse_Axes('No', 'No', 'Yes')

	Process.Stack__Morphology__Edge_Detect(edgedetectvalue, '2.0', '0.3', '30000')

	Process.Stack__Canvas__Reverse_Axes('No', 'No', 'Yes')

	Process.Stack__MultiStack__Swap_or_Copy_Stack_1_and_2('Main', '1 <- 2')

	Process.Stack__System__Set_Current_Stack('Both', '0')

	Process.Stack__MultiStack__Combine_Stacks('Product')

	Process.Stack__System__Set_Current_Stack('Work', '0')

	Process.Stack__Morphology__Dilate(dilatevalue, dilatevalue, "0", 'No', 'No')

	Process.Stack__System__Save(outputdir+"MGX_"+i[:-4]+'_e'+edgedetectvalue+'_d'+dilatevalue+'.tif', 'Work', '0', '0')

	Process.Stack__System__Clear_Work_Stack('0')

	Process.Stack__System__Clear_Main_Stack('0')

	Process.Stack__System__Clear_Main_Stack('1')

	Process.Stack__System__Clear_Work_Stack('1')

	Process.Stack__System__Clear_Main_Stack('2')

	Process.Stack__System__Clear_Work_Stack('2')

#!/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 -p "+outdir)

os.system("mkdir -p "+outdir+"/sections")

os.system("mkdir -p "+outdir+"/projections")

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

imagetype=$3

outdir=$indir'_oriented-sepals'

mkdir $indir'/'$outdir

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

	if test $imagetype=="grey"

	then