Laboratoire RDP » Biophysics Team » FLORIVAR

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