Statistiques
| Révision :

root / bin / image2geometry / measure_sepal.py @ 4

Historique | Voir | Annoter | Télécharger (3,81 ko)

1 1 akiss
#!/usr/bin/env python
2 1 akiss
3 1 akiss
'''
4 1 akiss
Usage:
5 4 akiss
measure_sepal.py filename outputdir
6 1 akiss

7 1 akiss
- filename = path and name of the stack file containing the object (background=0)
8 1 akiss
the object is supposed already oriented as follows:
9 1 akiss
- output = output directory
10 1 akiss
---------------------------------------------------------------------------------
11 1 akiss
- CM in the middle of the image
12 1 akiss
- first principal axis along y
13 1 akiss
- second principal axis along x
14 1 akiss

15 1 akiss
Output:
16 1 akiss

17 1 akiss
measured data is written into the file output.csv
18 1 akiss
flat projections are written in the output directory
19 1 akiss

20 1 akiss
'''
21 1 akiss
22 1 akiss
from sys import argv
23 1 akiss
import os
24 1 akiss
import numpy as np
25 1 akiss
from scipy.ndimage.morphology import binary_closing, binary_opening
26 1 akiss
27 1 akiss
from titk_tools.io import imread, imsave
28 1 akiss
from bib_sepals import *
29 1 akiss
30 1 akiss
31 1 akiss
print('============== measure_sepal.py =================')
32 1 akiss
33 1 akiss
filename = argv[1]
34 1 akiss
outdir = argv[2]
35 4 akiss
36 1 akiss
outputfile=outdir+"/"+outdir+".csv"
37 1 akiss
38 4 akiss
outputdir=outdir+"/sections/"
39 4 akiss
projectiondir=outdir+"/projections"
40 1 akiss
41 1 akiss
print(argv)
42 1 akiss
43 1 akiss
imname = filename.split('/')[-1]
44 1 akiss
imnameroot = imname.split('.')[0]
45 1 akiss
46 4 akiss
projected2Dimage=projectiondir+"/"+imnameroot+"2D.png"
47 4 akiss
#projected2Dimage_hight=outdir+"/hight_"+imnameroot+"2D.png"
48 1 akiss
49 1 akiss
os.system("mkdir -p "+outputdir)
50 1 akiss
51 1 akiss
outname = outputfile.split('/')[-1]
52 1 akiss
outroot = outname.split('.')[0]
53 1 akiss
54 1 akiss
print("----------------------")
55 1 akiss
print("reading the file ", filename)
56 1 akiss
print("----------------------")
57 1 akiss
sepale=imread(filename)
58 1 akiss
voxelsize = sepale.voxelsize
59 1 akiss
print("voxelsize=",voxelsize)
60 1 akiss
s = sepale.shape
61 1 akiss
print("shape=",s)
62 1 akiss
63 1 akiss
print("Projecting on 2D...")
64 1 akiss
65 1 akiss
def project2D(im3):
66 1 akiss
        s2=(s[1],s[0])
67 1 akiss
        im2=np.zeros(s2, dtype='uint8')
68 1 akiss
        for i in range(s[0]):
69 1 akiss
                im2[:,i]=np.array([max(im3[i,j,:]) for j in range(s[1])])
70 1 akiss
        return im2
71 1 akiss
72 1 akiss
def project2D_hightmap(im3):
73 1 akiss
        s2=(s[1],s[0])
74 1 akiss
        im2=np.zeros(s2, dtype='uint8')
75 1 akiss
        for i in range(s[0]):
76 1 akiss
                im2[:,i]=np.array([list(im3[i,j,:]).index(max(im3[i,j,:])) for j in range(s[1])])
77 1 akiss
        return im2
78 1 akiss
79 1 akiss
80 1 akiss
sepale2D=project2D(sepale)
81 1 akiss
imsave2D(projected2Dimage, sepale2D)
82 1 akiss
83 1 akiss
'''
84 1 akiss
sepale2D=project2D_hightmap(sepale)
85 1 akiss
imsave2D(projected2Dimage_hight, sepale2D)
86 1 akiss
'''
87 1 akiss
88 1 akiss
print('Extracting principal section curves...')
89 1 akiss
90 1 akiss
# longitudinal section
91 1 akiss
x_section = sepale[int(s[0]/2), :, :]
92 1 akiss
xup_filename = outputdir+imnameroot+'_section_xup.txt'
93 1 akiss
xdown_filename = outputdir+imnameroot+'_section_xdown.txt'
94 1 akiss
x_filename = outputdir+imnameroot+'_section_x.txt'
95 1 akiss
x_pixelsize = (voxelsize[1],voxelsize[2])
96 1 akiss
xup, xdown = extract_curves(x_section)
97 1 akiss
98 1 akiss
write_curve(xup, x_pixelsize, xup_filename)
99 1 akiss
100 1 akiss
101 1 akiss
# transversal section
102 1 akiss
y_section = sepale[:,int(s[1]/2), :]
103 1 akiss
yup_filename = outputdir+imnameroot+'_section_yup.txt'
104 1 akiss
ydown_filename = outputdir+imnameroot+'_section_ydown.txt'
105 1 akiss
y_filename = outputdir+imnameroot+'_section_y.txt'
106 1 akiss
y_pixelsize = (voxelsize[0],voxelsize[2])
107 1 akiss
yup, ydown = extract_curves(y_section)
108 1 akiss
109 1 akiss
write_curve(yup, y_pixelsize, yup_filename)
110 1 akiss
111 1 akiss
112 1 akiss
#
113 1 akiss
# ----------------------------------------
114 1 akiss
print('Analyzing principal section curves...')
115 1 akiss
116 1 akiss
if not(os.path.isfile(outputfile)):
117 1 akiss
        FILE = open(outputfile,"w")
118 1 akiss
        FILE.write('Sepal;CurvedLength;FlatLength;LHeight;LR;LR21;LR22;LR41;LR423;LR44;CurvedWidth;FlatWidth;THeight;TR;TR21;TR22;TR41;TR423;TR44\n')
119 1 akiss
        FILE.close()
120 1 akiss
121 1 akiss
FILE = open(outputfile,"a")
122 1 akiss
FILE.write(imnameroot+';')
123 1 akiss
124 1 akiss
outfile=outroot+'_up_longitudinal.csv'
125 1 akiss
graph_name= outputdir+imnameroot+'_up_longitudinal.png'
126 1 akiss
length_x, flatlength_x, height_x, R, R21, R22, R41, R423, R44=analyze_curve1(imnameroot, xup, x_pixelsize, outfile, graph_name)
127 4 akiss
os.system("rm "+outfile)
128 1 akiss
129 1 akiss
FILE.write(str(length_x)+';'+str(flatlength_x)+';'+str(height_x)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+';')
130 1 akiss
131 1 akiss
outfile=outroot+'_up_transversal.csv'
132 1 akiss
graph_name= outputdir+imnameroot+'_up_transversal.png'
133 1 akiss
length_y, flatlength_y, height_y, R, R21, R22, R41, R423, R44=analyze_curve1(imnameroot, yup, y_pixelsize, outfile, graph_name)
134 4 akiss
os.system("rm "+outfile)
135 1 akiss
136 1 akiss
FILE.write(str(length_y)+';'+str(flatlength_y)+';'+str(height_y)+';'+str(R)+';'+str(R21)+';'+str(R22)+';'+str(R41)+';'+str(R423)+';'+str(R44)+'\n')