diff --git a/README.md b/README.md
index 04558db..ef9463a 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,5 @@
 ## PhenoTrack3D: an automatic high-throughput phenotyping pipeline to track maize organs over time
 
-This work is based on our biorxiv report (TODO), by Benoit Daviet, Romain Fernandez, Llorenc Cabrera-Bosquet, Christophe Pradal and Christian Fournier.
 PhenoTrack3D can be used to track leaf organs over time, in a time-series of 3D maize plant reconstructions. 
 Such 3D reconstructions can be obtained from sets of RGB images with the Phenomenal pipeline (https://github.com/openalea/phenomenal)
 
@@ -15,6 +14,19 @@ Such 3D reconstructions can be obtained from sets of RGB images with the Phenome
 Our code is released under **Cecill-C** (https://cecill.info/licences/Licence_CeCILL_V1.1-US.txt) licence. 
 See LICENSE file for details.
 
+### Install
+
+install dependencies with conda:
+
+	conda create -n phenotrack -c conda-forge -c openalea3 openalea.phenomenal skan=0.9
+	conda activate phenotrack
+
+Clone repo and run setup
+
+    git clone https://github.com/openalea/phenotrack3d.git
+    cd phenotrack3d
+
+
 ### References
 
 If you use PhenoTrack3D to your research, cite:
diff --git a/paper/azimuth.py b/paper/azimuth.py
new file mode 100644
index 0000000..00ef331
--- /dev/null
+++ b/paper/azimuth.py
@@ -0,0 +1,52 @@
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from openalea.maizetrack.phenomenal_display import PALETTE
+import os
+
+plantid = 940
+
+# result of tracking (without 3D object) saved by trackedplant method (copied to cachezA17/local_benoit)
+dfid = pd.read_csv('data/tracking/{}.csv'.format(plantid))
+#dfid = dfid[~dfid['mature']]
+#from llorenc, currently copied to  modulor cache
+df_tt = pd.read_csv('data/TT_ZA17.csv')
+
+thermaltimes = np.array(df_tt['ThermalTime'])
+timestamps = np.array(df_tt['timestamp'])
+dfid['tt'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in dfid['timestamp']]
+
+fig, ax = plt.subplots(figsize =(10, 10), dpi=100, subplot_kw={'projection': 'polar'})
+fig.canvas.set_window_title(str(plantid))
+
+ranks = sorted([r for r in dfid['rank_tracking'].unique() if r != 0])
+for r in ranks[:9]:
+    dfr = dfid[dfid['rank_tracking'] == r].sort_values('tt').iloc[:25]
+    dfr = dfr[dfr['tt'] - min(dfr['tt']) < 40]
+    time = dfr['tt'] - min(dfr['tt'])
+    azimuth = dfr['a'] / 360 * 2*np.pi
+    ax.plot(azimuth, time, '-', color=PALETTE[r - 1] / 255.)
+
+ax.set_rlabel_position(-22.5)  # Move radial labels away from plotted line
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.legend()
+
+fig.savefig('paper/azimuth_growth', bbox_inches='tight')
+
+
+fig, ax = plt.subplots(figsize =(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+fig.canvas.set_window_title(str(plantid))
+dfid2 = dfid[dfid['rank_tracking'] != 0]
+profile = []
+ranks = dfid2['rank_tracking'].unique()
+for r in ranks:
+    dfr = dfid2[dfid2['rank_tracking'] == r].sort_values('tt').iloc[:15]
+    profile.append(np.median(dfr['a']))
+profile[-2] -= 180
+plt.plot(ranks, profile, 'k*-', markersize=10)
+plt.xlabel('Rank', fontsize=30)
+plt.ylabel('Azimuth (°)', fontsize=30)
+
+fig.savefig('paper/azimuth_profile', bbox_inches='tight')
+
diff --git a/paper/collar_detection_example.py b/paper/collar_detection_example.py
new file mode 100644
index 0000000..a622586
--- /dev/null
+++ b/paper/collar_detection_example.py
@@ -0,0 +1,111 @@
+import numpy as np
+import cv2
+import matplotlib.pyplot as plt
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths
+from openalea.maizetrack.utils import phm3d_to_px2d, simplify, get_rgb
+from openalea.maizetrack.phenomenal_display import plot_sk
+
+from openalea.deepcollar.predict_insertion_DL import detect_insertion
+
+import openalea.phenomenal.object as phm_obj
+import openalea.phenomenal.segmentation as phm_seg
+
+from skimage import io
+
+has_pgl_display = True
+try:
+    from openalea.plantgl import all as pgl
+except ImportError:
+    has_pgl_display = False
+
+
+def plot_skeleton2(vmsi, stem_path):
+
+    size = 15
+    col = pgl.Material(pgl.Color3(255, 0, 0))
+    col2 = pgl.Material(pgl.Color3(213, 100, 0))
+
+    shapes = []
+    h = 700  # - vmsi.get_stem().info['pm_z_base']
+    for leaf in vmsi.get_leafs():
+        segment = leaf.get_highest_polyline().polyline # for leaf
+
+        segment = segment[:-2]
+
+        for k in range(len(segment) - 1):
+            # arguments cylindre : centre1, centre2, rayon, nbre segments d'un cercle.
+            pos1 = np.array(segment[k]) + np.array([0, 0, h])
+            pos2 = np.array(segment[k + 1]) + np.array([0, 0, h])
+            cyl = pgl.Extrusion(pgl.Polyline([pos1, pos2]), pgl.Polyline2D.Circle(int(size * 0.8), 8))
+            cyl.solid = True  # rajoute 2 cercles aux extremites du cylindre
+            cyl = pgl.Shape(cyl, col)
+            shapes.append(cyl)
+
+    segment = stem_path.polyline
+    for k in range(len(segment) - 1):
+        # arguments cylindre : centre1, centre2, rayon, nbre segments d'un cercle.
+        pos1 = np.array(segment[k]) + np.array([0, 0, h])
+        pos2 = np.array(segment[k + 1]) + np.array([0, 0, h])
+        cyl = pgl.Extrusion(pgl.Polyline([pos1, pos2]), pgl.Polyline2D.Circle(int(size * 1.2), 8))
+        cyl.solid = True  # rajoute 2 cercles aux extremites du cylindre
+        cyl = pgl.Shape(cyl, col2)
+        shapes.append(cyl)
+
+    scene = pgl.Scene(shapes)
+    pgl.Viewer.display(scene)
+
+
+plantid = 1424
+metainfos = get_metainfos_ZA17(plantid)
+
+daydate = '2017-05-12'
+metainfo = next(m for m in metainfos if m.daydate == daydate)
+
+skeleton_path = metainfos_to_paths([metainfo], phm_parameters=(4, 1, 'notop', 4, 100), object='skeleton')[0]
+vmsi_path = metainfos_to_paths([metainfo], phm_parameters=(4, 1, 'notop', 4, 100), object='vmsi')[0]
+
+weights = 'deepcollar/examples/data/model/weights.weights'
+config = 'deepcollar/examples/data/model/config.cfg'
+net = cv2.dnn.readNetFromDarknet(config, weights)
+model = cv2.dnn_DetectionModel(net)
+
+angles = [a for a, v in zip(metainfo.camera_angle, metainfo.view_type) if v == 'side']
+images = {angle: get_rgb(metainfo, angle, main_folder='data/rgb_insertion_annotation/rgb',# used only to retrieve a copy of rgb image
+                         plant_folder=False, save=False, side=True)[0] for angle in angles}
+
+# =============================================================================================================
+
+skeleton = phm_obj.VoxelSkeleton.read_from_json_gz(skeleton_path)
+
+pl = phm_seg.get_highest_segment(skeleton.segments).polyline
+pl3d = simplify(pl, 30)
+
+stem_polylines = {angle: phm3d_to_px2d(pl3d, metainfo.shooting_frame, angle=angle) for angle in angles}
+
+for angle in [60]:
+
+    # predict x,y
+    res = detect_insertion(image=images[angle], stem_polyline=stem_polylines[angle], model=model, display=True,
+                           display_vignettes=True)
+
+    img = images[angle]
+    pl = stem_polylines[angle]
+    img = cv2.polylines(np.float32(img), [pl.astype(int).reshape((-1, 1, 2))], False, (213, 120, 0), 12)
+    plt.imshow(img.astype(np.uint8))
+
+    io.imsave('stem_path.png', img)
+
+stem_segment = phm_seg.get_highest_segment(skeleton.segments)
+vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(vmsi_path)
+plot_skeleton2(vmsi, stem_segment)
+
+
+# import os
+# folder = 'paper/method/'
+# for path in os.listdir(folder):
+#     img = io.imread(folder + path)
+#     img = img[:, :, :3]
+#     img[(img == np.array([255, 255, 255])).all(2)] = [120, 120, 120]
+#     img[(img == np.array([0, 0, 0])).all(2)] = [255, 255, 255]
+#     io.imsave(folder + 'v2' + path, img)
diff --git a/paper/dt_sensi.py b/paper/dt_sensi.py
new file mode 100644
index 0000000..b6082b0
--- /dev/null
+++ b/paper/dt_sensi.py
@@ -0,0 +1,92 @@
+import os
+import copy
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.trackedPlant import TrackedPlant
+
+folder = 'data/tracking/'
+plantids = [int(f.split('.')[0]) for f in os.listdir(folder) if os.path.isfile(folder + f)]
+
+df_list = []
+for plantid in plantids:
+    dfid = pd.read_csv(folder + '{}.csv'.format(plantid))
+    df_list.append(dfid)
+df = pd.concat(df_list)
+df = df[~df['rank_annotation'].isin([-1, 0])]
+# df = df[df['timestamp'] < 1496509351]
+# df = df[df['timestamp'] < 1495922400] # 05-28
+
+# ====== mature tracking, but with different end dates ============================================================*
+
+# selec = df
+n, accuracy = [], []
+dt = {}
+
+for dt_div in [1, 2, 4]:
+    dt[dt_div] = []
+
+    for plantid in df['plantid'].unique():
+
+        print(plantid)
+
+        metainfos = get_metainfos_ZA17(plantid)
+
+        metainfos = sorted(metainfos, key=lambda k: k.timestamp)
+
+        paths = metainfos_to_paths(metainfos, folder='local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking')
+        metainfos, paths = check_existence(metainfos, paths)
+
+        dt[dt_div].append(np.diff([m.timestamp / 3600 for m in metainfos]))
+
+        metainfos, paths = metainfos[::dt_div], paths[::dt_div]
+
+        vmsi_list = load_plant(metainfos, paths)
+
+        plant_ref = TrackedPlant.load_and_check(vmsi_list)
+        plant_ref.load_rank_annotation()
+
+        dmax = '2017-06-03'
+
+        plant = copy.deepcopy(plant_ref)
+
+        plant.snapshots = [s for s in plant.snapshots if s.metainfo.daydate < dmax]
+
+        plant.align_mature(direction=1, gap=12.365, w_h=0.03, w_l=0.002, gap_extremity_factor=0.2, n_previous=500,
+                           rank_attribution=True)
+        plant.align_growing()
+
+        df_res = plant.get_dataframe(load_anot=False)
+        print(dt_div, round(len(df_res[df_res['rank_tracking'] == df_res['rank_annotation']]) / len(df_res), 2))
+        df_res.to_csv('data/tracking/test_dt/{}_{}.csv'.format(plantid, dt_div))
+
+plt.figure()
+plt.hist(np.concatenate(dt[2]), 80)
+plt.xlabel('Δt between two consecutive images (h)')
+
+for dt_div in [1, 2, 4]:
+    acc1, acc2, n = [], [], []
+    for plantid in plantids:
+
+        df = pd.read_csv('data/tracking/test_dt/{}_{}.csv'.format(plantid, dt_div))
+
+        selecid = df[~df['rank_annotation'].isin([-1, 0])]
+        selec = selecid[(selecid['mature'] == True)]
+        a1 = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec)
+        selec = selecid[(selecid['mature'] == False)]
+        a2 = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec)
+        acc1.append(a1)
+        acc2.append(a2)
+        # print(plantid, round(a1, 3), round(a2, 3), len(selecid))
+        n.append(len(selec))
+    print('===== {} =========================='.format(dt_div))
+    print(min(acc1), np.mean(acc1))
+    print(min(acc2), np.mean(acc2))
+    print(min(n), np.mean(n))
+
+
+
+
+
diff --git a/paper/figure_poster.py b/paper/figure_poster.py
new file mode 100644
index 0000000..5b4da36
--- /dev/null
+++ b/paper/figure_poster.py
@@ -0,0 +1,82 @@
+import os
+import numpy as np
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.phenomenal_display import PALETTE, plot_vmsi_voxel
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.plantgl import all as pgl
+
+
+folder = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+plantids = [313, 316, 329, 330, 336, 348, 424, 439, 461, 474, 794, 832, 905, 907, 915, 925, 931, 940, 959, 1270,
+            1276, 1283, 1284, 1301, 1316, 1383, 1391, 1421, 1424, 1434]
+
+plantid = 348 # 959  # 1383
+
+metainfos = get_metainfos_ZA17(plantid)
+paths = metainfos_to_paths(metainfos, folder=folder)
+metainfos, paths = check_existence(metainfos, paths)
+vmsi_list = load_plant(metainfos, paths)
+
+plant = TrackedPlant.load_and_check(vmsi_list)
+plant.load_rank_annotation()
+
+# [2, 12, 22, 32, 45], cam = (0, 41, manuel bloqué, 112)
+
+phm_rank = True
+
+for i in [2, 12, 22, 32, 45]:
+
+    snapshot = plant.snapshots[i]
+
+    size = snapshot.voxels_size
+    shapes = []
+
+    organs = [[(0, 0, 0), snapshot.get_stem()]]
+    colors = np.random.choice(np.arange(18), len(snapshot.leaves), replace=False)
+    for k, (leaf, rank, i_col) in enumerate(zip(snapshot.leaves, snapshot.rank_annotation, colors)):
+        #rank_displayed = k if phm_rank else rank
+        #organs.append([tuple(PALETTE[rank_displayed]), leaf])
+        organs.append([tuple(PALETTE[i_col]), leaf])
+
+    ranks = None
+
+    for col, organ in organs:
+
+        c1, c2, c3 = col
+
+        m = pgl.Material(pgl.Color3(int(c1), int(c2), int(c3)))
+        for x, y, z in list(organ.voxels_position())[::1]:
+            m = pgl.Material(pgl.Color3(int(c1), int(c2), int(c3)))
+            b = pgl.Box(size, size, size)
+            vx = pgl.Translated(x, y, z, b)
+            vx = pgl.Shape(vx, m)
+            shapes.append(vx)
+
+    scene = pgl.Scene(shapes)
+
+    pgl.Viewer.display(scene)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/paper/growing_leaf_length.py b/paper/growing_leaf_length.py
new file mode 100644
index 0000000..4ecb93b
--- /dev/null
+++ b/paper/growing_leaf_length.py
@@ -0,0 +1,108 @@
+"""
+Test accuracy of leaf length of growing leaves after phenomenal + tracking
+"""
+
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import os
+from sklearn.metrics import r2_score
+from matplotlib.ticker import MaxNLocator
+
+# available on  NAsshare2
+df_anot = pd.read_csv('data/rgb_leaf_growing_annotation/leaf_growing_annotation.csv')
+
+plantids = df_anot['plantid'].unique()
+
+res = []
+n = 0
+for plantid in plantids:
+    # available on modulor local_benoit
+    df_pred = pd.read_csv('data/tracking/{}.csv'.format(plantid))
+    df_obs = df_anot[df_anot['plantid'] == plantid]
+
+    for _, row_obs in df_obs.iterrows():
+        selec = df_pred[(df_pred['rank_tracking'] == row_obs['rank']) & (df_pred['timestamp'] == row_obs['timestamp'])]
+        if not selec.empty:
+            n += ((selec['rank_tracking'] == selec['rank_annotation']).iloc[0])
+            row_pred = selec.iloc[0]
+            res.append([plantid, row_obs['rank'], row_obs['length_mm'], row_pred['l'], row_pred['l_extended']])
+
+res = pd.DataFrame(res, columns=['plantid', 'rank', 'obs', 'pred_phm', 'pred_ext'])
+
+
+for rank in [6, 9]:
+    selec = res[res['rank'] == rank]
+    fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+    ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+    plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+    plt.title('Length of leaf {} \nduring growing phase'.format(rank), fontsize=35)
+    plt.xlabel('observation (cm)', fontsize=30)
+    plt.ylabel('prediction (cm)', fontsize=30)
+    plt.plot([-20, 160], [-20, 160], '-', color='k')
+    plt.xlim((0, 130))
+    plt.ylim((0, 130))
+    x, y = selec['obs'] / 10, selec['pred_ext'] / 10
+    plt.plot(x, y, 'ro', markersize=8, alpha=0.5)
+    # x2, y2 = selec['obs'] / 10, selec['pred_phm'] / 10
+    # plt.plot(x2, y2, '^', color='grey', alpha=0.7, markersize=15)
+
+    a, b = np.polyfit(x, y, 1)
+    plt.plot([-10, 3000], a * np.array([-10, 3000]) + b, 'r--', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+    leg = plt.legend(prop={'size': 25}, loc='upper left', markerscale=2)
+    leg.get_frame().set_linewidth(0.0)
+
+    rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+    r2 = r2_score(x, y)
+    biais = np.mean(y - x)
+    mape = 100 * np.mean(np.abs((x - y) / x))
+    ax.text(0.60, 0.26, 'n = {} \nBias = {} cm \nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x), round(biais, 2), round(rmse, 2), round(mape, 2), round(r2, 3)),
+            transform=ax.transAxes, fontsize=25, verticalalignment='top')
+
+    fig.savefig('paper/results/growing_leaf_{}'.format(rank), bbox_inches='tight')
+
+# .format(len(x), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3))
+
+# ===== trait visu for one plant =================================================================================
+
+plantid = 931
+
+df_pred = pd.read_csv('data/tracking/{}.csv'.format(plantid))
+df_obs = df_anot[df_anot['plantid'] == plantid]
+
+df_tt = pd.read_csv('TT_ZA17.csv')
+thermaltimes = np.array(df_tt['ThermalTime'])
+timestamps = np.array(df_tt['timestamp'])
+df_obs['tt'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in df_obs['timestamp']]
+df_pred['tt'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in df_pred['timestamp']]
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+plt.title('Leaf growth curve \nexamples for one plant', fontsize=35)
+plt.ylabel('Leaf length (cm)', fontsize=30)
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+
+for rank in [6, 9]:
+    s_obs = df_obs[df_obs['rank'] == rank]
+    symbol = 'o' if rank == 6 else '*'
+    plt.plot(s_obs['tt'], s_obs['length_mm'] / 10, 'g--' + symbol, markersize=8, label=rank)
+    s_pred = df_pred[(df_pred['rank_tracking'] == rank) & (df_pred['timestamp'] <= max(s_obs['timestamp']))]
+    plt.plot(s_pred['tt'], s_pred['l_extended'] / 10, 'k-' + symbol, markersize=8, label=rank)
+
+fig.savefig('paper/results/growth_example', bbox_inches='tight')
+
+
+# ===== corresponding legend ============================================
+
+plt.plot(0, 0, 'ko', label='Leaf 6')
+plt.plot(0, 0, 'k*', label='Leaf 9')
+plt.plot(0, 0, 'g--', label='observation')
+plt.plot(0, 0, 'k-', label='prediction')
+plt.legend()
+leg = plt.legend(prop={'size': 17}, markerscale=2)
+for line in leg.get_lines():
+    line.set_linewidth(4.0)
+leg.get_frame().set_linewidth(0.0)
+
diff --git a/paper/gxe.py b/paper/gxe.py
new file mode 100644
index 0000000..b872b28
--- /dev/null
+++ b/paper/gxe.py
@@ -0,0 +1,507 @@
+import pandas as pd
+import numpy as np
+import os
+from skimage import io
+import matplotlib.pyplot as plt
+from openalea.maizetrack.utils import get_rgb
+from openalea.maizetrack.local_cache import get_metainfos_ZA17
+from openalea.maizetrack.stem_correction import smoothing_function, stem_height_smoothing
+from scipy.interpolate import interp1d
+import time
+import datetime
+from matplotlib.ticker import MaxNLocator
+
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+
+def mad(x):
+    """ median absolute deviation """
+    return np.median(np.abs(x - np.median(x)))
+
+
+# thermal time : available in modulor
+df_tt = pd.read_csv('TT_ZA17.csv')
+
+# collar detection data
+# available in modulor
+folder_dl = 'local_cache/cache_ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+all_files_dl = [folder_dl + '/' + rep + '/' + f for rep in os.listdir(folder_dl) for f in os.listdir(folder_dl + '/' + rep)]
+
+# ZA17 GxE data
+# available in modulor
+pheno_ZA17 = pd.read_csv('pheno_ZA17.csv')
+pheno_ZA17['timestamp'] = pheno_ZA17['resultdate'].map(date_to_timestamp)
+gxe = pheno_ZA17[['Pot', 'Variety_ID', 'Scenario', 'Rep', 'Sampling']].drop_duplicates()
+gxe = gxe[gxe['Sampling'] == 'Rep']
+gxe = gxe[gxe['Variety_ID'] != 'SolNu']
+genotypes = ['DZ_PG_34', 'DZ_PG_41', 'DZ_PG_68', 'DZ_PG_69', 'DZ_PG_72']
+gxe = gxe[gxe['Variety_ID'].isin(genotypes)]
+gxe = gxe[gxe['Pot'] != 1070]  # not in cache
+
+
+# ============
+
+df = pheno_ZA17[pheno_ZA17['Variety_ID'] == 'DZ_PG_41']
+df = df[df['observationcode'].isin(['ligul_f', 'visi_f'])]
+thermaltimes = np.array(df_tt['ThermalTime'])
+timestamps = np.array(df_tt['timestamp'])
+df['tt'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in df['timestamp']]
+plt.figure()
+for var, symb in zip(['ligul_f', 'visi_f'], ['o', '^']):
+    for scenario, color, decalage in zip(['WW', 'WD'], ['blue', 'red'], [0, 0.3]):
+        selec = df[(df['Scenario'] == scenario) & (df['observationcode'] == var)].sort_values('tt')
+        plt.plot(selec['tt'] + decalage, selec['observation'], symb, color=color, markersize=4)
+
+        f, _ = smoothing_function(np.array(selec['tt']), np.array(selec['observation']), dy=None, dw=4, polyorder=2, repet=2)
+        plt.plot(selec['tt'], [f(t) for t in np.array(selec['tt'])], '-', color=color)
+
+
+
+colors = {'WW': 'dodgerblue', 'WD': 'orangered', 'error': 'gainsboro'}
+
+# # get rgb
+# for g in gxe['Variety_ID'].unique():
+#     #for e in gxe['Scenario'].unique():
+#     for e in ['WW']:
+#         selec = gxe[(gxe['Variety_ID'] == g) & (gxe['Scenario'] == e)]
+#         plantid = selec.iloc[0]['Pot']
+#         metainfos = get_metainfos_ZA17(plantid)
+#         metainfo = next(m for m in metainfos if m.daydate == '2017-05-20')
+#         img, _ = get_rgb(metainfo, 60, main_folder='rgb', plant_folder=True, save=False, side=True)
+#         io.imsave('paper/ZA17_visu/{}_{}_2017-05-20_angle60.png'.format(g, e), img)
+
+genotype = 'DZ_PG_41'
+
+genotypes = ['DZ_PG_41', 'DZ_PG_69']
+genotypes = ['DZ_PG_34', 'DZ_PG_41', 'DZ_PG_68', 'DZ_PG_69', 'DZ_PG_72']
+for genotype in genotypes:
+
+    # ===== stem height = f(t) ==================================================================================
+
+    plt.figure()
+    res_stem = []
+    for scenario in ['WW', 'WD']:
+
+        selec = gxe[(gxe['Variety_ID'] == genotype) & (gxe['Scenario'] == scenario)]
+        for _, row in selec.iterrows():
+
+            dfid = []
+            plantid_files = [f for f in all_files_dl if int(f.split('/')[-1][:4]) == row['Pot']]
+            for f in plantid_files:
+                task = int(f.split('_')[-1].split('.csv')[0])
+                #timestamp = next(m for m in metainfos if m.task == task).timestamp
+                daydate = f.split('__')[-2]
+                tt = df_tt[df_tt['daydate'] == daydate].iloc[0]['ThermalTime']
+                dft = pd.read_csv(f)
+                dft['tt'] = tt
+                dfid.append(dft)
+            dfid = pd.concat(dfid)
+
+            # remove > 06-08
+            dfid = dfid[dfid['tt'] < 90.5]
+
+            #  remove big time gaps
+            times = list(dfid.sort_values('tt')['tt'].unique())
+            for i in range(1, len(times)):
+                if (times[i] - times[i - 1]) > 5 * np.mean(np.diff(times)):
+                    dfid = dfid[dfid['tt'] < times[i]]
+
+            # data preprocessing, stem extraction
+            dfid = dfid[dfid['score'] > 0.95]
+            df_stem = dfid.sort_values('z_phenomenal', ascending=False).drop_duplicates(['tt']).sort_values('tt')
+            f_smoothing = stem_height_smoothing(np.array(df_stem['tt']), np.array(df_stem['z_phenomenal']))
+            df_stem['z_smooth'] = [f_smoothing(t) for t in df_stem['tt']]
+
+            for _, row_stem in df_stem.iterrows():
+                res_stem.append([row['Pot'], scenario, row_stem['z_smooth'], row_stem['tt']])
+
+            plt.plot(df_stem['tt'], df_stem['z_phenomenal'], '-', color=colors[scenario])
+            #plt.plot(df_stem['tt'], df_stem['z_smooth'], '-', color=colors[scenario])
+
+    res_stem = pd.DataFrame(res_stem, columns=['plantid', 'scenario', 'z', 't'])
+
+    res_gxe = res_stem.copy()
+    res_gxe['z'] /= 10
+    res_gxe = res_gxe[res_gxe['t'] < np.min(res_gxe.groupby(['scenario']).max('t')['t'])] # same tmax for each scenario
+    res_gxe = res_gxe.groupby(['scenario', 't']).agg(['median', ('mad', mad), 'count']).reset_index()
+    res_gxe = res_gxe[res_gxe['z']['count'] > 2] # at least 3 rep
+    res_gxe = res_gxe[res_gxe['t'] <= min(res_gxe.groupby('scenario').max('t')['t'])]  # same end for each scenario
+
+    fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+    fig.canvas.set_window_title(genotype)
+    ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+    for scenario in ['WW', 'WD']:
+        selec = res_gxe[res_gxe['scenario'] == scenario]
+        plt.fill_between(selec['t'], selec['z']['median'] - selec['z']['mad']/1, selec['z']['median'] + selec['z']['mad']/1,
+                         color=colors['error'])
+        plt.plot(selec['t'], selec['z']['median'], '.-', color=colors[scenario], markersize=12)
+        plt.plot(selec['t'], selec['z']['median'], '.', color=colors[scenario], label=scenario, markersize=12)
+    # lgd = plt.legend(prop={'size': 20}, loc='upper left', markerscale=2, handletextpad=-0.5)
+    # lgd.get_frame().set_linewidth(0.0)
+    plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+    plt.ylabel('Stem height (cm)', fontsize=30)
+    plt.xlim([8, 93])
+    plt.ylim([-2, 202])
+
+    fig.savefig('paper/stem_{}'.format(genotype), bbox_inches='tight')
+    plt.close('all')
+
+    # ===== length profile ==================================================================================
+
+    res_profile_l = []
+    res_profile_h = []
+    for scenario in ['WW', 'WD']:
+
+        selec = gxe[(gxe['Variety_ID'] == genotype) & (gxe['Scenario'] == scenario)]
+        for _, row in selec.iterrows():
+
+            # available in modulor local_benoit
+            file = 'data/tracking/gxe/{}.csv'.format(row['Pot'])
+            if not os.path.isfile(file):
+                print('no file', row['Pot'])
+            else:
+                dfid = pd.read_csv(file)
+
+                dfid_mature = dfid[dfid['mature']]
+                ranks = sorted([r for r in dfid_mature['rank_tracking'].unique() if r != 0])
+                profile = []
+                for r in ranks:
+                    dfr = dfid_mature[dfid_mature['rank_tracking'] == r]
+                    dfr = dfr.sort_values('timestamp')[:15]
+
+                    res_profile_l.append([row['Pot'], scenario, np.median(dfr['l_extended']), r])
+                    # TODO : height correction
+                    res_profile_h.append([row['Pot'], scenario, 720 + np.median(dfr['h']), r])
+
+    res_profile_l = pd.DataFrame(res_profile_l, columns=['plantid', 'scenario', 'l', 'r'])
+    res_profile_h = pd.DataFrame(res_profile_h, columns=['plantid', 'scenario', 'h', 'r'])
+
+    res_gxe = res_profile_l.copy()
+    res_gxe['l'] /= 10
+
+    plt.figure('all reps {}'.format(genotype))
+    for plantid in res_gxe['plantid'].unique():
+        selec = res_gxe[res_gxe['plantid'] == plantid]
+        plt.plot(selec['r'], selec['l'], '-', color=colors[selec.iloc[0]['scenario']])
+
+    res_gxe = res_gxe.groupby(['scenario', 'r']).agg(['median', ('mad', mad), 'count']).reset_index()
+    #res_gxe = res_gxe[res_gxe['l']['count'] > 2] # at least 3 rep
+    #res_gxe = res_gxe[res_gxe['r'] <= min(res_gxe.groupby('scenario').max('r')['r'])]  # same end for each scenario
+
+    fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+    fig.canvas.set_window_title(genotype + '_profile')
+    ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+    ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+    for scenario in ['WW', 'WD']:
+        selec = res_gxe[res_gxe['scenario'] == scenario]
+        plt.fill_between(selec['r'], selec['l']['median'] - selec['l']['mad']/1, selec['l']['median'] + selec['l']['mad']/1,
+                         color=colors['error'])
+        plt.plot(selec['r'], selec['l']['median'], '.-', color=colors[scenario], markersize=12)
+        plt.plot(selec['r'], selec['l']['median'], '.', color=colors[scenario], label=scenario, markersize=12)
+    # lgd = plt.legend(prop={'size': 20}, loc='upper left', markerscale=2, handletextpad=-0.5)
+    # lgd.get_frame().set_linewidth(0.0)
+    plt.xlabel('Rank', fontsize=30)
+    plt.ylabel('Leaf length (cm)', fontsize=30)
+    plt.xlim([0.3, 18.5])
+    plt.ylim([-3, 120])
+
+    fig.savefig('paper/profile_{}'.format(genotype), bbox_inches='tight')
+    plt.close('all')
+
+    # ===== leaf stage =================================================================================
+
+    plt.figure(genotype)
+
+    res_emergence = []
+    for scenario in ['WW', 'WD']:
+
+        selec = gxe[(gxe['Variety_ID'] == genotype) & (gxe['Scenario'] == scenario)]
+        for _, row in selec.iterrows():
+
+            # ===== emergence visi ====================================
+
+            # available in modulor local_benoit
+            file = 'data/tracking/gxe/{}.csv'.format(row['Pot'])
+            if not os.path.isfile(file):
+                print('no file', row['Pot'])
+            else:
+                dfid = pd.read_csv(file)
+
+                # only reps continuing after july
+                if len(set(dfid[dfid['timestamp'] > 1495806823]['timestamp'])) < 3:
+                    short_rep = True
+                else:
+                    short_rep = False
+
+                print(row['Pot'], scenario)
+
+                # # stop at panicle emergence
+                # df_phenoid = pheno_ZA17[pheno_ZA17['Pot'] == row['Pot']]
+                # select = df_phenoid[
+                #     (df_phenoid['observationcode'] == 'panicule_visi') & (df_phenoid['observation'] == 1)]
+                # t_max = max(dfid['timestamp'])
+                # if not select.empty:
+                #     t_max = min(select['timestamp'])
+                # timestamps = [t for t in dfid['timestamp'].unique() if t < t_max]
+
+                timestamps = [t for t in dfid['timestamp'].unique()]
+
+                # # TODO : useful with new gxe data ?
+                # # remove last day
+                # timestamps = sorted(timestamps)[:-1]
+
+                df_visi = []
+                for t in timestamps:
+
+                    dfidt = dfid[dfid['timestamp'] == t]
+
+                    if short_rep:
+                        mature_ranks = dfidt[dfidt['mature']]['rank_tracking']
+                        if mature_ranks.empty:
+                            n_mature = 0
+                        else:
+                            n_mature = max(mature_ranks)
+                        n_growing = len(dfidt[dfidt['mature'] == False])
+                        n_leaves = n_mature + n_growing
+
+                    else:
+                        n_leaves = max(dfidt['rank_tracking'])
+
+                    df_visi.append([t, n_leaves])
+                df_visi = pd.DataFrame(df_visi, columns=['t', 'n'])
+
+                median_visi = df_visi.groupby('n').median('t').reset_index()
+
+                # add missing n values (interpolation)
+                for n in range(min(median_visi['n']) + 1, max(median_visi['n'])):
+                    if n not in list(median_visi['n']):
+                        t1 = median_visi[median_visi['n'] < n].sort_values('n').iloc[-1]['t']
+                        t2 = median_visi[median_visi['n'] > n].sort_values('n').iloc[0]['t']
+                        median_visi = median_visi.append({'n': n, 't': (t1 + t2) / 2}, ignore_index=True)
+
+                # compute timing of leaf appearance (except for first and last n)
+                emerg_visi = []
+                for n in sorted(list(median_visi['n'].unique()))[1:-1]:
+                    t = np.mean(median_visi[median_visi['n'].isin([n - 1, n])]['t'])
+                    emerg_visi.append([row['Pot'], scenario, t, n, 'visi'])
+                emerg_visi = pd.DataFrame(emerg_visi, columns=['plantid', 'scenario', 't', 'n', 'type'])
+
+                # force monotony : t = f(n) can only increase
+                emerg_visi = emerg_visi.sort_values('n')
+                for i_row, rowe in emerg_visi.iterrows():
+                    previous = emerg_visi[emerg_visi['n'] < rowe['n']]
+                    if not previous.empty:
+                        emerg_visi.iloc[i_row, emerg_visi.columns.get_loc('t')] = max(rowe['t'], max(previous['t']))
+
+                thermaltimes = np.array(df_tt['ThermalTime'])
+                timestamps = np.array(df_tt['timestamp'])
+                emerg_visi['t'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in emerg_visi['t']]
+
+                plt.plot(emerg_visi['t'], emerg_visi['n'], '-*', label=row['Pot'])
+
+                res_emergence.append(emerg_visi)
+
+                # ===== emergence ligu =======================================
+
+                df_stem = res_stem[res_stem['plantid'] == row['Pot']]
+                df_profile = res_profile_h[res_profile_h['plantid'] == row['Pot']]
+
+                # debug : some profiles may have no data for a rank (??)
+                for r in range(min(df_profile['r']) + 1, max(df_profile['r'])):
+                    if r not in list(df_profile['r']):
+                        h1 = df_profile.sort_values('r')[df_profile['r'] < r].iloc[-1]['h']
+                        h2 = df_profile.sort_values('r')[df_profile['r'] > r].iloc[0]['h']
+                        h = (h1 + h2) / 2
+                        df_profile = df_profile.append({'plantid':df_profile.iloc[0]['plantid'], 'scenario': scenario, 'h': h, 'r': r},
+                                                       ignore_index=True).sort_values('r')
+
+                T = np.linspace(min(df_stem['t']), max(df_stem['t']), 3000)
+                f = interp1d(df_stem['t'], df_stem['z'])
+                Z = np.array([f(ti) for ti in T])
+
+                emerg_ligu = []
+                ranks = sorted(df_profile['r'])
+                for r in ranks[1:]:
+                    t = T[np.argmin(np.abs(Z - df_profile[df_profile['r'] == r - 1].iloc[0]['h']))]
+                    emerg_ligu.append([row['Pot'], scenario, t, r, 'ligu'])
+                emerg_ligu = pd.DataFrame(emerg_ligu, columns=['plantid', 'scenario', 't', 'n', 'type'])
+
+                res_emergence.append(emerg_ligu)
+
+    res_emergence = pd.concat(res_emergence)
+
+    # ===== plot =======
+
+    linestyles = {'visi': 'solid', 'ligu': 'dashed'}
+
+    # brut
+    plt.figure('{} brut'.format(genotype))
+    for var in ['visi', 'ligu']:
+        for plantid in res_emergence['plantid'].unique():
+            selec = res_emergence[(res_emergence['type'] == var) & (res_emergence['plantid'] == plantid)]
+            plt.plot(selec['t'], selec['n'], '.', linestyle=linestyles[var], color=colors[selec.iloc[0]['scenario']])
+
+    # interpolation
+    res_stage = []
+    T = np.linspace(min(res_emergence['t']), max(res_emergence['t']), 20)
+    for var in ['visi', 'ligu']:
+        for plantid in res_emergence['plantid'].unique():
+            selec = res_emergence[(res_emergence['type'] == var) & (res_emergence['plantid'] == plantid)]
+
+            if var == 'visi':
+                f = interp1d(selec['t'], selec['n'])
+            elif var == 'ligu':
+                f = interp1d(selec['t'], selec['n'])
+
+            t_interpolation = [t for t in T if min(selec['t']) < t < max(selec['t'])]
+            for t in t_interpolation:
+                res_stage.append([plantid, selec.iloc[0]['scenario'], t, float(f(t)), var])
+
+    res_stage = pd.DataFrame(res_stage, columns=['plantid', 'scenario', 't', 'n', 'type'])
+
+    res_gxe = res_stage.groupby(['t', 'scenario', 'type']).agg(['median', 'count', ('mad', mad)]).reset_index()
+    #res_gxe = res_gxe[res_gxe['n']['count'] > 2]  # at least 3 rep
+
+    # tmax_ligu = min(res_gxe[res_gxe['type'] == 'ligu'].groupby('scenario').max('t')['t'])
+    # tmax_visi = min(res_gxe[res_gxe['type'] == 'visi'].groupby('scenario').max('t')['t'])
+    # # same end
+    # res_gxe = res_gxe[(res_gxe['type'] == 'visi') & (res_gxe['t'] <= tmax_visi) |
+    #                   (res_gxe['type'] == 'ligu') & (res_gxe['t'] <= tmax_ligu)]
+
+    fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+    fig.canvas.set_window_title(genotype + '_stage')
+    ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+    for var in ['visi', 'ligu']:
+        for scenario in res_gxe['scenario'].unique():
+            selec = res_gxe[(res_gxe['scenario'] == scenario) & (res_gxe['type'] == var)]
+            plt.plot(selec['t'], selec['n']['median'], '.', linestyle=linestyles[var], color=colors[scenario])
+            plt.fill_between(selec['t'],
+                             selec['n']['median'] - selec['n']['mad']/1,
+                             selec['n']['median'] + selec['n']['mad']/1,
+                             color=colors['error'])
+
+    handles = []
+    for label, linestyle in zip(['visible', 'ligulated'], ['solid', 'dashed']):
+        h, = plt.plot(-100, -100, linestyle=linestyle, markersize=12, label=label, color='black')
+        handles.append(h)
+    lgd1 = plt.legend(prop={'size': 20}, handles=handles, loc='lower right', markerscale=2)
+    lgd1.get_frame().set_linewidth(0.0)
+    plt.gca().add_artist(lgd1)
+    handles = []
+    for scenario in ['WW', 'WD']:
+        h, = plt.plot(-100, -100, '.', markersize=12, label=scenario, color=colors[scenario])
+        handles.append(h)
+    # lgd2 = plt.legend(prop={'size': 20}, handles=handles, loc='upper left', markerscale=2, handletextpad=-0.5)
+    # lgd2.get_frame().set_linewidth(0.0)
+
+    #xlim, ylim = ax.get_xlim(), ax.get_ylim()
+    plt.xlim((7, 93))
+    plt.ylim((2, 18))
+    plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+    plt.ylabel('Leaf stage', fontsize=30)
+
+    fig.savefig('paper/stage_{}'.format(genotype), bbox_inches='tight')
+    plt.close('all')
+
+
+    # ===== leaf growth =================================================================================
+    for RANK in [6, 9]:
+
+        plt.figure('{}_{}'.format(genotype, RANK))
+        res_growth = []
+        for scenario in ['WW', 'WD']:
+
+            selec = gxe[(gxe['Variety_ID'] == genotype) & (gxe['Scenario'] == scenario)]
+            for _, row in selec.iterrows():
+
+                # h = length_profile(rank = 8)
+                selec2 = res_profile_l[(res_profile_l['plantid'] == row['Pot']) & (res_profile_l['r'] == RANK)]
+                if selec2.empty:
+                    print('plantid {} no ligulated height'.format(row['Pot']))
+                else:
+                    l_max = selec2.iloc[0]['l']
+
+                    # available in modulor local_benoit
+                    file = 'data/tracking/gxe/{}.csv'.format(row['Pot'])
+                    if not os.path.isfile(file):
+                        print('plantid {} no tracking file'.format(row['Pot']))
+                    else:
+                        dfid = pd.read_csv(file)
+                        dfr = dfid[dfid['rank_tracking'] == RANK].sort_values('timestamp').iloc[:20]
+
+                        thermaltimes = np.array(df_tt['ThermalTime'])
+                        timestamps = np.array(df_tt['timestamp'])
+                        dfr['tt'] = [thermaltimes[np.argmin(np.abs(timestamps - t))] for t in dfr['timestamp']]
+
+                        # premiere date ou on atteint 70% longueur max
+                        t_start = dfr[dfr['l_extended'] > 0.70 * l_max].sort_values('timestamp').iloc[0]['tt']
+
+                        # method 2
+                        selec3 = res_emergence[(res_emergence['plantid'] == row['Pot']) &
+                                      (res_emergence['type'] == 'visi') &
+                                      (res_emergence['n'] == RANK)]
+                        if selec3.empty:
+                            print('no emergence date')
+                        else:
+                            t_start = selec3.iloc[0]['t']
+                            print(t_start)
+
+                        dfr['tt'] -= t_start
+
+                        plt.plot(dfr['tt'], dfr['l_extended'], '-.', color=colors[scenario])
+
+                        dfr = dfr[~dfr['mature']]
+                        for _, row in dfr.iterrows():
+                            res_growth.append([row['plantid'], scenario, row['l_extended'], row['tt']])
+
+        res_growth = pd.DataFrame(res_growth, columns=['plantid', 'scenario', 'l', 'tt'])
+
+        # interpolation
+        res_gxe = []
+        t_range = np.linspace(min(res_growth['tt']), max(res_growth['tt']), 20)
+        for plantid in res_growth['plantid'].unique():
+            selec = res_growth[res_growth['plantid'] == plantid]
+            f = interp1d(selec['tt'], selec['l'])
+            for t in [t for t in t_range if min(selec['tt']) <= t <= max(selec['tt'])]:
+                res_gxe.append([plantid, selec.iloc[0]['scenario'], float(f(t)), t])
+        res_gxe = pd.DataFrame(res_gxe, columns=['plantid', 'scenario', 'l', 'tt'])
+
+        res_gxe['l'] /= 10
+        res_gxe = res_gxe.groupby(['scenario', 'tt']).agg(['mean', 'median', 'count', ('mad', mad)]).reset_index()
+        res_gxe = res_gxe[res_gxe['l']['count'] > 2]  # at least 3 rep
+
+        #res_gxe = res_gxe[res_gxe['tt'] <= min(res_gxe.groupby('scenario').max('tt')['tt'])]  # same end for each scenario
+
+        fig, ax = plt.subplots(figsize=(5, 6), dpi=150)
+        fig.canvas.set_window_title(genotype + '_growth')
+        ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+        for scenario in ['WW', 'WD']:
+
+            selec = res_gxe[res_gxe['scenario'] == scenario]
+
+            plt.plot(selec['tt'], selec['l']['median'], '-', color=colors[scenario], markersize=12)
+            plt.plot(selec['tt'], selec['l']['median'], '.', color=colors[scenario], label=scenario, markersize=12)
+            plt.fill_between(selec['tt'],
+                             selec['l']['median'] - selec['l']['mad'] / 1,
+                             selec['l']['median'] + selec['l']['mad'] / 1,
+                             color=colors['error'])
+        # lgd = plt.legend(prop={'size': 20}, loc='upper left', markerscale=2, handletextpad=-0.5)
+        # lgd.get_frame().set_linewidth(0.0)
+        #plt.xlabel('Thermal time after reaching 70% of max length (°C)', fontsize=20)
+        plt.xlabel('Thermal time after leaf emergence $(day_{20°C})$', fontsize=25)
+        plt.ylabel('Leaf length (cm)', fontsize=30) #labelpad=-5)
+
+        ax.text(0.6, 0.13, 'Leaf {}'.format(RANK), transform=ax.transAxes,
+                fontsize=30,
+                verticalalignment='top')
+
+        plt.xlim((-2, 27))
+        plt.ylim((25, 125))
+
+        fig.savefig('paper/growth_{}_{}'.format(RANK, genotype), bbox_inches='tight')
+        plt.close('all')
+
+
+
diff --git a/paper/leaf_count.py b/paper/leaf_count.py
new file mode 100644
index 0000000..ea05e91
--- /dev/null
+++ b/paper/leaf_count.py
@@ -0,0 +1,382 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import time
+import datetime
+import os
+from matplotlib.ticker import MaxNLocator
+
+from sklearn.metrics import r2_score
+from scipy.interpolate import interp1d
+
+from openalea.maizetrack.stem_correction import savgol_smoothing_function
+
+from openalea.maizetrack.stem_correction import stem_height_smoothing
+
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+#from llorenc, currently copied to  modulor cache
+df_tt = pd.read_csv('TT_ZA17.csv')
+f_tt = interp1d(df_tt['timestamp'], df_tt['ThermalTime'])
+
+#from llorenc, currently copied to  modulor cache
+df_pheno = pd.read_csv('pheno_ZA17.csv')
+df_pheno['timestamp'] = df_pheno['resultdate'].map(date_to_timestamp)
+
+# available in modulor local_benoit
+folder = 'data/tracking/'
+folder = 'data/tracking/results_tracking_360/'
+#folder = 'data/tracking/gxe/'
+plantids = [int(f.split('.')[0]) for f in os.listdir(folder) if os.path.isfile(folder + f)]
+obs_visi, pred_visi, obs_ligul, pred_ligul, pred_visi_smooth, pred_ligul_smooth = [], [], [], [], [], []
+pred_visi_smooth2 = []
+pred_ligul2 = []
+
+res_emergence = []
+
+for plantid in plantids:
+
+    print(plantid)
+
+    # ===== data ===================================================================================
+
+    dfid = pd.read_csv(folder + '{}.csv'.format(plantid))
+    dfid_mature = dfid[dfid['mature']]
+
+    # ===========================================================================================
+    # ==== VISIBLE STAGE ========================================================================
+    # ===========================================================================================
+
+    # ===== brut visible stage ==================================================================
+
+    #select = df_phenoid[(df_phenoid['observationcode'] == 'panicule_deployee') & (df_phenoid['observation'] == 1)]
+    df_phenoid = df_pheno[df_pheno['Pot'] == plantid]
+    select = df_phenoid[(df_phenoid['observationcode'] == 'panicule_visi') & (df_phenoid['observation'] == 1)]
+    t_max = max(dfid['timestamp'])
+    if not select.empty:
+        t_max = min(select['timestamp'])
+
+    timestamps = [t for t in dfid['timestamp'].unique() if t < t_max]
+    df_visi = []
+    for t in timestamps:
+
+        dfidt = dfid[dfid['timestamp'] == t]
+
+        mature_ranks = dfidt[dfidt['mature']]['rank_tracking']
+        if mature_ranks.empty:
+            n_mature = 0
+        else:
+            n_mature = max(mature_ranks)
+        n_growing = len(dfidt[dfidt['mature'] == False])
+
+        df_visi.append([t, n_mature + n_growing])
+    df_visi = pd.DataFrame(df_visi, columns=['t', 'n'])
+
+    # ===== visible emergence timing (06/10 idee Christian) ===========================================
+
+    median_visi = df_visi.groupby('n').median('t').reset_index()
+
+    # add missing n values (interpolation)
+    for n in range(min(median_visi['n']) + 1, max(median_visi['n'])):
+        if n not in list(median_visi['n']):
+            t1 = median_visi[median_visi['n'] < n].sort_values('n').iloc[-1]['t']
+            t2 = median_visi[median_visi['n'] > n].sort_values('n').iloc[0]['t']
+            median_visi = median_visi.append({'n': n, 't': (t1 + t2) / 2}, ignore_index=True)
+
+    # compute timing of leaf appearance (except for first and last n)
+    emerg_visi = []
+    for n in sorted(list(median_visi['n'].unique()))[1:-1]:
+        t = np.mean(median_visi[median_visi['n'].isin([n - 1, n])]['t'])
+        emerg_visi.append([plantid, t, n, 'visi'])
+    emerg_visi = pd.DataFrame(emerg_visi, columns=['plantid', 't', 'n', 'type'])
+
+    # force monotony : t = f(n) can only increase
+    emerg_visi = emerg_visi.sort_values('n')
+    for i_row, row in emerg_visi.iterrows():
+        previous = emerg_visi[emerg_visi['n'] < row['n']]
+        if not previous.empty:
+            emerg_visi.iloc[i_row, emerg_visi.columns.get_loc('t')] = max(row['t'], max(previous['t']))
+
+    res_emergence.append(emerg_visi)
+
+    # ===========================================================================================
+    # ==== LIGULATED STAGE ======================================================================
+    # ===========================================================================================
+
+    # ===== profile height used for ligulated stage =================================================
+
+    ranks = sorted([r for r in dfid_mature['rank_tracking'].unique() if r != 0])
+    height_profile = {}
+    for r in ranks:
+        dfr = dfid_mature[dfid_mature['rank_tracking'] == r]
+        dfr = dfr.sort_values('timestamp')[:15]
+        height_profile[r] = np.median(dfr['h']) + 720
+
+    # ===== stem height used for ligulated stage =====================================================
+
+    dfs = []
+    folder_dl = 'local_cache/cache_ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+    #df_tt = pd.read_csv('TT_ZA17.csv')
+    all_files = [folder_dl + '/' + rep + '/' + f for rep in os.listdir(folder_dl) for f in os.listdir(folder_dl + '/' + rep)]
+    plantid_files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+    for f in plantid_files:
+        daydate = f.split('/')[-2]
+        timestamp = df_tt[df_tt['daydate'] == daydate].iloc[0]['timestamp']
+        dft = pd.read_csv(f)
+        dft['t'] = timestamp
+        dfs.append(dft)
+
+    if dfs == []:
+        print('no collar data {}'.format(plantid))
+    else:
+
+        df_collars = pd.concat(dfs)
+
+        df_stem = df_collars.sort_values('z_phenomenal', ascending=False).drop_duplicates(['t']).sort_values('t')
+
+        f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['z_phenomenal']))
+        df_stem['z_smooth'] = [f_smoothing(t) for t in df_stem['t']]
+
+        # ===== ligulated emergence timing =============================================================================
+
+        T = np.linspace(min(df_stem['t']), max(df_stem['t']), 3000)
+        Z = np.array([f_smoothing(ti) for ti in T])
+
+        emerg_ligu = []
+        ranks = sorted(height_profile.keys())
+        for r in ranks[1:]:
+            t = T[np.argmin(np.abs(Z - height_profile[r - 1]))]
+            emerg_ligu.append([plantid, t, r, 'ligu'])
+        emerg_ligu = pd.DataFrame(emerg_ligu, columns=['plantid', 't', 'n', 'type'])
+
+        res_emergence.append(emerg_ligu)
+
+    # ===================================================================================================
+    # ===================================================================================================
+    # ===================================================================================================
+
+    # ===== plot example for one plant ==================================================================
+
+    #if plantid == 461:
+    if False:
+
+        # VISI
+
+        fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+        ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+        plt.title(plantid)
+        plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+        plt.ylabel('Visible leaf stage', fontsize=30)
+
+        # pred
+        plt.plot([f_tt(t) for t in df_visi['t']], df_visi['n'], 'k.', markersize=15, label='raw prediction')
+        plt.plot([f_tt(t) for t in emerg_visi['t']], emerg_visi['n'], 'r.-', markersize=15, label='prediction after smoothing')
+
+        # anot
+        # df_phenoid = df_pheno[df_pheno['Pot'] == plantid]
+        # anot = df_phenoid[df_phenoid['observationcode'] == 'visi_f']
+        # plt.plot([f_tt(t) for t in anot['timestamp']], anot['observation'], 'g^', markersize=7, label='observation (ground-truth)')
+
+        plt.legend(prop={'size': 20}, loc='lower right', markerscale=2)
+
+
+    if False:
+
+        # LIGU
+
+        fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+        ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+        plt.title(plantid)
+        plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+        plt.ylabel('Ligulated leaf stage', fontsize=30)
+
+        # pred
+        plt.plot([f_tt(t) for t in emerg_ligu['t']], emerg_ligu['n'], 'k.', markersize=15, label='leaf emergence')
+        plt.step([f_tt(t) for t in emerg_ligu['t']], emerg_ligu['n'], 'r-', where='post', label='prediction (discrete interpolation)')
+        # plt.plot([f_tt(t) for t in emerg_ligu['t']], emerg_ligu['n'], 'b-', label='linear interpolation')
+
+        # anot
+        df_phenoid = df_pheno[df_pheno['Pot'] == plantid]
+        anot = df_phenoid[df_phenoid['observationcode'] == 'ligul_f']
+        plt.plot([f_tt(t) for t in anot['timestamp']], anot['observation'], 'g^', markersize=7, label='observation (ground-truth)')
+
+        plt.legend(prop={'size': 20}, loc='lower right', markerscale=2)
+
+
+    # ======================================================================================================
+
+res_emergence = pd.concat(res_emergence)
+
+
+# ===== from emergence timing to leaf stage + anot =========================================================
+
+res_stage = []
+
+for plantid in res_emergence['plantid'].unique():
+    for var in res_emergence['type'].unique():
+
+        df_phenoid = df_pheno[df_pheno['Pot'] == plantid]
+
+        selec = res_emergence[(res_emergence['plantid'] == plantid) & (res_emergence['type'] == var)]
+
+        if selec.empty:
+            print('no data {} _ {}'.format(plantid, var))
+        if not selec.empty:
+
+            if var == 'visi':
+                anot = df_phenoid[df_phenoid['observationcode'] == 'visi_f']
+                f = interp1d(selec['t'], selec['n'])
+            elif var == 'ligu':
+                anot = df_phenoid[df_phenoid['observationcode'] == 'ligul_f']
+                f = interp1d(selec['t'], selec['n'], kind='previous')
+
+            f_stage = lambda t: f(t) if min(selec['t']) < t < max(selec['t']) else None
+
+            for _, row in anot.iterrows():
+                obs = row['observation']
+                pred = f_stage(row['timestamp'])
+                #t_nearest = min(leaf_count.keys(), key=lambda x: abs(x - row['timestamp']))
+                res_stage.append([plantid, obs, pred, var])
+
+# # ===== add phenomenal brut result (2022) ==================================================
+
+res_stage = [k + ['new'] for k in res_stage]
+df_phm = pd.read_csv('paper/leaf_stage_phm2.csv')
+
+res_stage = pd.DataFrame(res_stage, columns=['plantid', 'obs', 'pred', 'type', 'method'])
+res_stage = pd.concat((df_phm, res_stage))
+
+res_stage.to_csv('paper/leaf_count.csv', index=False)
+
+# =========================================================================================
+
+res_stage = pd.read_csv('paper/leaf_count.csv')
+
+# ===== plot visible leaf stage ===========================================================
+
+s = res_stage[(res_stage['type'] == 'visi') & (~res_stage['pred'].isna())]
+s1 = s[s['method'] == 'new']
+x1, y1 = np.array(s1['obs']), np.array(s1['pred'])
+#x = np.ceil(x - 1)
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+#ax.tick_params(axis='both', which='major', labelsize=20)
+ax.xaxis.set_major_locator(MaxNLocator(integer=True)) # int axis
+plt.title('Visible leaf stage', fontsize=35)
+plt.xlabel('observation', fontsize=30)
+plt.ylabel('prediction', fontsize=30)
+plt.plot([-5, 25], [-5, 25], '-', color='black')
+plt.xlim((0, 18))
+plt.ylim((0, 18))
+
+s2 = s[s['method'] == 'phm']
+x2, y2 = np.array(s2['obs']), np.array(s2['pred'])
+plt.plot(x2, y2, '.', color='grey', markersize=8, alpha=0.2)
+
+plt.plot(x1, y1, 'o', color='red', markersize=8, alpha=0.1)
+
+# a0, a, b = np.polyfit(x, y, 2)
+# xu = np.array(sorted(np.unique(x)))
+# plt.plot(xu, a0 * np.array(xu**2) + a * np.array(xu) + b, 'r-')
+# x = x[[i for i in range(len(y)) if y[i] is not None]]
+# y = y[[i for i in range(len(y)) if y[i] is not None]]
+a, b = np.polyfit(x1, np.array([float(yi) for yi in y1]), 1)
+plt.plot([-10, 30], a*np.array([-10, 30]) + b, '--', color='red', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+a, b = np.polyfit(x2, np.array([float(yi) for yi in y2]), 1)
+plt.plot([-10, 30], a*np.array([-10, 30]) + b, '--', color='grey', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+# rmax = 10
+# x2, y2 = x[np.where(x <= rmax)], y[np.where(x <= rmax)]
+# a2, b2 = np.polyfit(x2, y2, 1)
+# plt.plot(x2, a2*x2 + b2, 'g-', label='linear regression (observation < {}) \n (a = {}, b = {})'.format(rmax, round(a2, 3), round(b2, 2)))
+leg = plt.legend(prop={'size': 22}, loc='upper left')
+leg.get_texts()[1].set_color('grey')
+leg.get_frame().set_linewidth(0.0)
+
+# set the linewidth of each legend object
+# for legobj in leg.legendHandles:
+#     legobj.set_linewidth(3)
+
+rmse = np.sqrt(np.sum((x1 - y1) ** 2) / len(x1))
+r2 = r2_score(x1, y1)
+biais = np.mean(y1 - x1)
+mape = 100 * np.mean(np.abs((x1 - y1) / x1))
+ax.text(0.60, 0.26, 'n = {} \nBias = {}\nRMSE = {}\nMAPE = {} % \nR² = {}'.format(
+    len(x1), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)),
+        transform=ax.transAxes, fontsize=25, verticalalignment='top')
+
+rmse = np.sqrt(np.sum((x2 - y2) ** 2) / len(x2))
+r2 = r2_score(x2, y2)
+biais = np.mean(y2 - x2)
+mape = 100 * np.mean(np.abs((x2 - y2) / x2))
+ax.text(0.05, 0.75, 'n = {} \nBias = {}\nRMSE = {} \nMAPE = {} % \nR² = {}'.format(
+    len(x2), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)), transform=ax.transAxes, fontsize=20,
+        verticalalignment='top', color='grey')
+
+fig.savefig('paper/results/nvisi_v2', bbox_inches='tight')
+
+# ===== plot ligulated leaf stage ============================================================
+
+s = res_stage[(res_stage['type'] == 'ligu') & (~res_stage['pred'].isna())]
+s1 = s[s['method'] == 'new']
+x1, y1 = np.array(s1['obs']), np.array(s1['pred'])
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+#ax.tick_params(axis='both', which='major', labelsize=20)
+ax.xaxis.set_major_locator(MaxNLocator(integer=True)) # int axis
+plt.title('Ligulated leaf stage', fontsize=35)
+plt.xlabel('observation', fontsize=30)
+plt.ylabel('prediction', fontsize=30)
+plt.plot([-5, 25], [-5, 25], '-', color='black')
+plt.xlim((0, 18))
+plt.ylim((0, 18))
+
+s2 = s[s['method'] == 'phm']
+x2, y2 = np.array(s2['obs']), np.array(s2['pred'])
+plt.plot(x2, y2, '.', color='grey', markersize=8, alpha=0.2)
+
+plt.plot(x1, y1, 'ro', markersize=8, alpha=0.05)
+
+a, b = np.polyfit(x1, np.array([float(yi) for yi in y1]), 1)
+plt.plot([-10, 30], a*np.array([-10, 30]) + b, 'r--', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+a, b = np.polyfit(x2, np.array([float(yi) for yi in y2]), 1)
+plt.plot([-10, 30], a*np.array([-10, 30]) + b, '--', color='grey', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+leg = plt.legend(prop={'size': 22}, loc='upper left', markerscale=2)
+leg.get_frame().set_linewidth(0.0)
+leg.get_texts()[1].set_color('grey')
+# set the linewidth of each legend object
+# for legobj in leg.legendHandles:
+#     legobj.set_linewidth(3)
+
+rmse = np.sqrt(np.sum((x1 - y1) ** 2) / len(x1))
+r2 = r2_score(x1, y1)
+biais = np.mean(y1 - x1)
+mape = 100 * np.mean(np.abs((x1 - y1) / x1))
+ax.text(0.60, 0.26, 'n = {} \nBias = {}\nRMSE = {} \nMAPE = {} % \nR² = {}'.format(
+    len(x1), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)), transform=ax.transAxes, fontsize=25,
+        verticalalignment='top')
+
+rmse = np.sqrt(np.sum((x2 - y2) ** 2) / len(x2))
+r2 = r2_score(x2, y2)
+biais = np.mean(y2 - x2)
+mape = 100 * np.mean(np.abs((x2 - y2) / x2))
+ax.text(0.05, 0.75, 'n = {} \nBias = {}\nRMSE = {} \nMAPE = {} % \nR² = {}'.format(
+    len(x2), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)), transform=ax.transAxes, fontsize=20,
+        verticalalignment='top', color='grey')
+
+fig.savefig('paper/results/nligu_v2', bbox_inches='tight')
+
+
+# histo
+
+fig, ax = plt.subplots(figsize=(10, 6), dpi=100)
+plt.hist(y - x, bins=len(set(y-x)), stacked=True, density=True, edgecolor='k', color='grey')
+ax.axes.xaxis.set_visible(False)
+ax.axes.yaxis.set_visible(False)
diff --git a/paper/leaf_extension_example.py b/paper/leaf_extension_example.py
new file mode 100644
index 0000000..632ccb7
--- /dev/null
+++ b/paper/leaf_extension_example.py
@@ -0,0 +1,43 @@
+from skimage import io
+import os
+
+from openalea.phenomenal import object as phm_obj
+import matplotlib.pyplot as plt
+
+from openalea.maizetrack.local_cache import metainfos_to_paths, get_metainfos_ZA17
+from openalea.maizetrack.leaf_extension import leaf_extension
+from openalea.maizetrack.phenomenal_display import plot_sk, plot_vmsi, plot_snapshot
+
+
+if __name__ == '__main__':
+
+    plantid = 1424
+
+    metainfos = get_metainfos_ZA17(plantid)
+
+    daydate = '2017-05-12'
+
+    metainfo = next(m for m in metainfos if m.daydate == daydate)
+
+    vmsi_path = metainfos_to_paths([metainfo], object='vmsi', stem_smoothing=True, old=True)[0]
+    vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(vmsi_path)
+
+    shooting_frame = metainfo.shooting_frame
+
+    # binary image
+    binaries = dict()
+    path = 'binary2/' + daydate
+    files = [f for f in os.listdir(path) if f[:4] == str(plantid).rjust(4, '0')]
+    files = [f for f in files if str(metainfo.task) in f]
+    angles = sorted([a for a, t in zip(metainfo.camera_angle, metainfo.view_type) if t == 'side'])
+    for angle in angles:
+        file = next(f for f in files if str(angle) in f[-7:])
+        binaries[angle] = io.imread(path + '/' + file)  # / 255. later
+
+    vmsi = leaf_extension(vmsi, binaries, shooting_frame, display_parameters=[60, True, False, False])
+
+    skeleton_path = metainfos_to_paths([metainfo], object='skeleton')[0]
+    skeleton = phm_obj.VoxelSkeleton.read_from_json_gz(skeleton_path)
+
+    plot_sk(skeleton)
+
diff --git a/paper/leaf_extension_validation.py b/paper/leaf_extension_validation.py
new file mode 100644
index 0000000..55ce8d1
--- /dev/null
+++ b/paper/leaf_extension_validation.py
@@ -0,0 +1,107 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+import time
+import datetime
+import os
+
+from sklearn.metrics import r2_score
+
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+# availble on NasShare2
+df_anot = pd.read_csv('data/rgb_leaf_annotation/all_anot.csv')  # attention 940
+
+# ====== leaf extension validation ==================================
+
+# mature data, 9 plants
+
+res = []
+
+for plantid in [int(p) for p in df_anot['plantid'].unique()]:
+    dfid = pd.read_csv('data/tracking/{}.csv'.format(plantid)) # modulor/local_benoit
+    df_anotid = df_anot[df_anot['plantid'] == plantid].sort_values('rank')
+
+    if plantid == 940:
+        df_anotid = df_anotid.groupby('rank').max('timestamp').reset_index()
+
+    for _, row in df_anotid.iterrows():
+        select = dfid[(dfid['rank_annotation'] == row['rank']) &
+                             (dfid['timestamp'] == row['timestamp'])]
+        if not select.empty:
+            select_row = select.iloc[0]
+            res.append([row['length'], select_row['l'], select_row['l_extended'], 'mature'])
+        else:
+            print('{} 1 row with no data'.format(plantid))
+
+# growing data, 10 plants (rank 6, rank 9)
+# availble on NasShare2
+df_anot2 = pd.read_csv('data/rgb_leaf_growing_annotation/leaf_growing_annotation.csv')
+
+for plantid in [int(p) for p in df_anot2['plantid'].unique()]:
+    dfid = pd.read_csv('data/tracking/{}.csv'.format(plantid))
+    df_anotid = df_anot2[df_anot2['plantid'] == plantid].sort_values('rank')
+
+    for _, row in df_anotid.iterrows():
+        select = dfid[(dfid['rank_annotation'] == row['rank']) &
+                             (dfid['timestamp'] == row['timestamp'])]
+        if not select.empty:
+            select_row = select.iloc[0]
+            res.append([row['length_mm'], select_row['l'], select_row['l_extended'], 'growing'])
+        else:
+            print('{} 1 row with no data'.format(plantid))
+
+
+res = pd.DataFrame(res, columns=['obs', 'pred_phm', 'pred_ext', 'var'])
+
+# ===== everything in the same graph ==================================================================
+
+# mature + growing
+# x = np.concatenate((obs, obs2))
+# y = np.concatenate((pred_phm, pred2_phm))
+# y_ext = np.concatenate((pred_ext, pred2_ext))
+
+# only_mature
+x = res['obs'] / 10
+y = res['pred_phm'] / 10
+y_ext = res['pred_ext'] / 10
+
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+r2 = r2_score(x, y)
+bias = np.mean(x - y)
+rmse_ext = np.sqrt(np.sum((x - y_ext) ** 2) / len(x))
+r2_ext = r2_score(x, y_ext)
+bias_ext = np.mean(x - y_ext)
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.xlim((0, 135))
+plt.ylim((0, 135))
+plt.plot([-10, 150], [-10, 150], '-', color='k')
+plt.gca().set_aspect('equal', adjustable='box')
+plt.xlabel('observation (cm)', fontsize=30)
+plt.ylabel('prediction (cm)', fontsize=30)
+plt.plot(x, y, '^', color='grey', markersize=6, label='Phenomenal') # \n RMSE = {} cm, R² = {}'.format(round(rmse, 1), round(r2, 3)))
+plt.plot(x, y_ext, 'o', color='black', markersize=6, label='Phenomenal with leaf extension') # \n RMSE = {} cm, R² = {}'.format(round(rmse_ext, 1), round(r2_ext, 3)))
+leg = plt.legend(prop={'size': 16}, loc='upper left', markerscale=2)
+leg.get_frame().set_linewidth(0.0)
+plt.title('Leaf length', fontsize=35)
+#plt.subtitle('Leaf length', fontsize=35, y=0.94)
+
+
+ax.text(0.3, 0.17, 'n = {} \nBias = {} cm \nRMSE = {} cm \nR² = {}'.format(len(x), round(bias_ext, 2), round(rmse_ext, 2), round(r2_ext, 3)), transform=ax.transAxes, fontsize=20,
+        verticalalignment='top', color='black')
+ax.text(0.67, 0.17, 'n = {} \nBias = {} cm \nRMSE = {} cm \nR² = {}'.format(len(x), round(bias, 2), round(rmse, 1), round(r2, 3)), transform=ax.transAxes, fontsize=20,
+        verticalalignment='top', color='grey')
+
+fig.savefig('paper/results/extension_validation', bbox_inches='tight')
+
+
+
+
+
+
+
diff --git a/paper/narcolepsie.py b/paper/narcolepsie.py
new file mode 100644
index 0000000..85dca85
--- /dev/null
+++ b/paper/narcolepsie.py
@@ -0,0 +1,150 @@
+from openalea.maizetrack.utils import phm3d_to_px2d
+from openalea.maizetrack.stem_correction import xyz_last_mature
+
+from openalea.maizetrack.local_cache import metainfos_to_paths, get_metainfos_ZA17, check_existence, load_plant
+from openalea.maizetrack.utils import phm3d_to_px2d
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.maizetrack.stem_correction import stem_height_smoothing
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import json
+import os
+from scipy.interpolate import interp1d
+from skimage import io
+import cv2
+
+df_tt = pd.read_csv('TT_ZA17.csv')
+f_tt = interp1d(df_tt['timestamp'], df_tt['ThermalTime'])
+
+# ============================================================================
+
+# for k in range(300, 600):
+#     img = io.imread('data/rgb_insertion_annotation/train2/img{}.png'.format(k))
+#     anot = np.loadtxt('data/rgb_insertion_annotation/train2/img{}.txt'.format(k))
+#     if anot.shape == (5,):
+#         anot = np.array([anot])
+#     for _, x, y, _, _ in anot:
+#
+#         img = cv2.circle(img, (int(x * 416), int(y * 416)), 6, (255, 0, 0), -1)
+#     io.imsave('Narcolepsie/dl/dl{}.png'.format(k), img)
+#
+
+# ============================================================================
+
+#df_tt = pd.read_csv('TT_ZA17.csv')
+
+plantid = 1429
+
+res = []
+
+# ===== phenomenal brut ======
+
+folder = 'local_cache/ZA17/1429_nosmooth/'
+
+metainfos = get_metainfos_ZA17(plantid)
+paths = metainfos_to_paths(metainfos, folder=folder)
+metainfos2, paths = check_existence(metainfos, paths)
+vmsi_list = load_plant(metainfos2, paths)
+
+for vmsi in vmsi_list:
+    if vmsi.get_mature_leafs():
+        i = np.argmax([l.info['pm_z_base'] for l in vmsi.get_mature_leafs()])
+        xyz = vmsi.get_mature_leafs()[i].info['pm_position_base']
+        xyz[2] = vmsi.get_mature_leafs()[i].info['pm_z_base_voxel']
+        y = phm3d_to_px2d(xyz, vmsi.metainfo.shooting_frame, angle=60)[0][1]
+        t = vmsi.metainfo.timestamp
+        res.append([t, 2448 - y, 'phm'])
+
+# ===== deep learning =====
+
+folder = 'local_cache/ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+dfs = []
+
+plantid_files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+for f in plantid_files:
+    daydate = f.split('/')[3]
+    dft = pd.read_csv(f)
+    #dft['t'] = df_tt[df_tt['daydate'] == daydate].iloc[0]['ThermalTime']
+    dft['t'] = next(m for m in metainfos2 if m.daydate == daydate).timestamp
+    dfs.append(dft)
+df = pd.concat(dfs)
+
+df = df[df['score'] > 0.95]
+df['y'] = 2448 - df['y']
+
+fig, ax = plt.subplots()
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.plot(df['t'], df['y'], 'k.', markersize=2)
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+plt.ylabel('collar height (pixels)', fontsize=30)
+
+df_stem = df.sort_values('y', ascending=False).drop_duplicates(['t']).sort_values('t')
+plt.plot(df_stem['t'], df_stem['y'], 'r-', label='stem height')
+f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['y']))
+plt.plot(df_stem['t'], [f_smoothing(t) for t in df_stem['t']], 'r--', label='stem height smoothing')
+plt.legend()
+
+for _, row in df_stem.iterrows():
+    res.append([row['t'], row['y'], 'dl'])
+
+# plt.figure(plantid)
+# plt.plot(df['t'], df['y'], 'k*', markersize=2)
+# plt.plot(df_stem['t'], df_stem['y'], '-b')
+# plt.plot(df_stem['t'], [f_smoothing(t) for t in df_stem['t']], '-', color='orange')
+
+# ====== stem height annotation =====================================
+
+with open('data/stem_annotation/stem_{}.json'.format(plantid)) as f:
+    d = json.load(f)
+
+    for name in d.keys():
+
+        task = int(name.split('t')[1].split('_')[0])
+        timestamp = next(m for m in metainfos if m.task == task).timestamp
+        angle = int(name.split('a')[1].split('.')[0])
+        sf = next(m for m in metainfos if m.task == task).shooting_frame
+
+        for shape in d[name]['regions']:
+
+                y = shape['shape_attributes']['cy']
+                #tt = df_tt.iloc[np.argmin(np.abs(df_tt['timestamp'] - timestamp))]['ThermalTime']
+                res.append([timestamp, 2448 - y, 'anot'])
+
+
+# ===== annotation =====
+res = pd.DataFrame(res, columns=['t', 'y', 'var'])
+res['tt'] = [f_tt(t) for t in res['t']]
+
+selec = res[res['var'] == 'dl'].sort_values('tt')
+plt.plot(selec['tt'], selec['y'], 'g-')
+selec = res[res['var'] == 'anot'].sort_values('tt')
+plt.plot(selec['tt'], selec['y'], 'r-')
+selec = res[res['var'] == 'phm'].sort_values('tt')
+plt.plot(selec['tt'], selec['y'], 'k-')
+
+
+
+selec = res[res['var'] == 'phm'].sort_values('t')
+plt.plot(selec['tt'], (selec['y'] - np.min(selec['y'])) / 10, 'k-')
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=20)
+plt.ylabel('Stem height (cm)', fontsize=20)
+
+
+
+
+
+df_stem['tt'] = [f_tt(t) for t in df_stem['t']]
+plt.plot(df_stem['tt'], np.array([f_smoothing(t) for t in df_stem['t']]) - 530, 'r-', label='stem height smoothing')
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=20)
+plt.ylabel('Stem height (cm)', fontsize=20)
+
+
+
+
+
+
+
diff --git a/paper/pipeline_overview.py b/paper/pipeline_overview.py
new file mode 100644
index 0000000..18fda3b
--- /dev/null
+++ b/paper/pipeline_overview.py
@@ -0,0 +1,155 @@
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.utils import get_rgb
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.phenomenal import object as phm_obj
+from openalea.maizetrack.phenomenal_display import plot_vmsi, plot_vmsi_voxel, plot_vg, plot_snapshot
+import copy
+
+import openalea.phenomenal.segmentation as phm_seg
+from openalea.maizetrack.trackedPlant import TrackedSnapshot
+
+# ===============================
+
+plantid = 474
+
+metainfos = get_metainfos_ZA17(plantid)
+
+metainfos = [m for m in metainfos if m.daydate != '2017-04-20']  # debug
+
+paths = metainfos_to_paths(metainfos, stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)
+
+metainfos, paths = check_existence(metainfos, paths)
+
+vmsi_list = load_plant(metainfos, paths)
+plant = TrackedPlant.load_and_check(vmsi_list)
+
+plant.align_mature(direction=1, gap=12.365, w_h=0.03, gap_extremity_factor=0.2, n_previous=500)
+
+plant_growing = copy.deepcopy(plant)
+plant_growing.align_growing()
+
+# for m in [m for m in metainfos if '2017-05-23' <= m.daydate <= '2017-05-23']:
+#     get_rgb(m, 150, main_folder='data/paper', plant_folder=True, save=True, side=True)
+#
+# for m in [m for m in metainfos if m.daydate in ['2017-04-14', '2017-05-02']]:
+#     get_rgb(m, 60, main_folder='data/paper', plant_folder=True, save=True, side=True)
+
+print([m.daydate for m in metainfos if m.daydate not in [s.metainfo.daydate for s in plant.snapshots]])
+path = [p for p in paths if '06-04' in p][0]
+vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(path)
+plot_vmsi([vmsi])
+
+plot_vmsi_voxel(vmsi)
+
+
+
+
+# ====== t=1 :
+
+d = '2017-04-14'
+s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+plot_snapshot(s, colored=False)  # mature vs growing
+plot_snapshot(s, colored=True)  # tracking mature
+plot_snapshot(sm, colored=True)  # tracking growing
+
+d = '2017-05-02'
+s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+plot_snapshot(s, colored=False)
+plot_snapshot(s, colored=True)
+plot_snapshot(sm, colored=True)
+
+# d = '2017-05-22'
+# s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+# sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+# plot_snapshot(s, colored=False)
+# plot_snapshot(s, colored=True)
+# plot_snapshot(sm, colored=True, ranks=[2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 10])
+
+# 05-23 : tige correcte, tige trop basse
+d = '2017-05-23'
+s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+plot_snapshot(s, colored=False)
+plot_snapshot(s, colored=True)
+plot_snapshot(s, colored=True, ranks=[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])
+
+# # 05-24 bug tige (artificial)
+# d = '2017-05-24'
+# m = next(m for m in metainfos if m.daydate == d)
+# path = \
+# metainfos_to_paths([m], object='voxelgrid', stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+# vx = phm_obj.VoxelGrid.read_from_csv(path)
+# path = \
+# metainfos_to_paths([m], object='skeleton', stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+# sk = phm_obj.VoxelSkeleton.read_from_json_gz(path)
+# graph = phm_seg.graph_from_voxel_grid(vx, connect_all_point=True)
+# vms = phm_seg.maize_segmentation(sk, graph, z_stem=100)
+# vmsi = phm_seg.maize_analysis(vms)
+# s = TrackedSnapshot(vmsi, m, order=None)
+# plot_snapshot(s, colored=False)
+#
+# d = '2017-05-24'
+# s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+# sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+# plot_snapshot(s, colored=False)
+# plot_snapshot(s, colored=True)
+# plot_snapshot(s, colored=True, ranks=[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])
+
+d = '2017-06-04'  # bug tige 2
+m = next(m for m in metainfos if m.daydate == d)
+path = metainfos_to_paths([m], stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(path)
+s = TrackedSnapshot(vmsi, m, order=None)
+plot_snapshot(s, colored=False)
+
+d = '2017-06-14'
+s = next(s for s in plant.snapshots if s.metainfo.daydate == d)
+sm = next(s for s in plant_growing.snapshots if s.metainfo.daydate == d)
+plot_snapshot(s, colored=False)
+plot_snapshot(s, colored=True)
+plot_snapshot(sm, colored=True)
+
+# ====== T : 06-04 (deformation), 06-03
+
+for k, col in enumerate(PALETTE[:20]):
+    plt.plot(0, k, '*', c=col / 255.)
+
+# ======== voxelgrid ######### -> 06-14
+d = '2017-04-14'
+m = [m for m in metainfos if m.daydate == d][0]
+print(len([m for m in metainfos if m.daydate == d]))
+path = \
+metainfos_to_paths([m], object='voxelgrid', stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+vx = phm_obj.VoxelGrid.read_from_csv(path)
+plot_vg(vx)
+
+# bug tige (artificial)
+d = '2017-05-02'
+m = [m for m in metainfos if m.daydate == d][1]
+path = \
+metainfos_to_paths([m], object='voxelgrid', stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+vx = phm_obj.VoxelGrid.read_from_csv(path)
+path = \
+metainfos_to_paths([m], object='skeleton', stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False)[0]
+sk = phm_obj.VoxelSkeleton.read_from_json_gz(path)
+graph = phm_seg.graph_from_voxel_grid(vx, connect_all_point=True)
+vms = phm_seg.maize_segmentation(sk, graph, z_stem=50)
+vmsi = phm_seg.maize_analysis(vms)
+s = TrackedSnapshot(vmsi, m, order=None)
+plot_snapshot(s, colored=False)
+
+# ===== stem shape anomaly example (plantid 330) ========================================================
+
+vx_path = 'local_cache/cache_ZA17/image3d_voxel4_tol1_notop/2017-05-19/0330_DZ_PG_01_ZM4381_WW_Rep_2_06_30_ARCH2017-03-30__2017-05-19__6726.csv'
+sk_path = 'local_cache/cache_ZA17/skeleton_voxel4_tol1_notop_vis4_minpix100/2017-05-19/0330_DZ_PG_01_ZM4381_WW_Rep_2_06_30_ARCH2017-03-30__2017-05-19__6726.json.gz'
+vx = phm_obj.VoxelGrid.read_from_csv(vx_path)
+sk = phm_obj.VoxelSkeleton.read_from_json_gz(sk_path)
+graph = phm_seg.graph_from_voxel_grid(vx, connect_all_point=True)
+
+vms = phm_seg.maize_segmentation(sk, graph, z_stem=260)
+vmsi = phm_seg.maize_analysis(vms)
+
+s = TrackedSnapshot(vmsi, m, order=None)
+plot_snapshot(s, colored=False)
\ No newline at end of file
diff --git a/paper/profiles.py b/paper/profiles.py
new file mode 100644
index 0000000..ba3ced4
--- /dev/null
+++ b/paper/profiles.py
@@ -0,0 +1,212 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+import os
+from sklearn.metrics import r2_score
+
+plantids = [int(f.split('.')[0]) for f in os.listdir('data/tracking') if os.path.isfile('data/tracking/' + f)]
+
+# available in NASHShare2
+df_anot = pd.read_csv('data/rgb_leaf_annotation/all_anot.csv')  # attention 940
+df_anot_h = pd.read_csv('data/rgb_insertion_profile_annotation/annotation/all_anot.csv')
+
+# ===== combine ligulated data of all plantids in a unique df =================================
+
+df = []
+for plantid in plantids:
+
+    # data plantid
+
+    # available in modulor local_benoit
+    dfid = pd.read_csv('data/tracking/{}.csv'.format(plantid))
+    dfid_mature = dfid[dfid['mature']]
+
+    dfid_mature = dfid_mature[dfid_mature['timestamp'] < 1496509351] # 06-03
+
+    ranks = sorted([r for r in dfid_mature['rank_tracking'].unique() if r != 0])
+    for r in ranks:
+        dfr = dfid_mature[(dfid_mature['mature']) & (dfid_mature['rank_tracking'] == r)]
+        dfr = dfr.sort_values('timestamp')[:15]
+        df.append(dfr)
+
+df = pd.concat(df)
+
+# ====== profile + merge with annotation ========================================================================
+
+res = []
+for plantid in df['plantid'].unique():
+
+    # length annotation
+    df_anotid = df_anot[df_anot['plantid'] == plantid].sort_values('rank')
+
+    # height annotation
+    df_anotid_h = df_anot_h[df_anot_h['plantid'] == plantid].sort_values('height_mm')
+    df_anotid_h['rank'] = np.arange(1, len(df_anotid_h) + 1)
+
+    selec = df[df['plantid'] == plantid]
+    for r in selec['rank_tracking'].unique():
+        if r in list(df_anotid['rank']):
+            dfr = selec[selec['rank_tracking'] == r]
+            obs = df_anotid[df_anotid['rank'] == r].sort_values('length').iloc[-1]['length']
+            pred = np.median(dfr['l_extended'])
+            res.append([plantid, obs, pred, 'l', r])
+
+        # height
+        if r in list(df_anotid_h['rank']):
+            dfr = selec[selec['rank_tracking'] == r]
+            obs = df_anotid_h[df_anotid_h['rank'] == r].iloc[0]['height_mm']
+            pred = np.median(dfr['h'])
+            res.append([plantid, obs, pred, 'h', r])
+
+res = pd.DataFrame(res, columns=['plantid', 'obs', 'pred', 'var', 'rank'])
+
+# ===== plot an example for 1 plant ===============================================================================
+
+# length
+#for plantid in [1276, 1301, 940]:
+plantid = 1301
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+
+# anot
+selec = res[(res['plantid'] == plantid) & (res['var'] == 'l')]
+plt.plot(selec['rank'], selec['obs'] / 10, 'g--', linewidth=3, label='observation')
+
+# pred
+selec = df[df['plantid'] == plantid]
+plt.plot(selec['rank_tracking'], selec['l_extended'] / 10, 'k.', markersize=8, label='prediction')
+selec = selec.groupby('rank_tracking').median('l_extended').reset_index()
+plt.plot(selec['rank_tracking'], selec['l_extended'] / 10, 'k-o', markersize=10, linewidth=3, label='Final median prediction')
+
+plt.legend()
+plt.title('Leaf length profile \n example for one plant', fontsize=35)
+plt.xlabel('Leaf rank', fontsize=30)
+plt.ylabel('Leaf length (cm)', fontsize=30)
+leg = plt.legend(prop={'size': 22}, loc='lower center', markerscale=2)
+leg.get_frame().set_linewidth(0.0)
+
+fig.savefig('paper/results/profil_example', bbox_inches='tight')
+
+# fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+# ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+# selec = df[df['plantid'] == plantid]
+# plt.plot(selec['rank_tracking'], (selec['h'] + 700) / 10, 'k.', markersize=10, label='Ligulated leaf')
+# selec = selec.groupby('rank_tracking').median('l_extended').reset_index()
+# plt.plot(selec['rank_tracking'], (selec['h'] + 700) / 10, 'r-*', markersize=10, label='Profile value (median)')
+# plt.legend()
+# plt.xlabel('Rank', fontsize=30)
+# plt.ylabel('Length (cm)', fontsize=30)
+# plt.legend(prop={'size': 22}, loc='lower right', markerscale=2)
+
+# insertion height
+for plantid in [439, 907]:
+
+    fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+    ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+
+    # anot
+    selec = res[(res['plantid'] == plantid) & (res['var'] == 'h')]
+    plt.plot(selec['rank'], selec['obs'] / 10, 'g-', linewidth=3, label='observation (ground-truth)')
+
+    # pred
+    selec = df[df['plantid'] == plantid]
+    selec['h'] += 700
+    plt.plot(selec['rank_tracking'], selec['h'] / 10, 'k.', markersize=10)
+    selec = selec.groupby('rank_tracking').median('l_extended').reset_index()
+    plt.plot(selec['rank_tracking'], selec['h'] / 10, 'r-', markersize=10, linewidth=3, label='prediction (median)')
+
+    plt.legend()
+    plt.title(plantid)
+    plt.xlabel('Rank', fontsize=30)
+    plt.ylabel('Insertion height (cm)', fontsize=30)
+    plt.legend(prop={'size': 22}, loc='lower right', markerscale=2)
+
+# ===== plot validation of profile length ============================================================================
+
+selec = res[res['var'] == 'l']
+x, y = np.array(selec['obs']) / 10, np.array(selec['pred']) / 10
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+plt.title('Length profile', fontsize=35)
+plt.xlabel('observation (cm)', fontsize=30)
+plt.ylabel('prediction (cm)', fontsize=30)
+plt.plot([-20, 300], [-20, 300], '-', color='black')
+plt.xlim((0, 130))
+plt.ylim((0, 130))
+plt.plot(x, y, 'ro', markersize=8, alpha=0.5)
+
+# add results without leaf extension
+df_ext = pd.read_csv('paper/leaf_ext.csv').iloc[:-1]
+x2, y2 = df_ext['obs'] / 10, df_ext['pred_phm'] / 10
+plt.plot(x2, y2, '.', color='grey', markersize=8, alpha=0.5)
+
+a, b = np.polyfit(x, y, 1)
+plt.plot([-10, 3000], a*np.array([-10, 3000]) + b, 'r--', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+a, b = np.polyfit(x2, y2, 1)
+plt.plot([-10, 3000], a*np.array([-10, 3000]) + b, '--', color='grey', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+leg = plt.legend(prop={'size': 22}, loc='upper left', markerscale=2)
+leg.get_frame().set_linewidth(0.0)
+leg.get_texts()[1].set_color('grey')
+
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+r2 = r2_score(x, y)
+biais = np.mean(y - x)
+mape = 100 * np.mean(np.abs((x - y) / x))
+ax.text(0.60, 0.26, 'n = {} \nBias = {} cm\nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)), transform=ax.transAxes, fontsize=25,
+        verticalalignment='top')
+rmse2 = np.sqrt(np.sum((x2 - y2) ** 2) / len(x2))
+r2_2 = r2_score(x2, y2)
+bias2 = np.mean(y2 - x2)
+mape2 = 100 * np.mean(np.abs((x2 - y2) / x2))
+ax.text(0.05, 0.75, 'n = {} \nBias = {} cm \nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x2), round(bias2, 2), round(rmse2, 1), round(mape2, 1), round(r2_2, 3)),
+        transform=ax.transAxes, fontsize=20, verticalalignment='top', color='grey')
+
+fig.savefig('paper/results/profil_l', bbox_inches='tight')
+
+# ===== plot validation of profile height ============================================================================
+
+selec = res[res['var'] == 'h']
+x, y = np.array(selec['obs']) / 10, (np.array(selec['pred']) + 700) / 10
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+plt.title('Insertion height profile', fontsize=35)
+plt.xlabel('observation (cm)', fontsize=30)
+plt.ylabel('prediction (cm)', fontsize=30)
+plt.plot([-20, 300], [-20, 300], '-', color='k')
+plt.xlim((-10, 210))
+plt.ylim((-10, 210))
+plt.plot(x, y, 'ro', markersize=8, alpha=0.5)
+
+a, b = np.polyfit(x, y, 1)
+plt.plot([-10, 3000], a*np.array([-10, 3000]) + b, 'r--', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+
+leg = plt.legend(prop={'size': 22}, loc='upper left', markerscale=2)
+leg.get_frame().set_linewidth(0.0)
+
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+r2 = r2_score(x, y)
+biais = np.mean(y - x)
+mape = 100 * np.mean(np.abs((x - y) / x))
+ax.text(0.60, 0.26, 'n = {} \nBias = {} cm\nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x), round(biais, 2), round(rmse, 2), round(mape, 1), round(r2, 3)), transform=ax.transAxes, fontsize=25,
+        verticalalignment='top')
+
+fig.savefig('paper/results/profil_h', bbox_inches='tight')
+
+
+
+# ===== for review, juste change fig title ========================================================================
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20) # axis number size
+# plt.gca().set_aspect('equal', adjustable='box') # same x y scale
+plt.title('Ligulated leaf length = f(rank)\n- example for one plant', fontsize=26)
+
+
diff --git a/paper/stem_height_px.py b/paper/stem_height_px.py
new file mode 100644
index 0000000..d16a382
--- /dev/null
+++ b/paper/stem_height_px.py
@@ -0,0 +1,59 @@
+import os
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from openalea.maizetrack.local_cache import get_metainfos_ZA17
+from openalea.maizetrack.stem_correction import stem_height_smoothing
+
+from openalea.maizetrack.local_cache import metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.stem_correction import get_median_polyline
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+import json
+
+from sklearn.metrics import r2_score
+from scipy.interpolate import interp1d
+
+# ================================================================
+
+# modulor
+df_tt = pd.read_csv('TT_ZA17.csv')
+
+# modulor
+folder = 'local_cache/cache_ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+#plantids = np.unique([int(f.split('/')[-1][:4]) for f in all_files])
+plantids = [1152,  475,  313, 1303, 1292,  803, 1430, 1424,  309,  958]
+
+dfs = []
+plantid_files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+for f in plantid_files:
+    task = int(f.split('_')[-1].split('.csv')[0])
+    dft = pd.read_csv(f)
+    daydate = f.split('/')[3]
+    tt = df_tt[df_tt['daydate'] == daydate].iloc[0]['ThermalTime']
+    dft['t'] = tt
+    dfs.append(dft)
+df = pd.concat(dfs)
+
+df = df[df['score'] > 0.95]
+df['z_phenomenal'] /= 10
+#df['y'] = 2448 - df['y']
+df_stem = df.sort_values('z_phenomenal', ascending=False).drop_duplicates(['t']).sort_values('t')
+f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['z_phenomenal']))
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+plt.title(plantid)
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+plt.ylabel('Stem height (cm)', fontsize=30)
+
+plt.plot(df_stem['t'], 2448 - df_stem['y'], 'k.', markersize=15)
+plt.legend(prop={'size': 20}, loc='lower right', markerscale=2)
+
+
+
+
+
+
diff --git a/paper/stem_height_validation.py b/paper/stem_height_validation.py
new file mode 100644
index 0000000..0b25e59
--- /dev/null
+++ b/paper/stem_height_validation.py
@@ -0,0 +1,334 @@
+import os
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+from openalea.maizetrack.local_cache import get_metainfos_ZA17
+from openalea.maizetrack.stem_correction import stem_height_smoothing
+
+from openalea.maizetrack.local_cache import metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.stem_correction import get_median_polyline
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+import json
+
+from sklearn.metrics import r2_score
+from scipy.interpolate import interp1d
+
+# ================================================================
+
+# modulor
+df_tt = pd.read_csv('TT_ZA17.csv')
+
+# modulor
+folder = 'local_cache/cache_ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+# available in NASHShare2
+folder2 = 'X:/lepseBinaries/TaggedImages/ARCH2017-03-30/stem_annotation/'
+plantids_anot = [int(f.split('_')[1].split('.')[0]) for f in os.listdir(folder2) if os.path.isfile(folder2 + f)]
+
+
+# ===== plot just one example ================================================================================
+
+plantid = 783
+
+dfs = []
+plantid_files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+for f in plantid_files:
+    task = int(f.split('_')[-1].split('.csv')[0])
+    dft = pd.read_csv(f)
+    daydate = f.split('/')[3]
+    tt = df_tt[df_tt['daydate'] == daydate].iloc[0]['ThermalTime']
+    dft['t'] = tt
+    dfs.append(dft)
+
+df = pd.concat(dfs)
+
+df = df[df['score'] > 0.95]
+df['z_phenomenal'] /= 10
+#df['y'] = 2448 - df['y']
+df_stem = df.sort_values('z_phenomenal', ascending=False).drop_duplicates(['t']).sort_values('t')
+f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['z_phenomenal']))
+
+plt.figure(plantid)
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+plt.title(plantid)
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=30)
+plt.ylabel('Stem height (cm)', fontsize=30)
+
+plt.plot(df_stem['t'], df_stem['z_phenomenal'], 'k.', markersize=15)
+plt.plot(df_stem['t'], [f_smoothing(t) for t in df_stem['t']], 'r-', linewidth=2, label='smoothing')
+
+plt.legend(prop={'size': 20}, loc='lower right', markerscale=2)
+
+
+# =======================================================================================================
+
+res = []
+plt.figure()
+for plantid in [p for p in plantids_anot if p != 452]: # TODO : plantid 452 not in cache !
+
+    print(plantid)
+    metainfos = get_metainfos_ZA17(plantid)
+
+    # ========== basic phenomenal ============================================
+
+    # available in modulor
+    folder_vmsi = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_no_stem_smooth_no_tracking'
+    paths = metainfos_to_paths(metainfos, phm_parameters=(4, 1, 'notop', 4, 100), folder=folder_vmsi)
+
+    metainfos2, paths = check_existence(metainfos, paths)
+
+    vmsi_list = load_plant(metainfos2, paths)
+
+    z_base = {sf: np.median([v.get_stem().get_highest_polyline().polyline[0][2] for v in vmsi_list
+                             if v.metainfo.shooting_frame == sf]) for sf in ['elcom_2_c1_wide', 'elcom_2_c2_wide']}
+
+    stem_phm = {}
+    for v in vmsi_list:
+        heights = [l.info['pm_z_base'] for l in v.get_mature_leafs()]
+        if len(heights) != 0:
+            stem_phm[v.metainfo.timestamp] = max(heights) - z_base[v.metainfo.shooting_frame]
+
+    # ========= median stem on new phenomenal ====================================
+
+    # available in modulor
+    folder_vmsi = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
+    paths = metainfos_to_paths(metainfos, stem_smoothing=True, phm_parameters=(4, 1, 'notop', 4, 100), old=False,
+                               folder=folder_vmsi)
+
+    metainfos2, paths = check_existence(metainfos, paths)
+    vmsi_list = load_plant(metainfos2, paths)
+
+    stem_polylines = [np.array(vmsi.get_stem().get_highest_polyline().polyline) for vmsi in vmsi_list]
+    median_stem = get_median_polyline(polylines=stem_polylines)
+
+    z_base = {sf: np.median([v.get_stem().get_highest_polyline().polyline[0][2] for v in vmsi_list
+                             if v.metainfo.shooting_frame == sf]) for sf in ['elcom_2_c1_wide', 'elcom_2_c2_wide']}
+
+    # ====== stem height annotation =====================================
+
+    t_anot = []
+    y_anot = []
+    z_anot = []
+
+    # available in NASHShare2
+    with open(folder2 + 'stem_{}.json'.format(plantid)) as f:
+        d = json.load(f)
+
+        for name in d.keys():
+
+            task = int(name.split('t')[1].split('_')[0])
+            timestamp = next(m for m in metainfos if m.task == task).timestamp
+            angle = int(name.split('a')[1].split('.')[0])
+            sf = next(m for m in metainfos if m.task == task).shooting_frame
+
+            for shape in d[name]['regions']:
+
+                    y = shape['shape_attributes']['cy']
+                    sf = next(m for m in metainfos if m.task == task).shooting_frame
+                    median_stem_2d = phm3d_to_px2d(median_stem, sf, angle=angle)
+
+                    if y > min(median_stem_2d[:, 1]):
+
+                        f_2dto3d = interp1d(median_stem_2d[:, 1], np.array(median_stem)[:, 2], fill_value="extrapolate")
+                        z = float(f_2dto3d(y))
+                        z -= z_base[sf]
+
+                        t_anot.append(timestamp)
+                        y_anot.append(2448 - y)
+                        z_anot.append(z)
+
+                    else:
+                        print('=============')
+
+    #======== getting full prediction dataset ================================
+
+    dfs = []
+    plantid_files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+    for f in plantid_files:
+        task = int(f.split('_')[-1].split('.csv')[0])
+
+        timestamp = next(m for m in metainfos if m.task == task).timestamp
+        dft = pd.read_csv(f)
+        dft['t'] = timestamp
+
+        dfs.append(dft)
+    df = pd.concat(dfs)
+
+    # ====== 2D ==============================================================
+
+    # data preprocessing, stem extraction
+    df = df[df['score'] > 0.95]
+    df['y'] = 2448 - df['y']
+    df_stem = df.sort_values('y', ascending=False).drop_duplicates(['t']).sort_values('t')
+    f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['y']))
+
+    # ====== 3D ===========================================================
+
+    # METHOD JUST FOR PLOT
+    df_stem['y'] = 2448 - df_stem['y'] # retour a la normale
+    z_row = []
+    for _, row in df_stem.iterrows():
+        sf = next(m for m in metainfos if m.timestamp == row['t']).shooting_frame
+        median_stem_2d = phm3d_to_px2d(median_stem, sf, angle=row['angle'])
+
+        f_2dto3d = interp1d(median_stem_2d[:, 1], np.array(median_stem)[:, 2], fill_value="extrapolate")
+        z = float(f_2dto3d(row['y']))
+        z -= z_base[sf]
+        z_row.append(z)
+
+    df_stem['z'] = z_row
+    df_stem['y'] = 2448 - df_stem['y']  # retour a la normale
+
+
+    # data preprocessing, stem extraction (NORMAL PIPELINE METHOD)
+    #df = df[df['score'] > 0.95]
+    #df_stem_3d = df.sort_values('z_phenomenal', ascending=False).drop_duplicates(['t']).sort_values('t')
+    f_smoothing = stem_height_smoothing(np.array(df_stem['t']), np.array(df_stem['z']))
+    df_stem['z_smooth'] = [f_smoothing(t) for t in df_stem['t']]
+
+    # plt.figure(str(plantid) + '3D')
+    # plt.plot(df['t'], df['z_phenomenal'], 'k*', markersize=2)
+    # plt.plot(df_stem['t'], df_stem['z_phenomenal'], '-g')
+    # plt.plot(df_stem['t'], [f_smoothing(t) for t in df_stem['t']], '-', color='orange')
+
+    # ============ fast plot =================================================
+
+    plt.figure('plantid {} 2D'.format(plantid))
+    plt.plot(t_anot, y_anot, 'g-', label=str(plantid))
+    plt.plot(df_stem['t'], df_stem['y'], 'k-', markersize=2)
+
+    plt.figure('plantid {} 3D'.format(plantid))
+    plt.plot(t_anot, z_anot, 'g-', label=str(plantid))
+    plt.plot(df_stem['t'], df_stem['z'], 'k-', markersize=2)
+    plt.plot(df_stem['t'], df_stem['z_smooth'], 'b-', markersize=2)
+
+
+    # ======= group couples of observation + prediction =======================
+
+    for t, y, z in zip(t_anot, y_anot, z_anot):
+        if t in list(df_stem['t']) and y < 2440:
+            if t in list(stem_phm.keys()):
+                y_pred = float(df_stem[df_stem['t'] == t]['y'])
+                z_pred_3d = float(df_stem[df_stem['t'] == t]['z'])
+                z_pred_3d_smooth = float(df_stem[df_stem['t'] == t]['z_smooth'])
+
+                res.append([y, y_pred, z, z_pred_3d, z_pred_3d_smooth, stem_phm[t]])
+
+            else:
+                print('no data')
+
+res = pd.DataFrame(res, columns=['obs', 'pred', 'obs3d', 'pred3d', 'pred3d_smooth', 'pred_phm'])
+res.to_csv('data/paper/stem_height_validation.csv')
+
+
+res = pd.read_csv('data/paper/stem_height_validation.csv')
+
+# pixel plot
+
+x, y = np.array(res['obs']), np.array(res['pred'])
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+mape = 100 * np.mean(np.abs((x - y) / x))
+r2 = r2_score(x, y)
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+plt.xlim((480, 2500))
+plt.ylim((480, 2500))
+plt.plot([-100, 3000], [-100, 3000], '-', color='grey')
+plt.plot(x, y, 'k.', markersize=4)
+plt.gca().set_aspect('equal', adjustable='box')
+plt.xlabel('observation (px)', fontsize=14)
+plt.ylabel('prediction (px)', fontsize=14)
+plt.title('pixel')
+ax.text(1950, 670, 'R² = {}'.format(round(r2, 3)), fontdict={'size': 12})
+ax.text(1950, 600, 'RMSE = {} px'.format(round(rmse, 2)), fontdict={'size': 12})
+ax.text(1950, 530, 'n = {}'.format(len(x)), fontdict={'size': 12})
+
+# 3D plot
+
+x, y = np.array(res['obs3d']), np.array(res['pred3d'])
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+mape = 100 * np.mean(np.abs((x - y) / x))
+
+for xi, yi in zip(x,y):
+    print( xi, yi, np.abs((xi - yi) / xi) )
+
+r2 = r2_score(x, y)
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+plt.xlim((-60, 2050))
+plt.ylim((-60, 2050))
+plt.plot([-100, 3000], [-100, 3000], '-', color='grey')
+plt.plot(x, y, 'k.', markersize=4)
+plt.gca().set_aspect('equal', adjustable='box')
+plt.xlabel('observation (mm)', fontsize=14)
+plt.ylabel('prediction (mm)', fontsize=14)
+plt.title('3D')
+ax.text(1950, 670, 'R² = {}'.format(round(r2, 3)), fontdict={'size': 12})
+ax.text(1950, 600, 'RMSE = {} mm'.format(round(rmse, 2)), fontdict={'size': 12})
+ax.text(1950, 530, 'n = {}'.format(len(x)), fontdict={'size': 12})
+
+# ===== Final graph, 3D z smooth deep learning vs 3D z phenomenal: ======================================
+
+x, y = np.array(res['obs3d']) / 10, np.array(res['pred3d_smooth']) / 10
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+mape = 100 * np.mean(np.abs((x - y) / x))
+r2 = r2_score(x, y)
+bias = np.mean(x - y)
+
+x2, y2 = np.array(res['obs3d']) / 10, np.array(res['pred_phm']) / 10
+rmse2 = np.sqrt(np.sum((x2 - y2) ** 2) / len(x2))
+mape2 = 100 * np.mean(np.abs((x2 - y2) / x2))
+r2_2 = r2_score(x2, y2)
+bias2 = np.mean(x2 - y2)
+
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.xlim((-6, 205))
+plt.ylim((-6, 205))
+plt.plot([-10, 300], [-10, 300], '-', color='black')
+plt.plot(x2, y2, '.', color='grey', markersize=8, alpha=0.5) # \n RMSE = {} cm, R² = {}'.format(round(rmse2, 1), round(r2_2, 3)))
+a, b = np.polyfit(x, y, 1)
+plt.plot([-10, 3000], a*np.array([-10, 3000]) + b, '--', color='red', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+a, b = np.polyfit(x2, y2, 1)
+plt.plot([-10, 3000], a*np.array([-10, 3000]) + b, '--', color='grey', label=f'y = {a:.{2}f}x {b:+.{2}f}')
+plt.plot(x, y, 'o', color='red', markersize=8, alpha=0.5) # \nstem detection \n RMSE = {} cm, R² = {}'.format(round(rmse, 2), round(r2, 3)))
+plt.gca().set_aspect('equal', adjustable='box')
+plt.xlabel('observation (cm)', fontsize=30)
+plt.ylabel('prediction (cm)', fontsize=30)
+plt.title('Stem height', fontsize=35)
+#ax.text(1950, 670, 'R² = {}'.format(round(r2, 3)), fontdict={'size': 12})
+#ax.text(1950, 600, 'RMSE = {} cm'.format(round(rmse, 2)), fontdict={'size': 12})
+#ax.text(1950, 530, 'n = {}'.format(len(x)), fontdict={'size': 12})
+leg = plt.legend(prop={'size': 22}, loc='upper left', markerscale=1.5)
+for lh in leg.legendHandles:
+    lh._legmarker.set_alpha(1)
+leg.get_texts()[1].set_color('grey')
+leg.get_frame().set_linewidth(0.0)
+leg.get_frame().set_alpha(None)
+
+ax.text(0.05, 0.75, 'n = {} \nBias = {} cm \nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x), round(bias2, 2), round(rmse2, 1), round(mape2, 1), round(r2_2, 3)),
+        transform=ax.transAxes, fontsize=20, verticalalignment='top', color='grey')
+ax.text(0.60, 0.26, 'n = {} \nBias = {} cm \nRMSE = {} cm \nMAPE = {} % \nR² = {}'.format(
+    len(x), round(bias, 2), round(rmse, 2), round(mape, 2), round(r2, 3)),
+        transform=ax.transAxes, fontsize=25, verticalalignment='top', color='black')
+
+
+fig.savefig('paper/results/stem_validation', bbox_inches='tight')
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/paper/tracking_accuracy.py b/paper/tracking_accuracy.py
new file mode 100644
index 0000000..78a9cf3
--- /dev/null
+++ b/paper/tracking_accuracy.py
@@ -0,0 +1,213 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator
+import os
+import time
+import datetime
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+# available in modulor local_benoit
+folder = 'data/tracking/'
+plantids = [int(f.split('.')[0]) for f in os.listdir(folder) if os.path.isfile(folder + f)]
+
+#df_pheno = pd.read_csv('pheno_ZA17.csv')
+#df_pheno[(df_pheno['Pot'] == plantid) & (df_pheno['observationcode'] == 'panicule_visi')].sort_values('resultdate').iloc[0]['resultdate']
+
+# ===== group all plantids in a same dataframe ===============================
+
+df_list = []
+for plantid in plantids:
+
+    dfid = pd.read_csv(folder + '{}.csv'.format(plantid))
+
+    df_list.append(dfid)
+df = pd.concat(df_list)
+
+df = df[~df['rank_annotation'].isin([-1, 0])]
+
+# available in modulor
+tt = pd.read_csv('TT_ZA17.csv')
+t = [tt[tt['timestamp'] < t].sort_values('timestamp').iloc[-1]['ThermalTime'] for t in df['timestamp']]
+df['tt'] = t
+
+#df = df[df['timestamp'] < 1496872800]  # 06-08
+df = df[df['timestamp'] < 1496509351]  # 06-03
+#df = df[df['timestamp'] < 1495310577]  # 05-20
+
+# ===== accuracy per plant ===================================================
+
+#selec = df
+n, accuracy = [], []
+for plantid in df['plantid'].unique():
+
+    selecid = df[df['plantid'] == plantid]
+
+    selec = selecid[(selecid['mature'] == True)]
+    a1 = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec)
+    selec = selecid[(selecid['mature'] == False)]
+    a2 = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec)
+    accuracy.append(a1)
+    print(plantid, round(a1, 3), round(a2, 3), len(selecid))
+    n.append(len(selecid))
+print('=====================================')
+print(min(accuracy), np.mean(accuracy))
+print(min(n), np.mean(n))
+
+# rmse
+df2 = df[~(df['rank_tracking'] == 0)]
+x, y = df2[df2['mature']]['rank_annotation'], df2[df2['mature']]['rank_tracking']
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+print('mature: RMSE = ', rmse)
+x, y = df2[~df2['mature']]['rank_annotation'], df2[~df2['mature']]['rank_tracking']
+rmse = np.sqrt(np.sum((x - y) ** 2) / len(x))
+print('growing: RMSE = ', rmse)
+
+# Brut : 94.9%. -> 06-03 : 97.7%.
+
+# ======= accuracy per plant but no senescence ===========================
+
+n, accuracy = [], []
+for plantid in df['plantid'].unique():
+    dfid = df[df['plantid'] == plantid]
+    ranks = [r for r in sorted(dfid['rank_annotation'].unique())]
+    total, correct = 0, 0
+    for rank in ranks:
+        selec = dfid[(dfid['rank_annotation'] == rank) & (dfid['mature'])]
+        selec = selec.sort_values('timestamp').iloc[:15]
+        total += len(selec)
+        correct += len(selec[selec['rank_annotation'] == selec['rank_tracking']])
+    accuracy.append(correct / total)
+    n.append(total)
+print(min(accuracy), np.mean(accuracy))
+print(min(n), np.mean(n))
+
+
+# ===== accuracy per rank ===================================================
+
+ranks = [r for r in sorted(df['rank_annotation'].unique())]
+res = []
+for m in [True, False]:
+    selec = df[(df['mature'] == m)]
+    n, accuracy = [], []
+    for rank in ranks:
+        selecrank = selec[selec['rank_annotation'] == rank]
+        if not selecrank.empty:
+            accuracy = len(selecrank[selecrank['rank_annotation'] == selecrank['rank_tracking']]) / len(selecrank) * 100
+            n = len(selecrank)
+            res.append([m, rank, accuracy, n])
+res = pd.DataFrame(res, columns=['mature', 'rank', 'accuracy', 'n'])
+
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.xaxis.set_major_locator(MaxNLocator(integer=True)) # int axis
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.xlabel('Leaf rank', fontsize=30)
+plt.ylabel('Tracking accuracy (%)', fontsize=30)
+plt.xlim([0.5, 17])
+plt.ylim([48, 102])
+plt.plot([-100, 100], [100, 100], 'k--')
+
+res = res[res['n'] > 30]
+resm = res[res['mature']]
+plt.plot(resm['rank'], resm['accuracy'], '^-', color='blue', label='Ligulated leaves', markersize=9)
+resg = res[res['mature'] == False]
+plt.plot(resg['rank'], resg['accuracy'], 'o-', color='orange', label='Growing leaves', markersize=9)
+
+leg = plt.legend(prop={'size': 20}, loc='lower right')
+leg.get_frame().set_linewidth(0.0)
+
+fig.savefig('paper/tracking_accuracy', bbox_inches='tight')
+
+
+# ===== accuracy per time period ===========================================
+
+#df['period'] = np.digitize(df['timestamp'], bins=np.linspace(min(df['timestamp']), max(df['timestamp']), 10))
+df['period'], limits = pd.qcut(df['tt'], 15, retbins=True)
+
+res2 = []
+for p in df['period'].unique():
+    dfp = df[df['period'] == p]
+    for m in [True, False]:
+        selec = dfp[(dfp['mature'] == m)]
+        n, accuracy = [], []
+        accuracy = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec) * 100
+        n = len(selec)
+        res2.append([p, m, accuracy, n])
+res2 = pd.DataFrame(res2, columns=['period', 'mature', 'accuracy', 'n'])
+
+fig, ax = plt.subplots()
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.xlabel('Thermal time', fontsize=30)
+plt.ylabel('Tracking accuracy (%)', fontsize=30)
+plt.xlim([min(df['tt']), max(df['tt'])])
+plt.ylim([-4, 104])
+plt.plot([-100, 100], [100, 100], 'k--')
+
+res2 = res2[res2['n'] > 30]
+resm = res2[res2['mature']]
+t = [float(str(x).split(',')[0][1:]) for x in resm['period']]
+plt.plot(t, resm['accuracy'], '^-', color='blue', label='Ligulated leaves', markersize=9)
+resg = res2[res2['mature'] == False]
+t = [float(str(x).split(',')[0][1:]) for x in resg['period']]
+plt.plot(t, resg['accuracy'], 'o-', color='orange', label='Growing leaves', markersize=9)
+plt.legend(prop={'size': 30}, loc='lower left')
+
+# ======== test heatmap ====================================
+
+dfm = df[df['mature']]
+
+dfm['r_interval'] = pd.qcut(dfm['rank_annotation'], 8)
+dfm['t_interval'] = pd.qcut(dfm['timestamp'], 8)
+
+m = np.zeros((len(dfm['r_interval'].unique()), len(dfm['t_interval'].unique())))
+for i, rint in enumerate(dfm['r_interval'].unique()):
+    for j, tint in enumerate(dfm['t_interval'].unique()):
+        selec = dfm[(dfm['r_interval'] == rint) & (dfm['t_interval'] == tint)]
+        if len(selec) > 50:
+            accuracy = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec) * 100
+            m[i, j] = accuracy
+        else:
+            m[i, j] = np.nan
+
+plt.imshow(m, vmin=85, vmax=100, cmap='RdYlGn')
+im_ratio = m.shape[0] / m.shape[1]
+plt.colorbar(fraction=0.046 * im_ratio, pad=0.04, label='Accuracy (%)')
+plt.xlabel('time interval')
+plt.ylabel('rank interval')
+
+
+# ===== growing leaf : jusqu'a quand on peut rembobiner =========================
+
+dfg = df[~df['mature']]
+
+rembobinage = []
+for _, row in dfg.iterrows():
+    selec = dfg[(dfg['plantid'] == row['plantid']) & (dfg['rank_annotation'] == row['rank_annotation'])]
+    selec = selec.sort_values('tt')
+    rembobinage.append(np.max(selec['tt']) - row['tt'])
+dfg['rembobinage'] = rembobinage
+
+dfg['rb_interval'] = pd.qcut(dfg['rembobinage'], 8)
+
+t = [float(str(x).split(',')[0][1:]) for x in dfg['rb_interval']]
+
+x, y = [], []
+for tint in sorted(dfg['rb_interval'].unique()):
+    x.append(float(str(tint).split(',')[0][1:]))
+    selec = dfg[dfg['rb_interval'] == tint]
+    accuracy = len(selec[selec['rank_annotation'] == selec['rank_tracking']]) / len(selec) * 100
+    y.append(accuracy)
+    print(len(selec))
+
+plt.plot(x, y, 'k-*')
+plt.xlabel('thermal time before ligulation')
+plt.ylabel('accuracy')
+
+
+
+
+
+
+
diff --git a/paper/tracking_details_2022.py b/paper/tracking_details_2022.py
new file mode 100644
index 0000000..c71a821
--- /dev/null
+++ b/paper/tracking_details_2022.py
@@ -0,0 +1,377 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import copy
+from matplotlib import colors
+import cv2
+import pandas as pd
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.maizetrack.phenomenal_display import PALETTE, plot_snapshot
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+from openalea.maizetrack.phenomenal_display import plot_vmsi, plot_leaves
+
+
+folder = 'data/tracking/'
+plantids = [int(f.split('.')[0]) for f in os.listdir(folder) if os.path.isfile(folder + f)]
+
+df_list = []
+for plantid in plantids:
+    dfid = pd.read_csv(folder + '{}.csv'.format(plantid))
+    df_list.append(dfid)
+df = pd.concat(df_list)
+df = df[~df['rank_annotation'].isin([-1, 0])]
+# df = df[df['timestamp'] < 1496509351]
+# df = df[df['timestamp'] < 1495922400] # 05-28
+
+# ====== mature tracking, but with different end dates ============================================================*
+
+# selec = df
+n, accuracy = [], []
+for plantid in df['plantid'].unique():
+
+    print(plantid)
+
+    metainfos = get_metainfos_ZA17(plantid)
+
+    paths = metainfos_to_paths(metainfos, folder='local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking')
+    metainfos, paths = check_existence(metainfos, paths)
+    vmsi_list = load_plant(metainfos, paths)
+
+    plant_ref = TrackedPlant.load_and_check(vmsi_list)
+    plant_ref.load_rank_annotation()
+
+    for dmax in ['2017-06-10', '2017-06-03', '2017-05-27', '2017-05-20']:
+
+        plant = copy.deepcopy(plant_ref)
+
+        plant.snapshots = [s for s in plant.snapshots if s.metainfo.daydate < dmax]
+
+        plant.align_mature(direction=1, gap=12.365, w_h=0.03, w_l=0.002, gap_extremity_factor=0.2, n_previous=500,
+                           rank_attribution=True)
+        snapshots = plant.snapshots
+
+        df_res = plant.get_dataframe(load_anot=False)
+        df_res = df_res[df_res['mature']]
+        print(dmax, round(len(df_res[df_res['rank_tracking'] == df_res['rank_annotation']]) / len(df_res), 2))
+        df_res.to_csv('data/tracking/tests_2022/{}_{}.csv'.format(plantid, dmax))
+
+# ===== growing tracking, but with/without ground-truth mature tracking ==============================================
+
+# selec = df
+n, accuracy = [], []
+for plantid in df['plantid'].unique():
+
+    print(plantid)
+
+    metainfos = get_metainfos_ZA17(plantid)
+
+    paths = metainfos_to_paths(metainfos, folder='local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking')
+    metainfos, paths = check_existence(metainfos, paths)
+    vmsi_list = load_plant(metainfos, paths)
+
+    plant_ref = TrackedPlant.load_and_check(vmsi_list)
+    plant_ref.load_rank_annotation()
+    plant_ref.snapshots = [s for s in plant_ref.snapshots if s.metainfo.daydate < '2017-06-03']
+
+    for method in ['normal', 'ground-truth']:
+
+        plant = copy.deepcopy(plant_ref)
+
+        # a) normal tracking
+        if method == 'normal':
+            plant.align_mature(direction=1, gap=12.365, w_h=0.03, w_l=0.002, gap_extremity_factor=0.2, n_previous=500,
+                               rank_attribution=True)
+            plant.align_growing()
+
+        # b) tracking starting from ground-truth
+        elif method == 'ground-truth':
+            for s in plant.snapshots:
+                rmax = 99
+                s.order = [-1 if i not in s.rank_annotation else s.rank_annotation.index(i) for i in range(rmax + 1)]
+                s.order = [i if (i != -1 and s.leaves[i].info['pm_label'] == 'mature_leaf') else -1 for i in s.order]
+            rmax2 = [all([s.order[i] == -1 for s in plant.snapshots]) for i in range(rmax)].index(True)
+            for s in plant.snapshots:
+                s.order = s.order[:rmax2]
+            plant.align_growing()
+
+        df_res = plant.get_dataframe(load_anot=False)
+        df_res = df_res[~df_res['mature']]
+        print(method, round(len(df_res[df_res['rank_tracking'] == df_res['rank_annotation']]) / len(df_res), 2))
+        df_res.to_csv('data/tracking/tests_2022/growing_{}_{}.csv'.format(plantid, method))
+
+# ===============================================================================================================
+# ===============================================================================================================
+# ===============================================================================================================
+# ===============================================================================================================
+# ===============================================================================================================
+# ===============================================================================================================
+
+df_res = []
+for plantid in df['plantid'].unique():
+    for dmax in ['2017-06-10', '2017-06-03', '2017-05-27', '2017-05-20']:
+        dfi = pd.read_csv('data/tracking/tests_2022/{}_{}.csv'.format(plantid, dmax))
+        dfi['dmax'] = dmax
+        dfi['method'] = None
+        dfi['mature'] = True
+        df_res.append(dfi)
+    for method in ['normal', 'ground-truth']:
+        dfi = pd.read_csv('data/tracking/tests_2022/growing_{}_{}.csv'.format(plantid, method))
+        dfi['dmax'] = None
+        dfi['method'] = method
+        dfi['mature'] = False
+        df_res.append(dfi)
+df_res = pd.concat(df_res)
+df_res = df_res[~df_res['rank_annotation'].isin([-1, 0])]
+df_res['i'] = np.arange(len(df_res))
+
+# remove growing leaves that didn't reach mature
+for plantid in df_res['plantid'].unique():
+    # mature ground-truth (max rank)
+    selec_m = df_res[(df_res['dmax'] == '2017-06-03') & (df_res['plantid'] == plantid)]
+    # growing ground-truth above
+    selec_g = df_res[(df_res['plantid'] == plantid) &
+                     (df_res['rank_annotation'] > max(selec_m['rank_annotation']))]
+    # remove above
+    df_res = df_res[~df_res['i'].isin(list(selec_g['i']))]
+
+ranks = [r for r in sorted(df_res['rank_annotation'].unique())]
+res = []
+
+# ===== A) MATURE LEAVES =======
+
+dmax = '2017-06-03'  # pipeline simulation end-date
+
+selec = df_res[(df_res['dmax'] == dmax)]
+for rank in ranks:
+    selecrank = selec[selec['rank_annotation'] == rank]
+    if not selecrank.empty:
+        acc_list = []
+        n_plants = 0
+        for plantid in selecrank['plantid'].unique():
+            s = selecrank[selecrank['plantid'] == plantid]
+            if not s.empty:
+                n_plants += 1
+                accuracy = len(s[s['rank_annotation'] == s['rank_tracking']]) / len(s) * 100
+                acc_list.append(accuracy)
+
+        # compute bootstrap confidence interval
+        n_bootstrap = 5000
+        conf = 0.95
+        bootstrap = sorted([np.mean(np.random.choice(acc_list, len(acc_list))) for k in range(n_bootstrap)])
+        err1, err2 = bootstrap[int(n_bootstrap * ((1 - conf)/2))], bootstrap[int(n_bootstrap * (conf + (1 - conf)/2))]
+
+        res.append([True, 'mature', rank, np.mean(acc_list), err1, err2, len(selecrank), n_plants])
+
+
+# ===== a) remove leaves too old
+to_keep = []
+for plantid in selec['plantid'].unique():
+    s = selec[selec['plantid'] == plantid]
+    for r in s['rank_annotation'].unique():
+        s2 = s[s['rank_annotation'] == r]
+        to_keep = to_keep + list(s2.sort_values('timestamp').iloc[:10]['i'])
+selec = selec[selec['i'].isin(to_keep)]
+
+# ===== b) remove timestamps with ear
+s = selec[selec['plantid'] == 461]
+s1 = s[s['timestamp'].isin(sorted(s['timestamp'].unique())[-2:])]
+s = selec[selec['plantid'] == 1391]
+s2 = s[s['timestamp'].isin(sorted(s['timestamp'].unique())[-2:])]
+s = selec[selec['plantid'] == 905]
+s3 = s[s['timestamp'].isin(sorted(s['timestamp'].unique())[-1:])]
+s = selec[selec['plantid'] == 832]
+s4 = s[s['timestamp'].isin(sorted(s['timestamp'].unique())[-1:])]
+selec = selec[~selec['i'].isin(list(s1['i']) + list(s2['i']) + list(s3['i']) + list(s4['i']))]
+
+# ===== c) remove timestamps with last ligul
+# to_keep = []
+# for plantid in selec['plantid'].unique():
+#     s = selec[selec['plantid'] == plantid]
+#     s2 = s[s['rank_annotation'] == max(s['rank_annotation'])]
+#     to_keep += list(s[s['timestamp'] < min(s2['timestamp'])]['i'])
+# selec = selec[selec['i'].isin(to_keep)]
+
+for rank in ranks:
+    selecrank = selec[selec['rank_annotation'] == rank]
+    if not selecrank.empty:
+        acc_list = []
+        n_plants = 0
+        for plantid in selecrank['plantid'].unique():
+            s = selecrank[selecrank['plantid'] == plantid]
+            if not s.empty:
+                n_plants += 1
+                accuracy = len(s[s['rank_annotation'] == s['rank_tracking']]) / len(s) * 100
+                acc_list.append(accuracy)
+
+        # compute bootstrap confidence interval
+        n_bootstrap = 5000
+        conf = 0.95
+        bootstrap = sorted([np.mean(np.random.choice(acc_list, len(acc_list))) for k in range(n_bootstrap)])
+        err1, err2 = bootstrap[int(n_bootstrap * ((1 - conf)/2))], bootstrap[int(n_bootstrap * (conf + (1 - conf)/2))]
+
+        res.append([True, 'mature2', rank, np.mean(acc_list), err1, err2, len(selecrank), n_plants])
+
+# ===== B) GROWING LEAVES =======
+
+for (criteria, method) in zip(['normal', 'ground-truth'], ['growing', 'growing2']):
+    selec = df_res[(df_res['method'] == criteria)]
+    for rank in df_res[(df_res['dmax'] == '2017-06-03')]['rank_annotation'].unique():
+        selecrank = selec[selec['rank_annotation'] == rank]
+        if not selecrank.empty:
+            acc_list = []
+            n_plants = 0
+            for plantid in selecrank['plantid'].unique():
+                s = selecrank[selecrank['plantid'] == plantid]
+                if not s.empty:
+                    n_plants += 1
+                    accuracy = len(s[s['rank_annotation'] == s['rank_tracking']]) / len(s) * 100
+                    acc_list.append(accuracy)
+
+            # compute bootstrap confidence interval
+            n_bootstrap = 5000
+            conf = 0.95
+            bootstrap = sorted([np.mean(np.random.choice(acc_list, len(acc_list))) for k in range(n_bootstrap)])
+            err1, err2 = bootstrap[int(n_bootstrap * ((1 - conf) / 2))], bootstrap[
+                int(n_bootstrap * (conf + (1 - conf) / 2))]
+
+            res.append([False, method, rank, np.mean(acc_list), err1, err2, len(selecrank), n_plants])
+
+res = pd.DataFrame(res, columns=['mature', 'option', 'rank', 'accuracy', 'err1', 'err2', 'n', 'n_plants'])
+
+# ===== C) PLOT =====
+
+from matplotlib.ticker import MaxNLocator
+fig, ax = plt.subplots(figsize=(10, 10), dpi=100)
+ax.xaxis.set_major_locator(MaxNLocator(integer=True)) # int axis
+ax.tick_params(axis='both', which='major', labelsize=20)
+plt.xlabel('Leaf rank', fontsize=30)
+plt.ylabel('Rank assignment accuracy (%)', fontsize=30)
+plt.xlim([0.5, 14.5])
+plt.ylim([65, 101])
+plt.plot([-100, 100], [100, 100], 'k--')
+
+res['pc'] = [row['n'] / sum(res[res['option'] == row['option']]['n']) for _, row in res.iterrows()]
+
+res_selec = res[res['n_plants'] > 15]
+#res_selec = res[res['n'] > 40]
+# res_selec = res[res['pc'] > 0.0000012]
+# res_selec = res_selec[res_selec['rank'] <= 14]
+
+label_dic = {'mature': 'Ligulated leaves',
+             'mature2': None,
+             'growing': 'Growing leaves',
+             'growing2': None}
+
+for option in label_dic.keys():
+
+    s = res_selec[res_selec['option'] == option]
+
+    color = 'blue' if 'mature' in option else 'orange'
+    symbol = '^' if 'mature' in option else 'o'
+    #dx = 0.02 if s.iloc[0]['mature'] else -0.02
+
+    linewidth = 3 if option in ['mature', 'growing'] else 1.2
+    markersize = 10 if option in ['mature', 'growing'] else 5
+    linestyle = '-' if option in ['mature', 'growing'] else '--'
+    linestyle = ':' if option == 'mature2' else linestyle
+
+    # color, linestyle = ('green', '-.') if option == 'old removed' else (color, linestyle)
+    # color, linestyle = ('purple', ':') if option == 'old removed + no ear' else (color, linestyle)
+    # color, linestyle = ('black', ':') if option == 'old removed + no ear + no last ligul' else (color, linestyle)
+    # color = 'red' if option == 'ground-truth' else color
+
+    plt.plot(s['rank'], s['accuracy'], symbol, color=color,
+             label=label_dic[option],
+             markersize=markersize, linewidth=linewidth)
+    plt.plot(s['rank'], s['accuracy'], linestyle + symbol, color=color,
+             markersize=markersize, linewidth=linewidth)
+
+    if option in ['mature', 'growing']:
+        plt.fill_between(s['rank'], s['err1'], s['err2'], color=color, alpha=0.15 if color == 'blue' else 0.25)
+
+leg = plt.legend(prop={'size': 30}, loc='lower left')
+leg.get_frame().set_linewidth(0.0)
+# remove "white non transparent rectangle" around legend
+leg.get_frame().set_alpha(None)
+leg.get_frame().set_facecolor((0, 0, 0, 0))
+
+
+# ===== visu some matrixes to understand =========================================================================
+
+# 1) plot errors
+
+res2 = []
+
+#for plantid in df['plantid'].unique()[:1]:
+for method in ['normal', 'ground-truth']:
+    dfm = pd.read_csv('data/tracking/tests_2022/growing_{}_{}.csv'.format(plantid, method))
+    for rank in dfm['rank_annotation'].unique():
+        selecrank = dfm[dfm['rank_annotation'] == rank]
+        accuracy = len(selecrank[selecrank['rank_annotation'] == selecrank['rank_tracking']]) / len(selecrank) * 100
+        res2.append([False, method, rank, accuracy, len(selecrank)])
+
+# plt.figure()
+# res2 = pd.DataFrame(res2, columns=['mature', 'option', 'rank', 'accuracy', 'n']).sort_values('rank')
+# s1 = res2[res2['option'] == 'normal']
+# plt.plot(s1['rank'], s1['accuracy'], 'r-*')
+# s2 = res2[res2['option'] == 'ground-truth']
+# plt.plot(s2['rank'] + 0.1, s2['accuracy'], 'g-*')
+
+# 2) plot matrix
+
+for plantid in plantids:
+
+    #df1 = pd.read_csv('data/tracking/tests_2022/growing_{}_{}.csv'.format(plantid, method))
+    df2 = pd.read_csv('data/tracking/tests_2022/{}_{}.csv'.format(plantid, '2017-06-03'))
+    #dfm = pd.concat((df1, df2))
+    dfm = df2
+
+    T = len(dfm['timestamp'].unique())
+    R = max(dfm['rank_tracking'].unique())
+    mat = np.zeros((T, R)) * np.NAN
+    for t_index, t in enumerate(dfm['timestamp'].unique()):
+        s = dfm[dfm['timestamp'] == t]
+        for _, row in s.iterrows():
+            mat[t_index, row['rank_tracking'] - 1] = row['rank_annotation'] - 1
+
+    # remove empty columns
+    mat = mat[:, ~np.isnan(mat).all(axis=0)]
+
+    fig, ax = plt.subplots()
+    fig.canvas.set_window_title(str(plantid))
+
+    #ax.set_xticks(np.arange(R) - 0.5, minor=True)
+    # ax.set_yticks(np.arange(T + 1) - 0.5, minor=True)
+    # ax.grid(which='minor', color='white', linewidth=12)
+
+    # # start axis at 1
+    # plt.xticks(np.arange(R), np.arange(R) + 1)
+    # plt.yticks(np.arange(T), np.arange(T) + 1)
+    #
+    # ax.set_xlabel('Leaf rank', fontsize=30)
+    # ax.xaxis.tick_top()
+    # ax.xaxis.set_label_position('top')
+    # ax.set_ylabel('Time', fontsize=30)
+    # plt.locator_params(nbins=4)
+    # ax.tick_params(axis='both', which='major', labelsize=25)  # axis number size
+
+    rgb = np.array(PALETTE) / 255.
+    rgb = np.concatenate((np.array([[0., 0., 0.]]), rgb))
+    cmap = colors.ListedColormap(rgb, "")
+    val = [k - 1.5 for k in range(50)]
+    norm = colors.BoundaryNorm(val, len(val)-1)
+    plt.imshow(mat, interpolation='nearest', cmap=cmap, norm=norm)
+
+    plt.title(plantid)
+
+    # for t_index, t in enumerate(dfm['timestamp'].unique()):
+    #     s = dfm[dfm['timestamp'] == t]
+    #     for _, row in s.iterrows():
+    #         if row['mature']:
+    #             plt.plot(row['rank_tracking'], t_index, 'k.')
+    #             #mat[t_index, row['rank_tracking']] = row['rank_annotation'] - 1
+    #
+    #
diff --git a/paper/tracking_matrix.py b/paper/tracking_matrix.py
new file mode 100644
index 0000000..cb7ed66
--- /dev/null
+++ b/paper/tracking_matrix.py
@@ -0,0 +1,139 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import copy
+from matplotlib import colors
+import cv2
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.maizetrack.phenomenal_display import PALETTE, plot_snapshot
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+
+def rgb_and_polylines2(snapshot, angle, ranks=None, ds=1):
+
+    # rgb image
+    img = snapshot.image[angle]
+
+    # adding polylines
+    polylines_px = [phm3d_to_px2d(leaf.real_pl, snapshot.metainfo.shooting_frame, angle) for leaf in snapshot.leaves]
+    matures = [leaf.info['pm_label'] == 'mature_leaf' for leaf in snapshot.leaves]
+
+    if ranks is None:
+        ranks = snapshot.rank_annotation
+
+    # plot leaves
+    for i, (pl, c, mature) in enumerate(zip(polylines_px, ranks, matures)):
+
+        if c == -2:
+            col = [0, 0, 0]
+        else:
+            col = [int(x) for x in PALETTE[c]]
+
+        # unknown rank vs known rank
+        if c == -1:
+            leaf_border_col = (0, 0, 0)
+        else:
+            leaf_border_col = (255, 255, 255)
+
+        #img = cv2.polylines(np.float32(img), [pl.astype(int).reshape((-1, 1, 2))], False, leaf_border_col, 10 * ds)
+        img = cv2.polylines(np.float32(img), [pl.astype(int).reshape((-1, 1, 2))], False, col, 7 * ds)
+
+    return img, polylines_px
+
+
+# available in modulor
+folder = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+plantids = [313, 316, 329, 330, 336, 348, 424, 439, 461, 474, 794, 832, 905, 907, 915, 925, 931, 940, 959, 1270,
+            1276, 1283, 1284, 1301, 1316, 1383, 1391, 1421, 1424, 1434]
+
+
+plantid = 313
+
+metainfos = get_metainfos_ZA17(plantid)
+paths = metainfos_to_paths(metainfos, folder=folder)
+metainfos, paths = check_existence(metainfos, paths)
+vmsi_list = load_plant(metainfos, paths)
+
+plant_ref = TrackedPlant.load_and_check(vmsi_list)
+plant_ref.load_rank_annotation()
+
+plant_aligned = copy.deepcopy(plant_ref)
+plant_aligned.align_mature(direction=1, gap=12.365, w_h=0.03, w_l=0.002, gap_extremity_factor=0.2, n_previous=500,
+                           rank_attribution=True)
+plant_aligned.align_growing()
+
+plant = plant_aligned
+only_mature = False
+
+snapshots = plant.snapshots[0::3]
+
+
+from skimage import io
+plant.load_images(60)
+plant.simplify_polylines()
+s = copy.deepcopy(snapshots[-2])
+s.leaves = [l for l in s.leaves if l.info['pm_label'] == 'mature_leaf']
+img, _ = rgb_and_polylines2(s, angle=60, ranks=[5, 6, 7, 8, 9, 10, 18, 11, 12, 13], ds=3)
+io.imsave('test.png', img)
+
+s = copy.deepcopy(snapshots[-8])
+s.leaves = [l for l in s.leaves if l.info['pm_label'] == 'growing_leaf']
+img, _ = rgb_and_polylines2(s, angle=60, ranks=[10, 12, 11, 13], ds=3)
+io.imsave('test2.png', img)
+
+#plot_snapshot(s, ranks=[5, 6, 7, 8, 9, 10, 18, 11, 12, 13], stem=False)
+
+
+T = len(snapshots)
+R = len(snapshots[0].order)
+mat = np.zeros((T, R)) * np.NAN
+for t, snapshot in enumerate(snapshots):
+    for r, index in enumerate(snapshot.order):
+        if index != -1:
+            if snapshot.leaves[index].info['pm_label'] == 'mature_leaf' or only_mature == False:
+                mat[t, r] = snapshot.leaves[index].rank_annotation
+
+# remove empty columns
+mat = mat[:, ~np.isnan(mat).all(axis=0)]
+
+fig, ax = plt.subplots()
+fig.canvas.set_window_title(str(plantid))
+
+#ax.set_xticks(np.arange(R) - 0.5, minor=True)
+ax.set_yticks(np.arange(T + 1) - 0.5, minor=True)
+ax.grid(which='minor', color='white', linewidth=12)
+
+# start axis at 1
+plt.xticks(np.arange(R), np.arange(R) + 1)
+plt.yticks(np.arange(T), np.arange(T) + 1)
+
+ax.set_xlabel('Leaf rank', fontsize=30)
+ax.xaxis.tick_top()
+ax.xaxis.set_label_position('top')
+ax.set_ylabel('Time', fontsize=30)
+plt.locator_params(nbins=4)
+ax.tick_params(axis='both', which='major', labelsize=25)  # axis number size
+
+
+rgb = np.array(PALETTE) / 255.
+rgb = np.concatenate((np.array([[0., 0., 0.]]), rgb))
+cmap = colors.ListedColormap(rgb, "")
+val = [k - 1.5 for k in range(50)]
+norm = colors.BoundaryNorm(val, len(val)-1)
+plt.imshow(mat, interpolation='nearest', cmap=cmap, norm=norm)
+
+
+
+# ==================================================================================================================
+
+
+
+
+
+
+
+
diff --git a/paper/video.py b/paper/video.py
new file mode 100644
index 0000000..e3b97d8
--- /dev/null
+++ b/paper/video.py
@@ -0,0 +1,113 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import copy
+from matplotlib import colors
+import cv2
+
+from openalea.maizetrack.local_cache import get_metainfos_ZA17, metainfos_to_paths, check_existence, load_plant
+from openalea.maizetrack.trackedPlant import TrackedPlant
+from openalea.maizetrack.phenomenal_display import PALETTE, plot_snapshot
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+
+def rgb_and_polylines2(snapshot, angle, ranks=None, ds=1):
+
+    # rgb image
+    img = snapshot.image[angle]
+
+    # adding polylines
+    polylines_px = [phm3d_to_px2d(leaf.real_pl, snapshot.metainfo.shooting_frame, angle) for leaf in snapshot.leaves]
+    matures = [leaf.info['pm_label'] == 'mature_leaf' for leaf in snapshot.leaves]
+
+    if ranks is None:
+        #ranks = snapshot.rank_annotation
+        ranks = snapshot.get_ranks()
+
+    # plot leaves
+    for i, (pl, c, mature) in enumerate(zip(polylines_px, ranks, matures)):
+
+        if c == -1:
+            col = [255, 255, 255]
+        else:
+            col = [int(x) for x in PALETTE[c]]
+
+        # unknown rank vs known rank
+        if c == -1:
+            leaf_border_col = (0, 0, 0)
+        else:
+            leaf_border_col = (255, 255, 255)
+
+        #img = cv2.polylines(np.float32(img), [pl.astype(int).reshape((-1, 1, 2))], False, leaf_border_col, 10 * ds)
+        img = cv2.polylines(np.float32(img), [pl.astype(int).reshape((-1, 1, 2))], False, col, 7 * ds)
+
+    # write date
+    id = str(int(snapshot.metainfo.plant[:4]))
+    text = 'plantid {} / task {} ({})'.format(id, snapshot.metainfo.task, snapshot.metainfo.daydate)
+    cv2.rectangle(img, (0, 2250), (2048, 2448), (0, 0, 0), -1)
+    cv2.putText(img, text, (50, 2380), cv2.FONT_HERSHEY_SIMPLEX, 3, (255, 255, 255), 9, cv2.LINE_AA)
+
+    return img, polylines_px
+
+
+# ===================================================================================================
+
+# available in modulor
+folder = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
+all_files = [folder + '/' + rep + '/' + f for rep in os.listdir(folder) for f in os.listdir(folder + '/' + rep)]
+
+plantids = [313, 316, 329, 330, 336, 348, 424, 439, 461, 474, 794, 832, 905, 907, 915, 925, 931, 940, 959, 1270,
+            1276, 1283, 1284, 1301, 1316, 1383, 1391, 1421, 1424, 1434]
+
+for plantid in plantids:
+
+    metainfos = get_metainfos_ZA17(plantid)
+
+    metainfos = [m for m in metainfos if m.daydate < '2017-06-08']
+
+    paths = metainfos_to_paths(metainfos, folder=folder)
+    metainfos, paths = check_existence(metainfos, paths)
+    vmsi_list = load_plant(metainfos, paths)
+
+    plant = TrackedPlant.load_and_check(vmsi_list)
+
+    plant.align_mature(direction=1, gap=12.365, w_h=0.03, w_l=0.004, gap_extremity_factor=0.2, n_previous=500)
+    plant.align_growing()
+
+    plant.load_images(60)
+
+
+    imgs = []
+    for s in plant.snapshots:
+        img, _ = rgb_and_polylines2(s, 60)
+        #plt.imshow(img.astype(np.uint8))
+        imgs.append(img.astype(np.uint8))
+
+    # ========================================================
+
+    width = 2048
+    height = 2448
+    channel = 3
+    fps = 3
+
+    fourcc = cv2.VideoWriter_fourcc(*'MP42')
+    video = cv2.VideoWriter('data/videos/mp4/{}.mp4'.format(plantid), fourcc, float(fps), (width, height))
+    for img in imgs:
+        img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)
+        video.write(img_bgr)
+
+    video.release()
+
+    # =============================================================
+
+    from PIL import Image
+
+    imgs_gif = imgs.copy()
+    imgs_gif = [Image.fromarray(np.uint8(img)) for img in imgs_gif]
+    fps = 3
+    imgs_gif[0].save('data/videos/gif/{}.gif'.format(plantid),
+                     save_all=True,
+                     append_images=imgs_gif[1:],
+                     optimize=True,
+                     duration=1000 / fps,
+                     loop=0)
\ No newline at end of file
diff --git a/paper/visu_ZA17.py b/paper/visu_ZA17.py
new file mode 100644
index 0000000..93c0aff
--- /dev/null
+++ b/paper/visu_ZA17.py
@@ -0,0 +1,196 @@
+import pandas as pd
+import numpy as np
+import os
+from skimage import io
+import matplotlib.pyplot as plt
+from openalea.maizetrack.utils import get_rgb
+from openalea.maizetrack.local_cache import get_metainfos_ZA17
+from openalea.maizetrack.stem_correction import smoothing_function, stem_height_smoothing
+from scipy.interpolate import interp1d
+import time
+import datetime
+from matplotlib.ticker import MaxNLocator
+
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+colors = {'WW': 'dodgerblue', 'WD': 'orangered', 'error': 'gainsboro'}
+
+# thermal time : available in modulor
+df_tt = pd.read_csv('data/TT_ZA17.csv')
+
+# collar detection data
+# available in modulor
+folder_dl = 'local_cache/cache_ZA17/collars_voxel4_tol1_notop_vis4_minpix100'
+all_files_dl = [folder_dl + '/' + rep + '/' + f for rep in os.listdir(folder_dl) for f in os.listdir(folder_dl + '/' + rep)]
+
+# ZA17 GxE data
+# available in modulor
+pheno_ZA17 = pd.read_csv('data/pheno_ZA17.csv')
+pheno_ZA17['timestamp'] = pheno_ZA17['resultdate'].map(date_to_timestamp)
+gxe = pheno_ZA17[['Pot', 'Variety_ID', 'Scenario', 'Rep', 'Sampling']].drop_duplicates()
+gxe = gxe[gxe['Sampling'] == 'Rep']
+gxe = gxe[gxe['Variety_ID'] != 'SolNu']
+
+
+# ===== stem =========================================================================================================
+
+res_stem = []
+for genotype in gxe['Variety_ID'].unique():
+
+    selec = gxe[gxe['Variety_ID'] == genotype]
+
+    for _, row in selec.iterrows():
+
+        plantid_files = [f for f in all_files_dl if int(f.split('/')[-1][:4]) == row['Pot']]
+        if len(plantid_files) > 45:
+
+            dfid = []
+            for f in plantid_files:
+                task = int(f.split('_')[-1].split('.csv')[0])
+                #timestamp = next(m for m in metainfos if m.task == task).timestamp
+                daydate = f.split('__')[-2]
+                tt = df_tt[df_tt['daydate'] == daydate].iloc[0]['ThermalTime']
+                dft = pd.read_csv(f)
+                dft['tt'] = tt
+                dfid.append(dft)
+            dfid = pd.concat(dfid)
+
+            # remove > 06-08
+            dfid = dfid[dfid['tt'] < 90.5]
+
+            #  remove big time gaps
+            times = list(dfid.sort_values('tt')['tt'].unique())
+            for i in range(1, len(times)):
+                if (times[i] - times[i - 1]) > 5 * np.mean(np.diff(times)):
+                    dfid = dfid[dfid['tt'] < times[i]]
+
+            # data preprocessing, stem extraction
+            dfid = dfid[dfid['score'] > 0.95]
+            df_stem = dfid.sort_values('z_phenomenal', ascending=False).drop_duplicates(['tt']).sort_values('tt')
+            f_smoothing = stem_height_smoothing(np.array(df_stem['tt']), np.array(df_stem['z_phenomenal']))
+            df_stem['z_smooth'] = [f_smoothing(t) for t in df_stem['tt']]
+
+            for _, row_stem in df_stem.iterrows():
+                res_stem.append([row['Variety_ID'], row['Pot'], row['Scenario'], row_stem['z_smooth'], row_stem['tt']])
+
+
+res_stem = pd.DataFrame(res_stem, columns=['genotype', 'plantid', 'scenario', 'z', 't'])
+
+res_gxe = res_stem.copy()
+res_gxe['z'] /= 10
+#res_gxe = res_gxe[res_gxe['t'] < np.min(res_gxe.groupby(['scenario', 'genotype']).max('t')['t'])] # same tmax for each scenario
+res_gxe = res_gxe.groupby(['genotype', 'scenario', 'plantid', 't']).agg(['median', 'count']).reset_index()
+#res_gxe = res_gxe[res_gxe['z']['count'] > 2] # at least 3 rep
+#res_gxe = res_gxe[res_gxe['t'] <= min(res_gxe.groupby('scenario').max('t')['t'])]  # same end for each scenario
+
+h_cabin = 300
+# res_gxe = res_gxe[res_gxe['z']['median'] < h_cabin] # cabin limit
+
+fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+fig.canvas.set_window_title(genotype)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+# for scenario in ['WW', 'WD']:
+#     for genotype in res_gxe['genotype'].unique():
+for plantid in res_gxe['plantid'].unique():
+    selec = res_gxe[(res_gxe['plantid'] == plantid)]
+    scenario = selec.iloc[0]['scenario'][0]
+    plt.plot(selec['t'], selec['z']['median'], '-', color=colors[scenario], linewidth=0.2)
+        #plt.plot(selec['t'], selec['z']['median'], '.', color=colors[scenario], label=scenario, markersize=12)
+plt.plot([-100, 1000], [h_cabin, h_cabin], 'k--', label='cabin limit')
+plt.xlabel('Thermal time $(day_{20°C})$', fontsize=25)
+plt.ylabel('Stem height (cm)', fontsize=25)
+plt.title('Stem growth', fontsize=35)
+plt.xlim([8, 93])
+plt.ylim([-2, 202])
+# lgd = plt.legend(prop={'size': 20}, loc='upper left', markerscale=2, handletextpad=-0.5)
+# lgd.get_frame().set_linewidth(0.0)
+
+fig.savefig('paper/stem_360'.format(genotype), bbox_inches='tight')
+# plt.close('all')
+
+fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+fig.canvas.set_window_title(genotype)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+dz = 4.5
+selec = res_gxe[res_gxe['t'] == res_gxe['t'].unique()[24]]
+print(selec.iloc[0]['t'])
+for scenario in ['WW', 'WD']:
+    s = selec[selec['scenario'] == scenario]
+    data = s['z']['median']
+    bins = int(np.ceil((data.max() - data.min()) / dz))
+    plt.hist(s['z']['median'], bins=bins, color=colors[scenario], alpha=0.8)
+plt.xlabel('Stem height (cm)', fontsize=30)
+
+
+fig.savefig('paper/stem_hist_360'.format(genotype), bbox_inches='tight')
+
+# ===== profile ======================================================================================================
+
+res_profile_l = []
+for genotype in gxe['Variety_ID'].unique():
+    print(genotype)
+    selec = gxe[gxe['Variety_ID'] == genotype]
+    for _, row in selec.iterrows():
+
+        # available in modulor local_benoit
+        file = 'data/tracking/results_tracking_360/{}.csv'.format(row['Pot'])
+        if os.path.isfile(file):
+
+            dfid = pd.read_csv(file)
+
+            dfid_mature = dfid[dfid['mature']]
+            ranks = sorted([r for r in dfid_mature['rank_tracking'].unique() if r != 0])
+            profile = []
+            for r in ranks:
+                dfr = dfid_mature[dfid_mature['rank_tracking'] == r]
+                dfr = dfr.sort_values('timestamp')[:15]
+
+                res_profile_l.append([row['Variety_ID'], row['Pot'], row['Scenario'], np.median(dfr['l_extended']), r])
+
+
+res_profile_l = pd.DataFrame(res_profile_l, columns=['genotype', 'plantid', 'scenario', 'l', 'r'])
+res_profile_l['l'] /= 10
+
+fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+fig.canvas.set_window_title('profile')
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+
+for scenario in ['WW', 'WD']:
+    for r in res_profile_l['r'].unique():
+        s = res_profile_l[(res_profile_l['scenario'] == scenario) & (res_profile_l['r'] == r)]
+        if len(s) > 10:
+            dr = 0.07
+            if scenario == 'WW':
+                dr *= -1
+            plt.errorbar(r + dr, np.median(s['l']), yerr=np.std(s['l']), marker='.', color=colors[scenario], capsize=3)
+
+plt.xlabel('Rank', fontsize=25)
+plt.ylabel('Ligulated leaf length (cm)', fontsize=25)
+plt.title('Leaf length profile', fontsize=35)
+
+fig.savefig('paper/profile_360'.format(genotype), bbox_inches='tight')
+
+
+# hist
+
+fig, ax = plt.subplots(figsize=(10, 6), dpi=150)
+ax.tick_params(axis='both', which='major', labelsize=20)  # axis number size
+#dz = 5
+selec = res_profile_l[res_profile_l['r'] == 10]
+for scenario in ['WW', 'WD']:
+    s = selec[selec['scenario'] == scenario]
+    data = s['l']
+    #bins = int(np.ceil((data.max() - data.min()) / dz))
+    bins = 30
+    plt.hist(data, bins=bins, range=[0, 150], color=colors[scenario], alpha=0.8)
+plt.xlabel('Leaf length (cm)', fontsize=30)
+
+fig.savefig('paper/profile_hist_360'.format(genotype), bbox_inches='tight')
+
+
+
+
+
diff --git a/scripts/dataset_ZA22.py b/scripts/dataset_ZA22.py
new file mode 100644
index 0000000..94327f8
--- /dev/null
+++ b/scripts/dataset_ZA22.py
@@ -0,0 +1,143 @@
+import pandas as pd
+import numpy as np
+import shutil
+import os
+
+# =================================================================================================
+
+
+
+# ===============================================================================================
+
+df_db = pd.read_csv('data/copy_from_database/images_ZA22.csv')
+
+img_files = [f for f in os.listdir('data/rgb_insertion_annotation/dataset_2022/images/') if 'ZA22' in f]
+sk_files = [f for f in os.listdir('data/rgb_insertion_annotation/dataset_2022/skeletons/') if 'ZA22' in f]
+
+img_files_pb = []
+for img_file in img_files:
+    plantid, task = img_file.split('.png')[0].split('_')[1:-1]
+    n = '{}_{}'.format(plantid, task)
+    if sum([n in f for f in sk_files]) == 0:
+        img_files_pb.append(img_file)
+
+import shutil
+for img_file in img_files_pb:
+    plantid, task, angle = np.array(img_file.split('.png')[0].split('_')[1:]).astype(int)
+    selec = df_db[(df_db['taskid'] == task) & (df_db['plantid'] == plantid)]
+    selec = selec[selec['imgangle'] == angle]
+    for _, row in selec.iterrows():
+        fd = 'data/rgb_insertion_annotation/dataset_2022/images/'
+        f = fd + 'ZA22_{}_{}_{}.png'.format(str(int(row['plantid'])).zfill(4), row['taskid'], row['imgangle'])
+        if os.path.isfile(f):
+            os.rename(f, fd + '0_ZA22_{}_{}_{}.png'.format(str(int(row['plantid'])).zfill(4), row['taskid'], row['imgangle']))
+        # path1 = 'X:/ARCH2022-01-10/{}/{}.png'.format(task, row['imgguid'])
+        # path2 = 'data/ZA22_visualisation/{}_{}.png'.format(img_file[:-4], row['imgangle'])
+        # shutil.copyfile(path1, path2)
+
+# ===== import some csv extracted with the R code from Llorenc ==================================
+
+df = pd.read_csv('data/copy_from_database/images_ZA22.csv')
+df['daydate'] = [d[5:10] for d in df['acquisitiondate']]
+df = df[(df['viewtypeid'] == 2) & df['imgangle'].isin([k*30 for k in range(12)])]  # 12 side image
+# x = df.groupby(['plantid', 'taskid']).size().reset_index()[0]
+
+df2 = pd.read_csv('data/copy_from_database/plants_ZA22.csv')
+df2 = df2[df2['plantcode'].str.contains('ARCH2022-01-10')]  # remove strange plant names
+df2['genotype'] = [p.split('/')[1] for p in df2['plantcode']]
+df2['scenario'] = [p.split('/')[5] for p in df2['plantcode']]
+
+# ===== deepcollar training dataset ==============================================================
+
+# plant set for deepcollar training
+
+genotypes = np.random.choice(np.unique(df2['genotype']), 250, replace=False)
+
+plants = []
+for g in genotypes:
+    for s in ['WW', 'WD']:
+        s = df2[(df2['genotype'] == g) & (df2['scenario'] == s)]
+        plants.append(np.random.choice(s['plantcode']))
+
+# load the corresponding images
+
+for plant in plants:
+    plantid = int(plant.split('/')[0])
+    if plantid > 425: # re-add old dates
+        s = df[df['plantid'] == plantid]
+        if np.random.random() < 0.1:
+            s = s[s['daydate'] == max(s['daydate'])]
+        else:
+            s = s[s['daydate'] != max(s['daydate'])]
+        row = s.sample().iloc[0]
+        name = 'ZA22_{}_{}_{}.png'.format(str(int(row['plantid'])).zfill(4), row['taskid'], row['imgangle'])
+        path1 = 'Y:/ARCH2022-01-10/{}/{}.png'.format(row['taskid'], row['imgguid'])
+        path2 = 'data/rgb_insertion_annotation/dataset_2022/images/' + name
+        #img = io.imread(path1)[:, :, :3]
+        shutil.copyfile(path1, path2)
+
+# plant set for tracking validation
+
+genotypes2 = np.random.choice([g for g in np.unique(df2['genotype']) if g not in genotypes], 100, replace=False)
+
+plants2 = []
+for g in genotypes2:
+    s = df2[df2['genotype'] == g]
+    plants2.append(np.random.choice(s['plantcode']))
+
+# ===== visu plant for a given plantid ======================================================================
+
+import pandas as pd
+import shutil
+df = pd.read_csv('data/copy_from_database/images_ZA22.csv')
+
+plantid = 1846
+angle = 60
+
+s = df[(df['plantid'] == plantid) & (df['imgangle'] == angle)]
+
+for _, row in s.iterrows():
+    name = 'ZA22_{}_{}_{}.png'.format(str(int(row['plantid'])).zfill(4), row['taskid'], row['imgangle'])
+    path1 = 'X:/ARCH2022-01-10/{}/{}.png'.format(row['taskid'], row['imgguid'])
+    path2 = 'data/ZA22_visualisation/' + name
+    shutil.copyfile(path1, path2)
+
+# ===== other data visu... ==================================================================================
+
+for g in s['genotype'].unique():
+    s2 = s[s['genotype'] == g]
+    for plantid in s2['plantid']:
+        s3 = df[df['plantid'] == plantid]
+
+        # tasks with 12 side views
+        gb = s3.groupby('taskid').agg('count').reset_index()
+        s3 = s3[s3['taskid'].isin(gb[gb['studyid'] == 12]['taskid'])]
+        # angle = 60
+        s3 = s3[s3['imgangle'] == 60]
+
+        print(g, len(s3[s3['imgangle'] == 60]), max(s3['daydate']))
+
+        print(np.array(sorted(s3['daydate'])))
+
+for d in sorted(df['daydate'].unique()):
+    s = df[df['daydate'] == d]
+    print(d, len(s['plantid'].unique()))
+
+
+import os
+p = 'X:/lepseBinaries/ARCH2022-01-10'
+tasks = [t for t in os.listdir(p) if '.' not in t]
+task_to_subdir = {int(t): next(s for s in subdirs if s in os.listdir(p + '/' + t)) for t in tasks}
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/scripts/pheno_multi_exp.py b/scripts/pheno_multi_exp.py
new file mode 100644
index 0000000..2babb91
--- /dev/null
+++ b/scripts/pheno_multi_exp.py
@@ -0,0 +1,64 @@
+import time
+import datetime
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+def date_to_timestamp(date):
+    return time.mktime(datetime.datetime.strptime(date[:10], "%Y-%m-%d").timetuple())
+
+# ===== generate a csv with phenology data from different maize experiments ===========================================
+
+df = []
+
+exp = 'ZB14'
+df_pheno = pd.read_csv('data/pheno_{}.csv'.format(exp))
+df_pheno['Date2'] = ['{}-{}-{}'.format(d.split('/')[2], d.split('/')[1], d.split('/')[0]) for d in df_pheno['Date']]
+df_pheno['timestamp'] = df_pheno['Date2'].map(date_to_timestamp)
+s = df_pheno[~df_pheno['Pheno'].isna()]
+for _, row in s.iterrows():
+    df.append([exp, int(row['ID code'][:4]), row['Date2'], row['timestamp'], 'visi_f', row['Pheno']])
+s = df_pheno[~df_pheno['Ligulee'].isna()]
+for _, row in s.iterrows():
+    df.append([exp, int(row['ID code'][:4]), row['Date2'], row['timestamp'], 'ligul_f', row['Ligulee']])
+
+exp = 'ZA17'
+df_pheno = pd.read_csv('data/pheno_{}.csv'.format(exp))
+df_pheno['timestamp'] = df_pheno['resultdate'].map(date_to_timestamp)
+s = df_pheno[df_pheno['observationcode'].isin(['visi_f', 'ligul_f'])]
+for _, row in s.iterrows():
+    df.append([exp, row['Pot'], row['resultdate'][:10], row['timestamp'], row['observationcode'], row['observation']])
+
+exp = 'ZA22'
+df_pheno = pd.read_csv('data/pheno_{}.csv'.format(exp))
+df_pheno['timestamp'] = df_pheno['resultdate'].map(date_to_timestamp)
+s = df_pheno[df_pheno['observationcode'].isin(['visi_f', 'ligul_f', 'leaf_tot'])]
+s = s[~((s['observationcode'] == 'visi_f') & (s['observation'].isin(['n', 'y', '12..5', '16..5'])))]
+s = s[~((s['observationcode'] == 'ligul_f') & (s['observation'].isin(['pas demar', '14.5', '10.5', '16.8', '610.5'])))]
+s = s[~((s['observationcode'] == 'leaf_tot') & (s['observation'].isin([k for k in s['observation'].unique() if not k.isdigit()])))]
+s['observation'] = pd.to_numeric(s['observation'])
+for _, row in s.iterrows():
+    df.append([exp, row['Pot'], row['resultdate'][:10], row['timestamp'], row['observationcode'], row['observation']])
+
+df = pd.DataFrame(df, columns=['exp', 'plantid', 'daydate', 'timestamp', 'observationcode', 'observation'])
+df.to_csv('data/pheno.csv', index=False)
+
+# ===== visualize data =============================================================================================
+
+df = pd.read_csv('data/pheno.csv')
+
+for code in ['visi_f', 'ligul_f']:
+    plt.figure()
+    for exp in df['exp'].unique():
+        s = df[(df['exp'] == exp) & (df['observationcode'] == code)]
+        plt.plot(s['timestamp'] - np.min(s['timestamp']), s['observation'], 'o', alpha=0.03, label=exp)
+    leg = plt.legend()
+    for lh in leg.legendHandles:
+        lh._legmarker.set_alpha(1)
+    plt.title(code)
+
+
+
+
+
diff --git a/scripts/phenomenal_gif.py b/scripts/phenomenal_gif.py
index b8f7397..439b446 100644
--- a/scripts/phenomenal_gif.py
+++ b/scripts/phenomenal_gif.py
@@ -12,8 +12,8 @@
 
 #local copy of phenoarch_cache/ZA17
 cache_dir = 'local_cache/cache_ZA17/'
-#folder =  cache_dir + 'segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
-folder = cache_dir + 'segmentation_voxel4_tol1_notop_vis4_minpix100_no_stem_smooth_no_tracking'
+folder =  cache_dir + 'segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking'
+#folder = cache_dir + 'segmentation_voxel4_tol1_notop_vis4_minpix100_no_stem_smooth_no_tracking'
 #folder = cache_dir + '1429_nosmooth'
 
 plantid = 1429
@@ -131,7 +131,7 @@
 #imgs_gif = [img[200:-200, 250:-250, :] for img in imgs_gif]
 imgs_gif = [Image.fromarray(np.uint8(img)) for img in imgs_gif]
 fps = 1
-imgs_gif[0].save('gif/animation_id{}_{}fps.gif'.format(plantid, fps),
+imgs_gif[0].save('gif/animation_id{}_{}fps_2022.gif'.format(plantid, fps),
               save_all=True,
               append_images=imgs_gif[1:],
               optimize=True,
diff --git a/scripts/stem_tracking_gif.py b/scripts/stem_tracking_gif.py
new file mode 100644
index 0000000..9ae5595
--- /dev/null
+++ b/scripts/stem_tracking_gif.py
@@ -0,0 +1,66 @@
+from openalea.maizetrack.local_cache import get_metainfos_ZA17
+from openalea.maizetrack.utils import get_rgb
+
+from openalea.phenomenal import object as phm_obj
+from openalea.maizetrack.utils import phm3d_to_px2d
+
+import os
+
+from PIL import Image
+
+
+vmsi_folder = 'local_cache/cache_ZA17/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking/'
+all_vmsis = [vmsi_folder + d + '/' + f for d in os.listdir(vmsi_folder) for f in os.listdir(vmsi_folder + d)]
+vmsi_files = [v for v in all_vmsis if int(v.split('/')[-1][:4]) == 1429]
+
+vmsis = {}
+for f in vmsi_files:
+    task = int(f.split('__')[-1].split('.')[0])
+    print(task)
+    vmsis[task] = phm_obj.VoxelSegmentation.read_from_json_gz(f)
+
+plantid = 1429
+metainfos = get_metainfos_ZA17(plantid)
+
+images = {}
+for task, sk in vmsis.items():
+    print(task)
+    m = next(m for m in metainfos if m.task == task)
+    img, _ = get_rgb(metainfo=m, angle=60, save=False, plant_folder=False)
+    images[task] = img
+
+#all_collars = pd.read_csv('data/visu_3D/stem_tracking/1429.csv')
+
+imgs = []
+T, H = [], []
+for task in sorted(vmsis.keys()):
+    print(task)
+    vmsi, img = vmsis[task], images[task]
+    m = next(m for m in metainfos if m.task == task)
+    #collars = all_collars[all_collars['t'] == m.timestamp]
+    xyz = vmsi.get_stem().info['pm_position_tip']
+    x, y = phm3d_to_px2d(xyz, sf=m.shooting_frame)[0]
+    T.append(m.timestamp)
+    H.append(y)
+    img = cv2.rectangle(np.float32(img), (int(x - 25), int(y - 25)), (int(x + 25), int(y + 25)), (255, 0, 0), 2)
+    img = cv2.circle(np.float32(img), (int(x), int(y)), 7, (255, 0, 0), -1)
+    #io.imsave('data/visu_3D/{}.png'.format(task), np.uint8(img))
+    imgs.append(img)
+
+imgs_gif = imgs[:-7].copy()
+#imgs_gif = [img[200:-200, 250:-250, :] for img in imgs_gif]
+imgs_gif = [Image.fromarray(np.uint8(img)) for img in imgs_gif]
+fps = 4
+imgs_gif[0].save('gif/stem_tracking_1429_{}fps.gif'.format(fps),
+              save_all=True,
+              append_images=imgs_gif[1:],
+              optimize=True,
+              duration=1000/fps,
+              loop=0)
+
+
+
+
+
+
+
diff --git a/scripts/tracking_develop.py b/scripts/tracking_develop.py
new file mode 100644
index 0000000..7dfe1b8
--- /dev/null
+++ b/scripts/tracking_develop.py
@@ -0,0 +1,50 @@
+import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from openalea.phenomenal import object as phm_obj
+from openalea.maizetrack.utils import missing_data
+
+from maizetrack.scripts.visu_3d import save_image
+from openalea.maizetrack.phenomenal_display import plot_vmsi
+
+PATH = 'data/copy_from_modulor/'
+
+exp = 'ZB14'
+
+seg_path = PATH + 'cache_{}/segmentation_voxel4_tol1_notop_vis4_minpix100_stem_smooth_tracking/'.format(exp)
+
+all_files = [seg_path + d + '/' + f for d in os.listdir(seg_path) for f in os.listdir(seg_path + d)]
+plantids = np.unique([int(f.split('/')[-1][:4]) for f in all_files])
+
+# contains daydate, task, timestamp, shooting_frame
+index = pd.read_csv(PATH + 'cache_{0}/snapshot_index_{0}.csv'.format(exp))
+
+for plantid in plantids:
+
+    print('=====================')
+
+    files = [f for f in all_files if int(f.split('/')[-1][:4]) == plantid]
+
+    for file in files:
+    #for file in np.random.choice(files, 5, replace=False):
+
+        plant = next(p for p in index['plant'] if int(p[:4]) == plantid)
+        task = int(file.split('/')[-1].split('.json')[0].split('_')[-1])
+        index_row = index[(index['plant'] == plant) & (index['task'] == task)].iloc[0]
+
+        print(plantid, index_row.daydate)
+        vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(file)
+
+        scene = plot_vmsi([vmsi], )
+        _ = save_image(scene, image_name='data/visu_3D/images_ZB14/{}_{}.png'.format(index_row['daydate'], plantid))
+
+
+
+
+
+
+
+
+
+
diff --git a/scripts/visu_3d.py b/scripts/visu_3d.py
new file mode 100644
index 0000000..391c467
--- /dev/null
+++ b/scripts/visu_3d.py
@@ -0,0 +1,144 @@
+import os
+import numpy as np
+import pandas as pd
+from skimage import io
+import matplotlib.pyplot as plt
+
+from openalea.phenomenal import object as phm_obj
+from openalea.maizetrack.phenomenal_display import plot_vmsi, plot_snapshot
+from openalea.plantgl.all import Viewer
+from openalea.maizetrack.trackedPlant import TrackedSnapshot
+from openalea.maizetrack.utils import missing_data, phm3d_to_px2d
+
+from openalea.phenomenal.calibration import Calibration
+
+
+# function from openalea-incubator/adel
+def save_image(scene, image_name='%s/img%04d.%s', directory='.', index=0, ext='png'):
+    '''
+    Save an image of a scene in a specific directory
+    Parameters
+    ----------
+        - scene: a PlantGL scene
+        - image_name: a string template
+            The format of the string is dir/img5.png
+        - directory (optional: ".") the directory where the images are written
+        - index: the index of the image
+        - ext : the image format
+    Example
+    -------
+        - Movie:
+            convert *.png movie.mpeg
+            convert *.png movie.gif
+            mencoder "mf://*.png" -mf type=png:fps=25 -ovc lavc -o output.avi
+            mencoder -mc 0 -noskip -skiplimit 0 -ovc lavc -lavcopts vcodec=msmpeg4v2:vhq "mf://*.png" -mf type=png:fps=18 -of avi  -o output.avi
+
+    '''
+
+    if not image_name:
+        image_name = '{directory}/img{index:0>4d}.{ext}'
+    filename = image_name.format(directory=directory, index=index, ext=ext)
+    Viewer.frameGL.saveImage(filename)
+    return scene,
+
+
+if __name__ == '__main__':
+
+
+    # ===== skeleton reprojected on rgb =================================================================================
+
+    df_db = pd.read_csv('data/copy_from_database/images_ZA22.csv')
+    df_sf = pd.read_csv('data/copy_from_modulor/snapshot_index_ZA22.csv')
+
+    rgb_path = 'X:/ARCH2022-01-10/'
+
+    sk_path = 'data/copy_from_modulor/skeleton_voxel4_tol1_notop_vis4_minpix100/'
+    sk_files = [sk_path + d + '/' + f for d in os.listdir(sk_path) for f in os.listdir(sk_path + d)]
+
+    sk_file = np.random.choice(sk_files)
+    sk = phm_obj.VoxelSkeleton.read_from_json_gz(sk_file)
+
+    plantid, task = np.array(sk_file.split('/')[-1].split('.json.gz')[0].split('_'))[[0, -1]].astype(int)
+
+    angle = 60
+    selec_db = df_db[(df_db['taskid'] == task) & (df_db['plantid'] == plantid) & (df_db['imgangle'] == angle)]
+    img = io.imread(rgb_path + str(task) + '/' + selec_db.iloc[0]['imgguid'] + '.png')[:, :, :3]
+
+    plant = next(p for p in df_sf['plant'].unique() if int(p.split('/')[0]) == plantid)
+    selec_sf = df_sf[(df_sf['task'] == task) & (df_sf['plant'] == plant)]
+    sf = selec_sf.iloc[0]['shooting_frame']
+
+    calib = Calibration.load('V:/lepseBinaries/Calibration/' + sf + '_calibration.json')
+    f_calib = calib.get_projection(id_camera='side', rotation=angle)
+
+    px_polylines = [f_calib(seg.polyline) for seg in sk.segments]
+
+    print(sf)
+    plt.imshow(img)
+    vxs = f_calib(sk.voxels_position())
+    plt.plot(vxs[:, 0], vxs[:, 1], 'w.', markersize=0.5)
+    for pl in px_polylines:
+        plt.plot(pl[:, 0], pl[:, 1], 'r-')
+
+    # ==========================================================================================
+
+    exp = 'ZB14'
+
+    PATH = 'data/visu_3D/'
+    vmsi_folder = 'set_{}_vmsi/'.format(exp)
+    image_folder = 'images_{}/'.format(exp)
+
+    if not os.path.isdir(PATH + image_folder):
+        os.mkdir(PATH + image_folder)
+
+    files = os.listdir(PATH + vmsi_folder)
+
+    for k, file in enumerate(files):
+        print(k)
+
+        vmsi = phm_obj.VoxelSegmentation.read_from_json_gz(PATH + vmsi_folder + file)
+
+        scene = plot_vmsi([vmsi], )
+
+        if not missing_data(vmsi):
+            rank_to_index = {l.info['pm_leaf_number_tracking']: k for k, l in enumerate(vmsi.get_leafs())}
+            mature_ranks = [-1 + vmsi.get_leaf_order(k).info['pm_leaf_number_tracking']
+                            if vmsi.get_leaf_order(k).info['pm_label'] == 'mature_leaf' else -1
+                            for k in range(1, 1 + vmsi.get_number_of_leaf())]
+            snapshot = TrackedSnapshot(vmsi, metainfo=None, order=None)
+            scene = plot_snapshot(snapshot, ranks=mature_ranks)
+
+            _ = save_image(scene, image_name=PATH + image_folder + '{}.png'.format(file))
+
+        else:
+            print('missing data', file)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/setup.py b/setup.py
index 4293c12..caea2e7 100644
--- a/setup.py
+++ b/setup.py
@@ -29,7 +29,7 @@
 
 setup(
     name="openalea.maizetrack",
-    version="0.0.1",
+    version="0.1.3",
     description="",
     long_description="",
 
diff --git a/src/openalea/maizetrack/alignment.py b/src/openalea/maizetrack/alignment.py
index 34333c1..c4b08af 100644
--- a/src/openalea/maizetrack/alignment.py
+++ b/src/openalea/maizetrack/alignment.py
@@ -297,7 +297,7 @@ def polylines_distance(pl1, pl2, n=20):
     return dist
 
 
-def phm_leaves_distance(leaf_ref, leaf_candidate, method):
+def phm_leaves_distance(leaf_ref, leaf_candidate):
     """ takes two leaf objects, deduct two polylines which start from a same point, compute the distance between
      these two polylines """
 
@@ -311,12 +311,12 @@ def phm_leaves_distance(leaf_ref, leaf_candidate, method):
     else:
         pl1 = polyline_until_z(leaf1.highest_pl, zbase2)
 
-    #normalize
+    # normalize
     len1 = np.sum([np.linalg.norm(np.array(pl1[k]) - np.array(pl1[k + 1])) for k in range(len(pl1) - 1)])
     pl1 = pl1 / np.max((len1, 0.0001))
     pl2 = pl2 / np.max((len1, 0.0001))
 
     # computing distance
-    d = polylines_distance(pl1, pl2, method)
+    d = polylines_distance(pl1, pl2)
 
     return d
\ No newline at end of file
diff --git a/src/openalea/maizetrack/phenomenal_display.py b/src/openalea/maizetrack/phenomenal_display.py
index 9878851..de20d54 100644
--- a/src/openalea/maizetrack/phenomenal_display.py
+++ b/src/openalea/maizetrack/phenomenal_display.py
@@ -160,9 +160,9 @@ def plot_vmsi(vmsi_list, only_mature=False):
         li = 'all'
         coli = 'same'
         shapes += vmsi_polylines_to_pgl(vmsi, li, coli, only_mature)
-    print(len(shapes), ' shapes to plot')
     scene = pgl.Scene(shapes)
     pgl.Viewer.display(scene)
+    return scene
 
 
 def plot_vmsi_voxel(vmsi, ranks=None):
@@ -345,4 +345,6 @@ def plot_snapshot(snapshot, colored=True, ranks=None, stem=True):
             shapes.append(cyl)
 
     scene = pgl.Scene(shapes)
-    pgl.Viewer.display(scene)
\ No newline at end of file
+    pgl.Viewer.display(scene)
+
+    return scene
\ No newline at end of file
diff --git a/src/openalea/maizetrack/stem_correction.py b/src/openalea/maizetrack/stem_correction.py
index 0159037..8d3fba2 100644
--- a/src/openalea/maizetrack/stem_correction.py
+++ b/src/openalea/maizetrack/stem_correction.py
@@ -94,6 +94,7 @@ def abnormal_stem(vmsi_list, dist_threshold=100):
 
 # ============================================================================================================
 
+
 def xyz_last_mature(vmsi):
 
     # xyz position of the last mature leaf insertion
@@ -117,7 +118,7 @@ def savgol_smoothing_function(x, y, dw, polyorder, repet):
     # x2, y2 = savgol_filter((x2, y2), window_length=w, polyorder=polyorder)
     x2, y2 = x, y
     for k in range(repet):
-        x2, y2 = savgol_filter((x2, y2), window_length=w, polyorder=polyorder)
+        x2, y2 = savgol_filter((x2, y2), window_length=max((w, polyorder + 1)), polyorder=polyorder)
 
     # TODO : bricolage ! trouver une fonction de lissage monotone
     # monotony
diff --git a/src/openalea/maizetrack/trackedPlant.py b/src/openalea/maizetrack/trackedPlant.py
index 6b1f123..3240cf7 100644
--- a/src/openalea/maizetrack/trackedPlant.py
+++ b/src/openalea/maizetrack/trackedPlant.py
@@ -192,7 +192,7 @@ def simplify_polylines(self, new_length=50):
          computation of growing leaves tracking """
         for snapshot in self.snapshots:
             for leaf in snapshot.leaves:
-                # polyline starting from insertion point
+                # polyline starting from insertion point # TODO : or the contrary ?
                 leaf.highest_pl = simplify(leaf.get_highest_polyline().polyline, new_length)
                 # polyline starting from stem base
                 leaf.real_pl = simplify(leaf.real_longest_polyline(), new_length)
@@ -225,14 +225,14 @@ def get_ref_skeleton(self, display=False, nmax=15):
             # remove old leaves (that could have a different shape)
             leaves = leaves[:nmax]
 
-            vec_mean = np.mean(np.array([l.vec for l in leaves]), axis=0)
-            dists = [np.sum(abs(l.vec - vec_mean)) for l in leaves]
-
-            ref_skeleton[rank] = leaves[np.argmin(dists)]
+            if leaves:
+                vec_mean = np.mean(np.array([l.vec for l in leaves]), axis=0)
+                dists = [np.sum(abs(l.vec - vec_mean)) for l in leaves]
+                ref_skeleton[rank] = leaves[np.argmin(dists)]
 
         if display:
             if has_pgl_display:
-                plot_leaves([ref_skeleton[r] for r in ranks], ranks)
+                plot_leaves(list(ref_skeleton.values()), list(ref_skeleton.keys()))
             else:
                 warnings.warn('PlantGL not available')
 
@@ -251,7 +251,7 @@ def tracking_info_update(self):
                 leaf.info['pm_leaf_number_tracking'] = rank + 1
 
     def align_mature(self, gap=12.35, gap_extremity_factor=0.2, direction=1, n_previous=5000, w_h=0.03, w_l=0.004,
-                     rank_attribution=True):
+                     rank_attribution=True, remove_senescence=False):
         """
         alignment and rank attributions in a time-series of sequences of leaves.
         Step 1 : use a multiple sequence alignment algorithm to align the sequences.
@@ -289,7 +289,13 @@ def align_mature(self, gap=12.35, gap_extremity_factor=0.2, direction=1, n_previ
 
         # init order attribute of each snapshot, with only mature leaves:
         for snapshot in self.snapshots:
-            snapshot.order = [i for i, l in enumerate(snapshot.leaves) if l.info['pm_label'] == 'mature_leaf']
+            # TODO maj 15/06/2022, only for ZA17
+            if remove_senescence:
+                zbase = {'elcom_2_c1_wide': -730, 'elcom_2_c2_wide': -690}
+                snapshot.order = [i for i, l in enumerate(snapshot.leaves) if l.info['pm_label'] == 'mature_leaf' and
+                                  l.real_longest_polyline()[-1][2] - zbase[snapshot.metainfo.shooting_frame] > 10]
+            else:
+                snapshot.order = [i for i, l in enumerate(snapshot.leaves) if l.info['pm_label'] == 'mature_leaf']
 
         # time-series of sequences of vectors (sequences have different lengths, all vectors have the same size)
         sequences = [np.array([snapshot.leaves[i].vec for i in snapshot.order]) for snapshot in self.snapshots]
@@ -299,7 +305,9 @@ def align_mature(self, gap=12.35, gap_extremity_factor=0.2, direction=1, n_previ
 
         # update order attributes
         for i, order in enumerate(alignment_matrix):
-            self.snapshots[i].order = list(order)
+            #self.snapshots[i].order = list(order)
+            # TODO maj 15/06/2022 /!\
+            self.snapshots[i].order = [o if o == -1 else self.snapshots[i].order[o] for o in order]
 
         # Step 2 - abnormal ranks removing
         # ==============================================
diff --git a/src/openalea/maizetrack/utils.py b/src/openalea/maizetrack/utils.py
index e83ef9b..c50bf76 100644
--- a/src/openalea/maizetrack/utils.py
+++ b/src/openalea/maizetrack/utils.py
@@ -8,8 +8,10 @@
 from openalea.maizetrack.phenomenal_display import PALETTE
 from skimage import io
 
+from openalea.phenomenal.calibration import Calibration
 
-##### based on phenomenal #####
+
+# ===== based on Phenomenal =========================================================================================
 
 def shooting_frame_conversion(s_frame_name):
     # returns 2 parameters needed for this conversion :
@@ -35,7 +37,15 @@ def phm3d_to_px2d(xyz, sf, angle=60):
     if xyz.ndim == 1:
         xyz = np.array([xyz])
 
-    f = get_shooting_frame(sf).get_calibration('side').get_projection(angle)
+    if sf.startswith('ARCH'):
+        # TODO : temporary
+        # phenomenal function
+        calib = Calibration.load('V:/lepseBinaries/Calibration/' + sf + '_calibration.json')
+        f = calib.get_projection(id_camera='side', rotation=angle)
+    else:
+        # phenoarch function
+        f = get_shooting_frame(sf).get_calibration('side').get_projection(angle)
+
     return f(xyz)
 
 
@@ -104,7 +114,7 @@ def rgb_and_polylines(snapshot, angle, selected=None, ranks=None):
     return img, polylines_px
 
 
-##### polyline analysis #####
+# ===== polyline analysis =====================================================================================
 
 def quantile_point(pl, q):
     pl = np.array(pl)
@@ -218,7 +228,8 @@ def missing_data(vmsi):
 
 def dataset_mean_distance(w_h=0.03, w_l=0.004, step=1):
     """
-    mean distance between consecutive leaves (spatially) in a small dataset
+    mean distance between consecutive leaves (spatially) in a small dataset.
+    file leaf_vectors.npy generated using 30 random plants (available in modulor local_benoit)
     w_h=0.03, w_l=0.004 => d = 4.23
     """
     v = np.load('leaf_vectors.npy', allow_pickle=True)