Laboratoire RDP » Biophysics Team » FLORIVAR

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

os.system("mkdir -p "+symdir)

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+';')

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)

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)+';')

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)

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')