Commit 5ec92d10 authored by Jonathan Barés's avatar Jonathan Barés

Fully updated code adding inputs (stored in yaml) in 00_, fixing bugs,...

Fully updated code adding inputs (stored in yaml) in 00_, fixing bugs, improving graphical interface, improving method, adding outputs
parent 18d6f1b9
# This code permits to give the inputs of the image post processing
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
# Load libraries:
import yaml
import io
from glob import glob
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
# Set inputs:
## File organization:
### Name of the image folder [string]:
nam_img = 'picture'
### Image extension [string]:
ext_img = '.jpg'
### Number step to treat [int]:
nb_stp = len(glob(nam_img+'/*'+ext_img))
## Image preprocessing:
### Image rotation [int, degree]:
img_rot = 270
### Need to correct picture shift [bool]:
corr_shift = True
## Displaying:
### Number of frame per second for movies [int]:
fps_mov = 15
### Amplitude of the root tip direction vector [int, px]:
amp_dir = 200
## Root tip detection and tracking:
### Dilation of the obstacle mask for obstcale removing [odd int, px]:
dil_msk = 7
### Proximity threshold for root tip connection [int, px]:
prox_conn_th = 70
### Proximity threshold for root tip detection [int, px]:
prox_detec_th = 15
### Minimum area of the root [int, px^2]:
min_area_root = 40
### Root icture smothening intensity [int, px]:
pict_smth = 6
### Root contour smothening intensity [int, px]:
cont_smth = 4
### Root contour curvature smoothening [int, step]:
curv_smth = 10
### Threshold on the curvature for root tips detection [float,<0]:
curv_th = -0.12
### Root tip speed smoothening intensity [int, step]:
spd_smth = 5
## Main root tip analysis:
### Degree of the Chaikin's corner smoothening of the main root tip trajectory [int]:
traj_smth = 3
### Main root tip speed smothening intensity [int, step]:
main_spd_smth = 10
### Mask smoothing dimensions [odd int, px]:
msk_smth = 11
### Root length smoothing intensity [int, step]:
lgth_smth = 10
### Targetted error for the polynomial fitting [float]:
err_pol = 0.2
### Distance over wich to compute root tip direction [int, px]:
d_dir = 20
## Contact detection intputs:
### Contact detection threshold [int, px]:
cnt_th = 12
### Radius area considered to compute contact direction [int, px]:
cnt_rad = 15
### Contact detection threshold for root back [int, px]:
cnt_th_b = 8
## Computational:
### Number of processor for post-processing parallelization:
nb_proc = 6
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Save inputs:
## Make the dictionnary:
data_input={
'general': { 'image folder name': nam_img ,
'image extension': ext_img ,
'number of step': nb_stp} ,
'image preprocessing': { 'image rotation': img_rot ,
'is shift correction': corr_shift} ,
'root tip': { 'dilation obstacle': dil_msk ,
'proximity connect': prox_conn_th ,
'proximity detect': prox_detec_th ,
'minimum root area': min_area_root ,
'dilate erode kernel': pict_smth ,
'smooth contour': cont_smth ,
'smooth curvature': curv_smth ,
'curvature threshold': curv_th ,
'smooth speed': spd_smth} ,
'main root': { 'smooth speed': main_spd_smth ,
'trajectory smooth': traj_smth ,
'mask smooth': msk_smth ,
'length smooth': lgth_smth ,
'error polynomial': err_pol ,
'distance direction': d_dir} ,
'displaying': { 'frame per second': fps_mov ,
'amplitude direction': amp_dir } ,
'contact': { 'contact detection front': cnt_th ,
'contact detection back': cnt_th_b ,
'radius contact': cnt_rad } ,
'computational': { 'number processor': nb_proc }
}
## Save it in a YAML file:
with io.open('input_parameters.yaml', 'w', encoding='utf8') as outfile:
yaml.dump(data_input, outfile, default_flow_style=False, allow_unicode=True)
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Intialisation:
#Load libraries:
##Load libraries:
import io
import yaml
import numpy as np
import math as m
from scipy import misc
from PIL import Image
from scipy import signal as sg
import fnmatch
import os
import matplotlib.pyplot as plt
from matplotlib.path import Path
import cv2
from cv2 import Canny
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
# Input parameter loading:
## Load dictionnary:
with io.open('input_parameters.yaml', 'r') as infile:
data_input = yaml.load(infile, Loader=yaml.FullLoader)
##Input parameter:
###Contact detection threshold [px]:
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
## Number of frame per second for movies:
fps_mov = data_input['displaying']['frame per second']
## Number of steps:
nb_stp = data_input['general']['number of step']
##Contact detection threshold:
cnt_th = data_input['contact']['contact detection front']
##Radius area considered to compute contact direction:
cnt_rad = data_input['contact']['radius contact']
##Contact detection threshold for root back:
cnt_th_b = data_input['contact']['contact detection back']
## Distance over wich to compute root tip direction:
d_dir = data_input['main root']['distance direction']
##Count the number of pictures in the folder to get the number of keeped steps:
NbStp=len(fnmatch.filter(os.listdir('pictureWh'),'*.jpg'))
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
# Get contour of already found obstacles:
## Load obstacle map:
map_obst = np.array(Image.open('obstacle/obstacleMask.png'))
### Loop around contour
obst_list = []
for it_obst in range(map_obst.max()):
#### Make mask:
mask_obst_cur = np.zeros(map_obst.shape)
mask_obst_cur[np.where(map_obst==it_obst+1)]=1
mask_obst_cur = mask_obst_cur.astype('uint8')
#### Get contour:
cont_cur,h = cv2.findContours(mask_obst_cur,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
cont_cur = cont_cur[0]
x_cur = cont_cur[:,0,0].astype('float')
y_cur = cont_cur[:,0,1].astype('float')
#### Close the contour:
x_cur = np.append(x_cur,x_cur[0])
y_cur = np.append(y_cur,y_cur[0])
#!Debug!
# ~ plt.imshow(mask_obst_cur)
# ~ plt.plot(x_cur,y_cur,'-o')
# ~ plt.show()
#### Store it:
obst_list.append(np.vstack((x_cur,y_cur)).T)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Get obstacles position on the picture:
# Get missing obstacles position on the picture:
##Load the picture on which obstacle are specified:
pict0=misc.imread('pictureWh/'+'%04d'%(NbStp-1)+'.jpg')
## Load the picture on which obstacle are specified:
pict0 = np.array(Image.open('picture_processed/'+'%04d'%(0)+'.jpg'))
##Select obstacles:
###Define a class to draw obstacles:
## Select obstacles:
### Define a class to draw obstacles:
class LineBuilder:
def __init__(self, line):
self.line=line
self.xs=list(line.get_xdata())
self.ys=list(line.get_ydata())
self.cid=line.figure.canvas.mpl_connect('button_press_event', self)
self.flag=True
self.xs = list(line.get_xdata())
self.ys = list(line.get_ydata())
self.cid = line.figure.canvas.mpl_connect('button_press_event', self)
self.flag = True
def __call__(self,event):
print('click',event)
if event.button==3 and event.inaxes:
if event.inaxes!=self.line.axes: return
if self.flag:
self.xs=[event.xdata]
self.ys=[event.ydata]
self.flag=False
self.xs = [event.xdata]
self.ys = [event.ydata]
self.flag = False
else:
self.xs.append(event.xdata)
self.ys.append(event.ydata)
self.line.set_data(self.xs, self.ys)
self.line.figure.canvas.draw()
###Loop over obstacles:
obst_list=[]
flag0=True
while flag0:
####Select an obstacle or close to finish:
fig=plt.figure()
mng=plt.get_current_fig_manager()
mng.resize(*mng.window.maxsize())
ax=fig.add_subplot(111)
### Loop over obstacles:
while True:
#### Select an obstacle or close to finish:
fig = plt.figure()
mng = plt.get_current_fig_manager()
ax = fig.add_subplot(111)
ax.imshow(pict0)
for it_l in range(len(obst_list)):
plt.plot(obst_list[it_l][:,0],obst_list[it_l][:,1],'-b',linewidth=1)
ax.set_title('Click around an obstacles clicks\n -right click to select, left click to zoom, close to end-')
line,=ax.plot([0],[0])
linebuilder=LineBuilder(line)
ax.set_title('Add missing obstacle (if any)\n Click around an obstacles clicks\n -right click to select, left click to zoom, close to end-')
line, = ax.plot([0],[0])
linebuilder = LineBuilder(line)
plt.show()
if len(linebuilder.xs)>1:
####Store the obstacle border:
x_cur=np.array(linebuilder.xs); x_cur=np.append(x_cur,x_cur[0])
y_cur=np.array(linebuilder.ys); y_cur=np.append(y_cur,y_cur[0])
#### Store the obstacle border:
x_cur = np.array(linebuilder.xs); x_cur = np.append(x_cur,x_cur[0])
y_cur = np.array(linebuilder.ys); y_cur = np.append(y_cur,y_cur[0])
obst_list.append(np.vstack((x_cur,y_cur)).T)
else:
####Close:
flag0=False
#obst_list.pop(6)
break
###Display obstacles:
####Make storage folder:
### Display obstacles:
#### Make storage folder:
os.system('mkdir obstacle')
####Create figure:
#### Create figure:
plt.imshow(pict0)
for it_l in range(len(obst_list)):
plt.plot(obst_list[it_l][:,0],obst_list[it_l][:,1],'-',linewidth=1)
......@@ -99,147 +132,140 @@ plt.axis('off')
plt.savefig('obstacle/numbering.png',bbox_inches='tight',dpi=200)
plt.close()
##Save obstacle coordinates:
mat_save=[]
## Save obstacle coordinates:
mat_save = []
for it_l in range(len(obst_list)):
for it_m in range(len(obst_list[it_l][:,0])):
mat_save.append(np.array([it_l+1,obst_list[it_l][it_m,0],obst_list[it_l][it_m,0]]))
np.savetxt('obstacle/coordinate.txt',np.vstack(mat_save))
##Make obstacle mask:
###Initialise mask:
obst_mask=np.zeros((pict0.shape[0],pict0.shape[1]))
I0,J0=np.meshgrid(np.arange(obst_mask.shape[0]),np.arange(obst_mask.shape[1]))
I0=I0.flatten(); J0=J0.flatten()
pict_point=np.vstack((J0,I0)).T
###Loop over obstacles:
## Make obstacle mask:
### Initialise mask:
obst_mask = np.zeros((pict0.shape[0],pict0.shape[1]))
I0,J0 = np.meshgrid(np.arange(obst_mask.shape[0]),np.arange(obst_mask.shape[1]))
I0 = I0.flatten(); J0 = J0.flatten()
pict_point = np.vstack((J0,I0)).T
### Loop over obstacles:
for it_ob in range(len(obst_list)):
####Load the current obstacle
obst_c=obst_list[it_ob]
####Make the path:
path_c=Path(obst_c)
####Get points inside:
grid_c=path_c.contains_points(pict_point)
####Update the mask:
mask_c=grid_c.reshape(pict0.shape[1],pict0.shape[0]).T
obst_mask+=(it_ob+1)*mask_c
###Save the mask:
misc.imsave('obstacle/bulkMask.png',obst_mask)
##Get obstacle edges:
###Get contour:
edge_raw=Canny(obst_mask.astype('uint8'),1,1)
###Number edges:
edge_raw=sg.convolve2d(edge_raw,np.ones((2,2)),mode='same')>100
obst_edge=edge_raw*obst_mask
###Save it:
misc.imsave('obstacle/edgeMask.png',obst_edge)
#### Load the current obstacle
obst_c = obst_list[it_ob]
#### Make the path:
path_c = Path(obst_c)
#### Get points inside:
grid_c = path_c.contains_points(pict_point)
#### Update the mask:
mask_c = grid_c.reshape(pict0.shape[1],pict0.shape[0]).T
obst_mask += (it_ob+1)*mask_c
### Save the mask:
im = Image.fromarray(obst_mask).convert('L')
im.save('obstacle/bulkMask.png')
## Get obstacle edges:
### Get contour:
edge_raw = Canny(obst_mask.astype('uint8'),1,1)
### Number edges:
edge_raw = sg.convolve2d(edge_raw,np.ones((2,2)),mode='same')>100
obst_edge = edge_raw*obst_mask
### Save it:
im = Image.fromarray(obst_edge).convert('L')
im.save('obstacle/edgeMask.png')
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Detect contact between main root tip and obstacles:
##Load obstacle:
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:
pos_tip=np.loadtxt('mainRoot/rootTip/position.txt')
time_tip=np.loadtxt('mainRoot/rootTip/time.txt')
##Load the time and root tip direction:
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)
##Get contact, obstacle number and contact direction:
###Initialise storage vectors:
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:
# Detect contact between main root tip and obstacles:
## Load obstacle:
obst_edge = np.array(Image.open('obstacle/edgeMask.png'))
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_smoothed.txt')
time_tip = np.loadtxt('mainRoot/rootTip/time_smoothed.txt')
## Load the time and root tip direction:
dir0 = np.loadtxt('mainRoot/direction.txt')
time0 = np.loadtxt('time[min].txt')
## Interpolate direction:
dir_tip = np.interp(time_tip,time0,dir0)
## Get contact, obstacle number and contact direction:
### Initialise storage vectors:
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)
N_o = obst_edge[I_o,J_o]
### Loop over positions:
for it_t in range(len(time_tip)):
####Get current position:
Jc,Ic=pos_tip[it_t,0],pos_tip[it_t,1],
####Get distance to obstacle edges:
dc=np.sqrt((Ic-I_o)**2.+(Jc-J_o)**2.)
####Contact if distance is smaller than thresold:
#### Get current position:
Jc,Ic = pos_tip[it_t,0],pos_tip[it_t,1],
#### Get distance to obstacle edges:
dc = np.sqrt((Ic-I_o)**2.+(Jc-J_o)**2.)
#### Contact if distance is smaller than thresold:
if np.min(dc)<cnt_th:
####Get obstacle number:
II=np.where(dc==np.min(dc))[0][0]
cnt_num[it_t]=N_o[II]
####Get contact direction:
#####Measure obstacle direction:
II=np.where(dc<rad_px)[0]
I_d=I_o[II]; J_d=J_o[II]
Do=-m.atan(np.polyfit(I_d,J_d,1)[0])+m.pi/2.
#####Load root tip direction:
Dc=-dir_tip[it_t]
#####Compute and store angle difference:
diffD=abs(m.atan2(m.sin(Dc-Do),m.cos(Dc-Do)))
cnt_dir[it_t]=diffD
#### Get obstacle number:
II = np.where(dc==np.min(dc))[0][0]
cnt_num[it_t] = N_o[II]
#### Get contact direction:
##### Measure obstacle direction:
II = np.where(dc<cnt_rad)[0]
I_d = I_o[II]; J_d = J_o[II]
Do = -m.atan(np.polyfit(I_d,J_d,1)[0])+m.pi/2.
##### Load root tip direction:
Dc = -dir_tip[it_t]
##### Compute and store angle difference:
diffD = abs(m.atan2(m.sin(Dc-Do),m.cos(Dc-Do)))
cnt_dir[it_t] = diffD
# !Debug!
#~ plt.plot(J_o,I_o,'ok')
#~ plt.plot(J_d,I_d,'xg')
#~ plt.plot(pos_tip[0:it_t,0],pos_tip[0:it_t,1],'-')
#~ plt.plot(Jc,Ic,'+r')
#~ nn=10.
#~ plt.plot(np.array([Jc,Jc+nn*m.cos(Dc)]),np.array([Ic,Ic+nn*m.sin(Dc)]),'-r')
#~ plt.plot(np.array([Jc,Jc+nn*m.cos(Do)]),np.array([Ic,Ic+nn*m.sin(Do)]),'-g')
#~ plt.axis('equal')
#~ plt.xlim(Jc-50,Jc+50)
#~ plt.ylim(Ic-50,Ic+50)
#~ plt.title(str(int(diffD/m.pi*180.)))
#~ plt.savefig('tmp/'+str(it_t).zfill(5)+'.png')
#~ plt.close()
###Save contact information:
# ~ plt.plot(J_o,I_o,'ok')
# ~ plt.plot(J_d,I_d,'xg')
# ~ plt.plot(pos_tip[0:it_t,0],pos_tip[0:it_t,1],'-m')
# ~ plt.plot(pos_tip[max(0,it_t-10):it_t,0],pos_tip[max(0,it_t-10):it_t,1],'.k')
# ~ plt.plot(Jc,Ic,'+r')
# ~ nn=10.
# ~ plt.plot(np.array([Jc,Jc+nn*m.cos(Dc)]),np.array([Ic,Ic+nn*m.sin(Dc)]),'-r')
# ~ plt.plot(np.array([Jc,Jc+nn*m.cos(Do)]),np.array([Ic,Ic+nn*m.sin(Do)]),'-g')
# ~ plt.axis('equal')
# ~ plt.xlim(Jc-50,Jc+50)
# ~ plt.ylim(Ic-50,Ic+50)
# ~ plt.title(str(int(diffD/m.pi*180.)))
# ~ 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:
# 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 obstacle mask:
obst_edge = np.array(Image.open('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 time and position of the tip:
time0 = np.loadtxt('time[min].txt')
##Load image edges:
IJ_edge=np.loadtxt('pictureSkeleton/IJ_edge.txt')
##Make storage folder:
## 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:
## Get contact and obstacle number:
### Prepare obstacle:
I_o,J_o = np.where(obst_edge>0)
N_o = obst_edge[I_o,J_o]
### 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]
#### 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 = np.array(Image.open('pictureSkeleton/weightedSkeleton/'+str(it_t).zfill(4)+'.png'))
skl_cur = skl_cur/np.sort(np.unique(skl_cur))[1]
##~DEBUG~##
# plt.imshow(skl_cur)
......@@ -247,12 +273,12 @@ for it_t in range(len(time0)):
# plt.show()
# plt.close()
####Loop along path to detect contact:
cnt_num=np.zeros(len(Ic))
#### 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:
##### 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)
......@@ -263,73 +289,71 @@ for it_t in range(len(time0)):
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]
cnt_num[it] = N_o[II]
####Save the result:
#### Save the result:
np.savetxt('mainRoot/rootContact/'+str(it_t).zfill(4)+'.txt',cnt_num)
#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#~/"\~#
#Display results:
# Display results:
##Make temporary storage:
## 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 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')
## Load root tip directon:
dir0 = np.loadtxt('mainRoot/direction.txt')
##Loop over time step:
## 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: