Merge pull request #2 from openalea/v0_1_3

V0 1 3

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:

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')
+

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)

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))
+
+
+
+
+

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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+