Commit f0b5d4c8 authored by Jonathan Barés's avatar Jonathan Barés
Browse files

Code is converted to work with Python3.

In 03 a piece of code is added to compute the local orientation of the root.
parent a9c0950a
......@@ -18,9 +18,6 @@ Nl=965
##Number of processor for post-processing parallelization:
NbProc=5
#Define global variables:
global row, col, timePic
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Rename file and create time vectors:
......@@ -36,7 +33,7 @@ null=os.system('mkdir picturePl')
null=os.system('mkdir pictureWh')
##Loop over picture names:
cnt=0
for it1 in xrange(Ni,Nl+1):
for it1 in range(Ni,Nl+1):
try:
###Find picture names:
namCurPl=glob.glob('%04d'%it1+'_*_Pl.jpg')[0]
......@@ -61,14 +58,14 @@ print('... is rotating pictures')
##Function to rotate pictures:
def RotatePicture(iPct):
print iPct
print(iPct)
null=os.system('convert picturePl/'+'%04d'%iPct+'.jpg -rotate 90 picturePl/'+'%04d'%iPct+'.jpg')
null=os.system('convert pictureWh/'+'%04d'%iPct+'.jpg -rotate 90 pictureWh/'+'%04d'%iPct+'.jpg')
return 0
##Parallelization of the rotation:
p=Pool(NbProc)
null=p.map(RotatePicture,xrange(NbStp))
null=p.map(RotatePicture,range(NbStp))
p.close()
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
......@@ -81,7 +78,7 @@ row,col=pict1.shape
##Function to measure the shift:
def MeasurePictureShift(iPct):
global row, col
print iPct
print(iPct)
###Load pictures:
pict1=misc.imread('picturePl/'+'%04d'%iPct+'.jpg')[:,:,1]
pict2=misc.imread('picturePl/'+'%04d'%(iPct+1)+'.jpg')[:,:,1]
......@@ -104,7 +101,7 @@ def MeasurePictureShift(iPct):
##Parallelization of the shift computation:
print('... is centering pictures (measure displacement)')
p=Pool(NbProc)
xShiftVec0,yShiftVec0,test0=zip(*p.map(MeasurePictureShift,xrange(NbStp-1)))
xShiftVec0,yShiftVec0,test0=zip(*p.map(MeasurePictureShift,range(NbStp-1)))
xShiftVec0=np.array(xShiftVec0); yShiftVec0=np.array(yShiftVec0); test0=np.array(test0)
p.close()
......@@ -150,7 +147,7 @@ if len(II)>1:
##Add the shifts:
xShiftVec=np.zeros(NbStp)
yShiftVec=np.zeros(NbStp)
for it in xrange(1,NbStp):
for it in range(1,NbStp):
xShiftVec[it]=xShiftVec[it-1]+xShiftVec0[it-1]
yShiftVec[it]=yShiftVec[it-1]+yShiftVec0[it-1]
......@@ -165,7 +162,7 @@ xShiftVecN=(MaxSft-xShiftVec).astype(int)
##Function to reframe the pictures:
def ReframePicture(iPct):
print iPct
print(iPct)
###Load pictures:
pict0=misc.imread('picturePl/'+'%04d'%iPct+'.jpg')
pict1=misc.imread('pictureWh/'+'%04d'%iPct+'.jpg')
......@@ -179,12 +176,10 @@ def ReframePicture(iPct):
##Parallelization of the cropping:
print('... is centering pictures (crop)')
p=Pool(NbProc)
null=p.map(ReframePicture,xrange(NbStp))
null=p.map(ReframePicture,range(NbStp))
p.close()
##Select the area of interest:
###Waiting:
print raw_input('Is waiting. Press ENTER to move...')
###Open the last picture:
pict1=mpimg.imread('pictureWh/'+'%04d'%(NbStp-1)+'.jpg')
###Define selection functions:
......@@ -215,7 +210,7 @@ plt.close()
##Function to crop the pictures:
def CropPicture(iPct):
global x1, y1, x2, y2
print iPct
print(iPct)
###Load pictures:
pict0=misc.imread('picturePl/'+'%04d'%iPct+'.jpg')
pict1=misc.imread('pictureWh/'+'%04d'%iPct+'.jpg')
......@@ -229,7 +224,7 @@ def CropPicture(iPct):
##Parallelization of the cropping:
print('... is cropping pictures')
p=Pool(NbProc)
null=p.map(CropPicture,xrange(NbStp))
null=p.map(CropPicture,range(NbStp))
p.close()
......@@ -243,7 +238,7 @@ null=os.system('mkdir tmpWh')
##function to plot the pictures with the time:
def PlotPicture(iPct):
global timePic
print iPct
print(iPct)
###Load the time:
t0=timePic[iPct]
h0,m0=divmod(t0,60)
......@@ -266,7 +261,7 @@ def PlotPicture(iPct):
##Parallelization of the cropping:
print('... is plotting pictures')
p=Pool(1)
null=p.map(PlotPicture,xrange(NbStp))
null=p.map(PlotPicture,range(NbStp))
p.close()
##Make movies:
......
......@@ -134,7 +134,6 @@ print('Is extracting root mask...')
os.system('mkdir pictureSkeleton/mask')
###Loop over the pictures:
for iPct in range(NbStp):
global Imin,Imax,Jmin,Jmax,timePic,JIseed,ThrdPict,pict_mask
print(str(iPct)+' / '+str(NbStp))
####Load pictures:
pictA=misc.imread('pictureWh/'+'%04d'%iPct+'.jpg')[:,:,0]
......@@ -162,7 +161,7 @@ for iPct in range(NbStp):
#####Get current pixels
I=np.where(pictA==itArea)
#####If current area is big enough to be considered:
if (len(I[0])>1000):
if (len(I[0])>500):
#####If the seed point is in the current area this is the root:
if np.min(np.sqrt((I[0]-JIseed[0][1])**2.+(I[1]-JIseed[0][0])**2.))<2:
flag1=0
......@@ -206,7 +205,7 @@ os.system('mkdir rootTip/byStep')
timePic=np.loadtxt('time.txt')
##Loop over the steps:
for iPct in xrange(NbStp):
for iPct in range(NbStp):
###Load the mask:
pict4=misc.imread('pictureSkeleton/mask/'+'%04d'%iPct+'.png')
###Computation of the mask contour:
......@@ -230,7 +229,7 @@ for iPct in xrange(NbStp):
####Computation of the cross-product:
Curv=np.zeros(len(contour))
for itCurv in xrange(1,len(contour)-1):
for itCurv in range(1,len(contour)-1):
x0=X0[itCurv-1]; y0=Y0[itCurv-1]
x1=X0[itCurv]; y1=Y0[itCurv]
x2=X0[itCurv+1]; y2=Y0[itCurv+1]
......@@ -289,7 +288,7 @@ for iPct in xrange(NbStp):
##Plot and save root tip pathes in color time evolution:
###Loop over the steps:
X=[]; Y=[]; T=[]
for iPct in xrange(NbStp):
for iPct in range(NbStp):
if os.path.isfile('rootTip/byStep/'+'%04d'%iPct+'.txt'):
XY_tmp=np.loadtxt('rootTip/byStep/'+'%04d'%iPct+'.txt')
if XY_tmp.ndim>1:
......@@ -324,7 +323,7 @@ plt.close()
os.system('mkdir rootTip/byTip')
###Loop over the steps
for iPct in xrange(NbStp):
for iPct in range(NbStp):
###Load current tips position vector:
A=np.loadtxt('rootTip/byStep/'+'%04d'%iPct+'.txt')
if (A.ndim>0):
......@@ -340,7 +339,7 @@ for iPct in xrange(NbStp):
tc=timePic[iPct]
if (iPct==0):
###Initialize storage matrix:
for iPt in xrange(len(Xc)):
for iPt in range(len(Xc)):
if (iPt==0):
MatRoot=[np.array([Xc[iPt],Yc[iPt],tc])]
else:
......@@ -354,7 +353,7 @@ for iPct in xrange(NbStp):
Xc=np.delete(Xc,0); Yc=np.delete(Yc,0)
####Measure the distance to the previously detected tips:
Dc=np.zeros(len(MatRoot))
for iPt in xrange(len(MatRoot)):
for iPt in range(len(MatRoot)):
#####Load the last position:
if (MatRoot[iPt].ndim==1):
xcc=MatRoot[iPt][0]; ycc=MatRoot[iPt][1]
......@@ -376,7 +375,7 @@ for iPct in xrange(NbStp):
MatRoot0=MatRoot
MatRoot=[]
###Loop over detected roots:
for iMrt in xrange(len(MatRoot0)):
for iMrt in range(len(MatRoot0)):
if (MatRoot0[iMrt].ndim>1):
####Load the current root:
curRoot=MatRoot0[iMrt]
......@@ -390,7 +389,7 @@ for iMrt in xrange(len(MatRoot0)):
MatRoot=MatRoot+[MatRoot0[iMrt]]
# ! Debug !:
#~ for iMrt in xrange(len(MatRoot)):
#~ for iMrt in range(len(MatRoot)):
#~ curRoot=MatRoot0[iMrt]
#~ plt.plot(curRoot[:,0],curRoot[:,1],'-')
......@@ -401,7 +400,7 @@ for iMrt in xrange(len(MatRoot0)):
###Prepare to store the speed amplitude:
Vmin=float("inf"); Vmax=0.
###Loop over root trajectories:
for iMrt in xrange(len(MatRoot)):
for iMrt in range(len(MatRoot)):
####Get the current trajectory:
Xc=MatRoot[iMrt][:,0]; Yc=MatRoot[iMrt][:,1]; Tc=MatRoot[iMrt][:,2]
####Compute the speed:
......@@ -423,7 +422,7 @@ for iMrt in xrange(len(MatRoot)):
##Storage of the data:
for iMrt in xrange(len(MatRoot)):
for iMrt in range(len(MatRoot)):
np.savetxt('rootTip/byTip/'+'%04d'%iMrt+'.txt',MatRoot[iMrt],delimiter=' ')
##Plot by speed:
......@@ -431,7 +430,7 @@ for iMrt in xrange(len(MatRoot)):
pict4=misc.imread('pictureSkeleton/mask/'+'%04d'%(NbStp-1)+'.png')
plt.imshow(pict4,plt.cm.Greys)
###Loop over roots:
for iMrt in xrange(len(MatRoot)):
for iMrt in range(len(MatRoot)):
Xc=MatRoot[iMrt][:,0]; Yc=MatRoot[iMrt][:,1]; Tc=MatRoot[iMrt][:,2]; Vc=MatRoot[iMrt][:,3]
if iMrt==0:
sc=plt.scatter(Xc,Yc,c=Vc,s=5,cmap=plt.cm.get_cmap('jet'))
......@@ -454,7 +453,7 @@ os.system('mkdir tmp')
def PlotSkeleton(iPct):
global timePic
print iPct
print(iPct)
###Load the time:
t0=timePic[iPct]
h0,m0=divmod(t0,60)
......@@ -470,7 +469,7 @@ def PlotSkeleton(iPct):
##Cropping:
print('... is plotting root skeleton')
for itStp in xrange(NbStp):
for itStp in range(NbStp):
PlotSkeleton(itStp)
##Make movie:
......@@ -485,7 +484,7 @@ null=os.system('rm -f -R tmp')
os.system('mkdir tmp')
fig = plt.figure()
ax = fig.gca(projection='3d')
for iMrt in xrange(len(MatRoot)):
for iMrt in range(len(MatRoot)):
Xc=MatRoot[iMrt][:,0]; Yc=MatRoot[iMrt][:,1]; Tc=MatRoot[iMrt][:,2];
ax.plot(Tc/(24*60),Xc,-Yc,'k-',linewidth=3)
......@@ -494,7 +493,7 @@ ax.set_xlabel('time (d)')
ax.set_ylabel('x (px)')
ax.set_zlabel('y (px)')
ang=np.linspace(70,70+360,360)
for iPic in xrange(len(ang)):
for iPic in range(len(ang)):
ax.view_init(5,ang[iPic])
plt.savefig('tmp/'+'%04d'%iPic+'.png',dpi=150)
......
......@@ -45,7 +45,7 @@ while flag0:
X0=np.vstack(MatRoot)[:,0]
Y0=-np.vstack(MatRoot)[:,1]
####Plot them
for iMrt in xrange(len(MatRoot)):
for iMrt in range(len(MatRoot)):
Xc=MatRoot[iMrt][:,0]; Yc=MatRoot[iMrt][:,1]
plt.plot(Xc,-Yc,'-',markersize=5,linewidth=2)
......@@ -126,7 +126,7 @@ for NU in range(3):
#~ plt.show()
#~ plt.close()
##Compute the speed from the smothed signal: !!! attention certain diff T nuls a verifier!!!
##Compute the speed from the smothed signal: !!! attention parfois certain diff T nuls a verifier!!!
###Direct computation:
V=np.zeros(len(T))
V[1:]=np.sqrt((P[1:,0]-P[0:-1,0])**2.+(P[1:,1]-P[0:-1,1])**2.)/(T[1:]-T[0:-1])
......@@ -140,12 +140,19 @@ Vn=np.convolve(V,np.ones(20)/20.,mode='same')
#~ plt.show()
#~ plt.close()
##Plot result:
plt.plot(T,Vn,'-k',linewidth=3)
plt.xlabel('time (min)')
plt.ylabel('tip speed (px/min)')
plt.savefig('mainRoot/rootTip/speedTip.png',dpi=400)
plt.savefig('mainRoot/rootTip/speedTip.svg')
##Saved data:
os.system('mkdir mainRoot')
os.system('mkdir mainRoot/rootTip')
np.savetxt('mainRoot/rootTip/position.txt',P)
np.savetxt('mainRoot/rootTip/time.txt',T)
np.savetxt('mainRoot/rootTip/speedTip.txt',V)
np.savetxt('mainRoot/rootTip/speedTip.txt',Vn)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
......@@ -181,7 +188,7 @@ def GetWeightSkel(iPct):
###Parallelization of the skeleton:
p=Pool(NbProc)
null=p.map(GetWeightSkel,xrange(NbStp))
null=p.map(GetWeightSkel,range(NbStp))
p.close()
##Function to follow the main root:
......@@ -374,7 +381,6 @@ a=plt.ginput(1,timeout=-1)
plt.close()
S_max=int(a[0][0])
##Fix pathes:
###Loop over root path:
for iPct in range(1,NbStp):
......@@ -421,6 +427,8 @@ for iPct in range(NbStp):
sc=plt.scatter(np.zeros(len(time0)),np.zeros(len(time0)),c=time0,s=0,cmap=cm)
plt.colorbar(sc)
plt.axis('equal')
plt.xlabel('px')
plt.ylabel('px')
plt.savefig('mainRoot/mainRoot.png',dpi=400)
plt.savefig('mainRoot/mainRoot.svg')
plt.show()
......@@ -478,15 +486,15 @@ vec_V[1:]=(vec_L_smth[1:]-vec_L_smth[0:-1])/(time0[1:]-time0[0:-1])
plt.plot(time0,vec_V,'-k',linewidth=3)
plt.xlabel('time (min)')
plt.ylabel('main root growth speed (px/min)')
plt.savefig('mainRoot/mainRootSpeed.png',dpi=400)
plt.savefig('mainRoot/mainRootSpeed.svg')
plt.savefig('mainRoot/mainRootLengthSpeed.png',dpi=400)
plt.savefig('mainRoot/mainRootLengthSpeed.svg')
plt.show()
plt.close()
###Save data:
np.savetxt('mainRoot/speedLength.txt',vec_V)
##Computation of the main root curvature evolution:
###Function to compute curvature:
##Computation of the main root curvature and orientation evolution:
###Function to compute curvature and local orientation:
def GetCurvature(iPct):
print(str(iPct)+' / '+str(NbStp))
####Load root path coordinates:
......@@ -525,12 +533,20 @@ def GetCurvature(iPct):
####Store curvature:
np.savetxt('mainRoot/rootCurvature/'+str(iPct).zfill(4)+'.txt',C0)
####Local orientation computation:
O0=np.arctan2(J1(s[1:])-J1(s[0:-1]),-I1(s[1:])+I1(s[0:-1]))*180./m.pi
O0=np.append(O0,O0[-1])
####Store curvature:
np.savetxt('mainRoot/rootOrientation/'+str(iPct).zfill(4)+'.txt',O0)
###Make storage folder:
os.system('mkdir mainRoot/rootCurvature')
os.system('mkdir mainRoot/rootOrientation')
###Parallelization of the main root curvature computation:
p=Pool(NbProc)
null=p.map(GetCurvature,xrange(NbStp))
null=p.map(GetCurvature,range(NbStp))
p.close()
###Display the curvature:
......@@ -558,8 +574,9 @@ plt.imshow(matDisp,vmin=Cm,vmax=CM,cmap=plt.cm.get_cmap('jet'),aspect='auto')
plt.colorbar()
plt.title('root curvature (1/px)')
plt.xlabel('curvilinear coordinate (px)')
plt.ylabel('time')
ax.set_yticks([])
plt.ylabel('time (min)')
ax.set_yticks([0,0.2*matDisp.shape[0],0.4*matDisp.shape[0],0.6*matDisp.shape[0],0.8*matDisp.shape[0],matDisp.shape[0]])
ax.set_yticklabels([str(int(time0[-1])),str(int(time0[int(len(time0)*0.8)])),str(int(time0[int(len(time0)*0.6)])),str(int(time0[int(len(time0)*0.4)])),str(int(time0[int(len(time0)*0.2)])),'0'])
plt.savefig('mainRoot/mainRootCurvature1.png',dpi=400)
plt.savefig('mainRoot/mainRootCurvature1.svg')
plt.show()
......@@ -579,11 +596,63 @@ for iPct in range(NbStp):
plt.title('root curvature magnitude (1/px)')
plt.axis('equal')
plt.xlabel('px')
plt.ylabel('px')
plt.savefig('mainRoot/mainRootCurvature2.png',dpi=400)
#~ plt.savefig('mainRoot/mainRootCurvature2.svg')
plt.show()
plt.close()
###Display the local orientation:
####Color amplitude:
N0=0
for it in range(NbStp):
n0=len(np.loadtxt('mainRoot/rootOrientation/'+str(it).zfill(4)+'.txt'))
if N0<n0:
N0=n0
####In matrix form:
matDisp=float('nan')*np.ones((len(time0),N0))
for iPct in range(NbStp):
if os.path.isfile('mainRoot/rootOrientation/'+str(iPct).zfill(4)+'.txt'):
O0=np.loadtxt('mainRoot/rootOrientation/'+str(iPct).zfill(4)+'.txt')
O0=O0[np.arange(len(O0)-1,0,-1)]
matDisp[-1-iPct,0:len(O0)]=O0
fig, ax = plt.subplots()
plt.imshow(matDisp,vmin=np.nanmin(matDisp),vmax=np.nanmax(matDisp),cmap=plt.cm.get_cmap('jet'),aspect='auto')
plt.colorbar()
plt.title('root orientation (degree)')
plt.xlabel('curvilinear coordinate (px)')
plt.ylabel('time (min)')
ax.set_yticks([0,0.2*matDisp.shape[0],0.4*matDisp.shape[0],0.6*matDisp.shape[0],0.8*matDisp.shape[0],matDisp.shape[0]])
ax.set_yticklabels([str(int(time0[-1])),str(int(time0[int(len(time0)*0.8)])),str(int(time0[int(len(time0)*0.6)])),str(int(time0[int(len(time0)*0.4)])),str(int(time0[int(len(time0)*0.2)])),'0'])
plt.savefig('mainRoot/mainRootOrientation1.png',dpi=400)
plt.savefig('mainRoot/mainRootOrientation1.svg')
plt.show()
plt.close()
###Along the root:
for iPct in range(NbStp):
IJ=np.loadtxt('mainRoot/rootPath/'+str(iPct).zfill(4)+'.txt')
II=np.where(IJ[:,0]>Imax)[0]
I0=-IJ[II,0]
J0=IJ[II,1]
O0=np.loadtxt('mainRoot/rootOrientation/'+str(iPct).zfill(4)+'.txt')
if iPct==0:
sc=plt.scatter(J0,I0,c=O0,s=0.03,cmap=plt.cm.get_cmap('jet'),vmin=np.nanmin(matDisp),vmax=np.nanmax(matDisp))
plt.colorbar(sc)
else:
plt.scatter(J0,I0,c=O0,s=0.03,cmap=plt.cm.get_cmap('jet'),vmin=np.nanmin(matDisp),vmax=np.nanmax(matDisp))
plt.title('root orientation (degree)')
plt.axis('equal')
plt.xlabel('px')
plt.ylabel('px')
plt.savefig('mainRoot/mainRootOrientation2.png',dpi=400)
#~ plt.savefig('mainRoot/mainRootOrientation2.svg')
plt.show()
plt.close()
##Computation of the main root tortuosity evolution:
###Make storage vectors:
vec_Tor1=np.zeros(len(vec_L))
......@@ -665,7 +734,7 @@ def FractalDim(iPct):
###Parallelization of the main root fractal dimension computation:
p=Pool(NbProc)
vec_FD=p.map(FractalDim,xrange(NbStp))
vec_FD=p.map(FractalDim,range(NbStp))
p.close()
###Display the evolution of the main root tortuosities:
plt.plot(time0,vec_FD,'-k',linewidth=3)
......
comme C sortir l'orientation locale de la racine principale
axe vertical de la sortie graphique de de C sous forme de matrice, mettre le temps.
problem affichage temps (en rouge) dans le sortie graphique combinée
remonter dans le chapeau des codes toutes les variables
sortie graphique pour les contacts contact de l'arrière de racine
sortie graphique v_length et v_tip
problem mesure de C bruité sur la dernière manip (comprendre)
sortie graphique pour les contacts contact de 'l'arrière' de racine
Utiliser un widget pour le choix de la valeur de seuillage
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment