Commit dc83e35e authored by Jonathan Barés's avatar Jonathan Barés

Improvement of the contact detection and contact detection all along the root. Graphical output

parent f0b5d4c8
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Load libraries:
import numpy as np
import os
import glob
import time
import datetime
from PIL import Image
from scipy import misc
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from multiprocessing import Pool
from matplotlib.widgets import RectangleSelector
#Input parameter:
##Number of processor for post-processing parallelization:
NbProc=6
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Rename file and create time vectors:
print('... is sorting pictures')
##Get picutre names:
nam_pict=glob.glob('picture/*.jpg')
##Count the number of pictures in the folder to get the number of kept steps:
NbStp=len(nam_pict)
##Initialise time vector:
timePic=np.zeros(NbStp)
##Create picture vector:
null=os.system('mkdir pictureWh')
##Loop over picture to get their times:
for it in range(len(nam_pict)):
tmp=Image.open(nam_pict[it])._getexif()[36867]
dt_tmp=datetime.datetime.strptime(tmp, '%Y:%m:%d %H:%M:%S')
timePic[it]=time.mktime(dt_tmp.timetuple())
##Sort time vector convert in min and save it:
timePic-=np.min(timePic)
Is=np.argsort(timePic)
timePic=timePic[Is]
timePic/=60.
np.savetxt('time.txt',timePic,delimiter='\n')
##Loop over picture to sort and rename them:
for it in range(len(nam_pict)):
os.system('cp '+nam_pict[Is[it]]+' pictureWh/'+str(it).zfill(4)+'.jpg')
##Remove picture folder:
os.system('rm -rf picture')
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Rotate pictures:
print('... is rotating pictures')
##Function to rotate pictures:
def RotatePicture(iPct):
print(iPct)
null=os.system('convert pictureWh/'+'%04d'%iPct+'.jpg -rotate 270 pictureWh/'+'%04d'%iPct+'.jpg')
return 0
##Parallelization of the rotation:
p=Pool(NbProc)
null=p.map(RotatePicture,range(NbStp))
p.close()
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Reframe pictures:
##Measure picture size:
pict1=misc.imread('pictureWh/'+'%04d'%0+'.jpg')[:,:,1]
row,col=pict1.shape
##Select the ROI and compute FFT:
pict1=pict1[int(2*row/3):-1,:]
pict1FFT=np.fft.fft2(pict1)
##Function to measure the shift:
def MeasurePictureShift(iPct):
global row, col, pict1FFT
print(iPct)
###Load pictures an select ROI:
pict2=misc.imread('pictureWh/'+'%04d'%(iPct)+'.jpg')[:,:,1]
pict2=pict2[int(2*pict2.shape[0]/3):-1,:]
###Compute Fourier transform:
pict2FFT=np.conjugate(np.fft.fft2(pict2))
###Convolute:
pictCCor=np.real(np.fft.ifft2((pict1FFT*pict2FFT)))
###Compute the shift:
pictCCorShift=np.fft.fftshift(pictCCor)
yShift,xShift=np.unravel_index(np.argmax(pictCCorShift),(row,col))
yShift=yShift-int(row/2.)
xShift=xShift-int(col/2.)
a=np.argmax(pictCCorShift)
###return the shifts:
return xShift, yShift, a
##Parallelization of the shift computation:
print('... is centering pictures (measure displacement)')
p=Pool(NbProc)
xShiftVec0,yShiftVec0,test0=zip(*p.map(MeasurePictureShift,range(NbStp)))
xShiftVec0=np.array(xShiftVec0); yShiftVec0=np.array(yShiftVec0); test0=np.array(test0)
p.close()
##Detect picture problems and fix it:
###Peak a threshold on image correlation if needed:
plt.plot(np.abs(test0-np.mean(test0)),'o',markersize=5)
plt.xlabel('step')
plt.ylabel('correlation')
plt.title('Select a maximum acceptable value')
X=plt.ginput(1,timeout=-1)
plt.close()
II1=np.hstack((np.where(np.abs(test0-np.mean(test0))>X[0][1])[0],10000000))
###Peak a threshold on image shift if needed:
plt.plot(np.abs(xShiftVec0),'o',markersize=5)
plt.xlabel('step')
plt.ylabel('image shift')
plt.title('Select a maximum acceptable value')
X=plt.ginput(1,timeout=-1)
plt.close()
II2=np.hstack((np.where(np.abs(xShiftVec0)>X[0][1])[0],10000000))
II=np.unique(np.hstack((II1,II2)))
###Solve problems if necessary:
if len(II)>1:
print('Problem with following pictures: '+str(II)[1:-2])
III=np.hstack((np.array([-1]),np.where(np.diff(II)>1)[0]))
for it_III in range(len(III)-1):
it_ref=II[III[it_III]+1]
I_mod=np.array(range(it_ref+1,II[III[it_III+1]]+1))
for it_mod in I_mod:
os.system('rm -rf pictureWh/'+'%04d'%it_mod+'.jpg')
os.system('cp pictureWh/'+'%04d'%it_ref+'.jpg pictureWh/'+'%04d'%it_mod+'.jpg')
os.system('rm -rf picturePl/'+'%04d'%it_mod+'.jpg')
os.system('cp picturePl/'+'%04d'%it_ref+'.jpg picturePl/'+'%04d'%it_mod+'.jpg')
for it_II in II[0:-1]:
A0=MeasurePictureShift(it_II)
xShiftVec0[it_II]=A0[0]
yShiftVec0[it_II]=A0[1]
test0[it_II]=A0[2]
yShiftVec0=-yShiftVec0
##Computation of the new image size:
MaxSft_x=np.max(xShiftVec0)+1
MinSft_x=np.min(xShiftVec0)-1
MaxSft_y=np.max(yShiftVec0)+1
MinSft_y=np.min(yShiftVec0)-1
colN=int(col-MaxSft_x+MinSft_x)
rowN=int(row-MaxSft_y+MinSft_y)
##Rescale vectors:
xShiftVecN=(MaxSft_x-xShiftVec0).astype(int)
yShiftVecN=(MaxSft_y-yShiftVec0).astype(int)
##Function to reframe the pictures:
def ReframePicture(iPct):
print(iPct)
###Load pictures:
pict1=misc.imread('pictureWh/'+'%04d'%(iPct)+'.jpg')
###Crop pictures:
pict2=pict1[:,xShiftVecN[iPct]:xShiftVecN[iPct]+colN,:]
pict3=pict2[yShiftVecN[iPct]:yShiftVecN[iPct]+rowN,:,:]
###Save picture:
misc.imsave('pictureWh/'+'%04d'%iPct+'.jpg',pict3)
##Parallelization of the cropping:
print('... is centering pictures (crop)')
p=Pool(NbProc)
null=p.map(ReframePicture,range(NbStp-1))
p.close()
##Select the area of interest:
###Open the last picture:
pict1=mpimg.imread('pictureWh/'+'%04d'%(NbStp-1)+'.jpg')
###Define selection functions:
def line_select_callback(eclick,erelease):
global x1, y1, x2, y2
x1,y1=eclick.xdata,eclick.ydata
x2,y2=erelease.xdata,erelease.ydata
def toggle_selector(event):
print(' Key pressed.')
if event.key in ['Q', 'q'] and toggle_selector.RS.active:
print(' RectangleSelector deactivated.')
toggle_selector.RS.set_active(False)
if event.key in ['A', 'a'] and not toggle_selector.RS.active:
print(' RectangleSelector activated.')
toggle_selector.RS.set_active(True)
###Select:
fig,current_ax=plt.subplots()
plt.imshow(pict1,cmap=plt.cm.Greys)
plt.title('select area to consider, close to validate')
toggle_selector.RS=RectangleSelector(current_ax,line_select_callback,drawtype='box',useblit=True,button=[1,3],minspanx=5,minspany=5,spancoords='pixels',interactive=True)
plt.connect('key_press_event',toggle_selector)
plt.show()
plt.close()
##Function to crop the pictures:
def CropPicture(iPct):
global x1, y1, x2, y2
print(iPct)
###Load pictures:
pict1=misc.imread('pictureWh/'+'%04d'%iPct+'.jpg')
###Crop pictures:
pictF1=pict1[int(y1):int(y2),int(x1):int(x2),:]
###Save picture:
misc.imsave('pictureWh/'+'%04d'%iPct+'.jpg',pictF1)
##Parallelization of the cropping:
print('... is cropping pictures')
p=Pool(NbProc)
null=p.map(CropPicture,range(NbStp))
p.close()
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Create movies:
##Create picture vector:
null=os.system('mkdir tmpWh')
##function to plot the pictures with the time:
def PlotPicture(iPct):
global timePic
print(iPct)
###Load the time:
t0=timePic[iPct]
h0,m0=divmod(t0,60)
j0,h0=divmod(h0,24)
###Load pictures:
pict1=mpimg.imread('pictureWh/'+'%04d'%iPct+'.jpg')
###Plot picture:
plt.imshow(pict1)
plt.axis('off')
plt.title('%02d'%j0+'days '+'%02d'%h0+'hours - '+str(iPct))
plt.savefig('tmpWh/'+'%04d'%iPct+'.png',dpi=200)
plt.close()
##Parallelization of the cropping:
print('... is plotting pictures')
p=Pool(1)
null=p.map(PlotPicture,range(NbStp))
p.close()
##Make movies:
print('... is making the movie')
null=os.system('ffmpeg -r 15 -y -f image2 -i tmpWh/%04d.png -qscale 1 Wh.avi')
##Remove folders:
null=os.system('rm -f -R tmpWh')
......@@ -19,6 +19,9 @@ from cv2 import Canny
cnt_th=12
###Radius area considered to compute contact direction [px]:
rad_px=15
###Contact detection threshold for root back [px]:
cnt_th_b=8
##Count the number of pictures in the folder to get the number of keeped steps:
NbStp=len(fnmatch.filter(os.listdir('pictureWh'),'*.jpg'))
......@@ -104,7 +107,6 @@ for it_l in range(len(obst_list)):
np.savetxt('obstacle/coordinate.txt',np.vstack(mat_save))
##Make obstacle mask:
###Initialise mask:
obst_mask=np.zeros((pict0.shape[0],pict0.shape[1]))
......@@ -137,12 +139,11 @@ misc.imsave('obstacle/edgeMask.png',obst_edge)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Detect contact with obstacles:
#Detect contact between main root tip and obstacles:
##Load obstacle:
obst_edge=misc.imread('obstacle/edgeMask.png')
np.sort(np.unique(obst_edge))[1]
obst_edge=obst_edge/np.sort(np.unique(obst_edge))[1]
obst_edge=np.round(obst_edge/np.sort(np.unique(obst_edge))[1])
##Load time and position of the tip:
pos_tip=np.loadtxt('mainRoot/rootTip/position.txt')
......@@ -152,6 +153,9 @@ time_tip=np.loadtxt('mainRoot/rootTip/time.txt')
dir0=np.loadtxt('mainRoot/direction.txt')
time0=np.loadtxt('time.txt')
##Load image edges:
IJ_edge=np.loadtxt('pictureSkeleton/IJ_edge.txt')
##Interpolate direction:
dir_tip=np.interp(time_tip,time0,dir0)
......@@ -161,6 +165,8 @@ cnt_num=np.zeros(len(pos_tip))
cnt_dir=-1.*np.ones(len(pos_tip))
###Prepare obstacle:
I_o,J_o=np.where(obst_edge>0)
J_o=J_o-IJ_edge[2]
I_o=I_o-IJ_edge[0]
N_o=obst_edge[I_o,J_o]
###Loop over positions:
for it_t in range(len(time_tip)):
......@@ -199,26 +205,131 @@ for it_t in range(len(time_tip)):
#~ plt.savefig('tmp/'+str(it_t).zfill(5)+'.png')
#~ plt.close()
###Save contact information:
np.savetxt('mainRoot/rootTip/contactObstacle.txt',cnt_num)
np.savetxt('mainRoot/rootTip/contactAngle.txt',cnt_dir)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Detect contact between main root body and obstacles:
##Load obstacle mask:
obst_edge=misc.imread('obstacle/edgeMask.png')
obst_edge=np.round(obst_edge/np.sort(np.unique(obst_edge))[1])
##Load time and position of the tip:
time0=np.loadtxt('time.txt')
##Load image edges:
IJ_edge=np.loadtxt('pictureSkeleton/IJ_edge.txt')
##Make storage folder:
os.system('mkdir mainRoot/rootContact')
##Get contact and obstacle number:
###Prepare obstacle:
I_o,J_o=np.where(obst_edge>0)
N_o=obst_edge[I_o,J_o]
J_o=J_o-IJ_edge[2]
I_o=I_o-IJ_edge[0]
###Loop over positions:
for it_t in range(len(time0)):
####Get current root path:
Ic_Jc=np.loadtxt('mainRoot/rootPath/'+str(it_t).zfill(4)+'.txt')
Ic=Ic_Jc[:,0]; Jc=Ic_Jc[:,1]
####Get current weighted skeleton:
skl_cur=misc.imread('pictureSkeleton/weightedSkeleton/'+str(it_t).zfill(4)+'.png')
skl_cur=skl_cur/np.sort(np.unique(skl_cur))[1]
##~DEBUG~##
# plt.imshow(skl_cur)
# plt.plot(Jc,Ic,'-r')
# plt.show()
# plt.close()
####Loop along path to detect contact:
cnt_num=np.zeros(len(Ic))
for it in range(len(Ic)):
#####Measure distance to obstacles:
d_cur=np.sqrt((Ic[it]-I_o)**2.+(Jc[it]-J_o)**2.)
#####If close enough consider the contact:
##~DEBUG~##
# plt.imshow(skl_cur)
# plt.plot(Jc[it],Ic[it],'or')
# plt.plot(J_o,I_o,'.r')
# plt.show()
# plt.close()
if np.min(d_cur)<(skl_cur[int(Ic[it]),int(Jc[it])]+cnt_th_b):
II=np.where(d_cur==np.min(d_cur))[0][0]
cnt_num[it]=N_o[II]
####Save the result:
np.savetxt('mainRoot/rootContact/'+str(it_t).zfill(4)+'.txt',cnt_num)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Display results:
##Make temporary storage:
os.system('mkdir mainRoot/tmp_obst')
##Load root tip contact and contact direction:
cnt_num_0=np.loadtxt('mainRoot/rootTip/contactObstacle.txt')
time_0=np.loadtxt('mainRoot/rootTip/time.txt')
cnt_num0=np.interp(time0,time_0,cnt_num_0)
#Load root tip directon:
dir0=np.loadtxt('mainRoot/direction.txt')
##Loop over time step:
for it_t in range(len(time0)):
###Load picture:
pict0=misc.imread('pictureWh/'+str(it_t).zfill(4)+'.jpg')
###Load root path:
Ic_Jc=np.loadtxt('mainRoot/rootPath/'+str(it_t).zfill(4)+'.txt')
Ic=Ic_Jc[:,0]+IJ_edge[0]; Jc=Ic_Jc[:,1]+IJ_edge[2]
###Load contacts:
cnt_num=np.loadtxt('mainRoot/rootContact/'+str(it_t).zfill(4)+'.txt')
###Plot background:
plt.imshow(pict0)
###Plot obstacles:
for it_l in range(len(obst_list)):
plt.plot(obst_list[it_l][:,0],obst_list[it_l][:,1],'-',linewidth=2,color=plt.cm.jet(float(it_l)/len(obst_list)))
###Plot path:
plt.plot(Jc,Ic,'-w',linewidth=2)
###Plot contacts:
cnt_num_u=np.unique(cnt_num)
for it_o in range(len(cnt_num_u)-1):
N_o_c=int(cnt_num_u[it_o+1])
II=np.where(cnt_num==N_o_c)[0]
plt.plot(Jc[II],Ic[II],'.',markersize=5,color=plt.cm.jet(float(N_o_c-1)/len(obst_list)))
###Plot tip and contact orientation:
if cnt_num0[it_t]>0:
Jc=Jc[0]
Ic=Ic[0]
Dc=-dir0[it_t]
plt.plot(Jc,Ic,'+r',color=plt.cm.jet(float(cnt_num0[it_t]-1)/len(obst_list)))
nn=300.
plt.plot(np.array([Jc,Jc+nn*m.cos(Dc)]),np.array([Ic,Ic+nn*m.sin(Dc)]),'-',linewidth=2,color=plt.cm.jet(float(cnt_num0[it_t]-1)/len(obst_list)))
plt.axis('equal')
plt.axis('off')
t0=time0[it_t]
h0=int(t0/60.)
m0=int(t0-60.*h0)
plt.title(str(h0).zfill(2)+' : '+str(m0).zfill(2)+' - '+str(it_t).zfill(4))
plt.savefig('mainRoot/tmp_obst/'+str(it_t).zfill(3)+'.png',bbox_inches='tight',dpi=200)
# plt.show()
plt.close()
##Make movies:
null=os.system('ffmpeg -r 10 -y -f image2 -i mainRoot/tmp_obst/%03d.png -qscale 1 mainRoot/MainRootContact.avi')
##Remove folders:
null=os.system('rm -f -R mainRoot/tmp_obst')
problem affichage temps (en rouge) dans le sortie graphique combinée
sortie graphique pour les contacts contact de 'l'arrière' de racine
Faire tourner sur un set d'images test
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