{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthetic Validation\n",
    "\n",
    "## Import Statement :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cv2\n",
    "import collections\n",
    "import numpy\n",
    "import math\n",
    "import pandas\n",
    "import scipy.integrate\n",
    "\n",
    "import openalea.phenomenal.data as phm_data\n",
    "import openalea.phenomenal.object as phm_obj\n",
    "import openalea.phenomenal.multi_view_reconstruction as phm_mvr\n",
    "import openalea.phenomenal.mesh as phm_mesh\n",
    "import openalea.phenomenal.display as phm_display\n",
    "import openalea.phenomenal.display.notebook as phm_display_notebook\n",
    "import openalea.phenomenal.segmentation as phm_seg\n",
    "from openalea.phenotyping_data.fetch import fetch_all_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Import synthetic data\n",
    "\n",
    "### 1.1 Select plant_number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibration_dir = fetch_all_data('plant_1/calibration').parent\n",
    "data_dir = fetch_all_data('synthetic_plant_1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vertices, faces, meta_data = phm_data.synthetic_plant(data_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Viewing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phm_display_notebook.show_synthetic_plant(vertices, faces, meta_data=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Phenomenal Measurements\n",
    "\n",
    "### 2.1 Projection & binairies images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_shape = (2448, 2048)\n",
    "angles = range(0, 360, 30)\n",
    "calibrations = phm_data.calibrations(calibration_dir)\n",
    "\n",
    "selected = [(\"side\", angle, (2448, 2048)) for angle in range(0, 360, 30)]\n",
    "selected.append((\"top\", 0, (2048, 2448)))\n",
    "\n",
    "\n",
    "bin_images = collections.defaultdict(dict)\n",
    "for id_camera, angle, image_shape in selected:\n",
    "    projection = calibrations[id_camera].get_projection(angle)\n",
    "    image = phm_mesh.project_mesh_on_image(vertices, faces, image_shape, projection)\n",
    "    bin_images[id_camera][angle] = image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "phm_display.show_images(\n",
    "    list(bin_images[\"side\"].values()) + list(bin_images[\"top\"].values())\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Multi-view reconstruction\n",
    "\n",
    "#### 2.2.1 Associate images and projection function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def routine_select_ref_angle(bin_side_images):\n",
    "    max_len = 0\n",
    "    max_angle = None\n",
    "\n",
    "    for angle in bin_side_images:\n",
    "        x_pos, y_pos, x_len, y_len = cv2.boundingRect(\n",
    "            cv2.findNonZero(bin_side_images[angle])\n",
    "        )\n",
    "\n",
    "        if x_len > max_len:\n",
    "            max_len = x_len\n",
    "            max_angle = angle\n",
    "\n",
    "    return max_angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "refs_angle_list = [routine_select_ref_angle(bin_images[\"side\"])]\n",
    "\n",
    "image_views = list()\n",
    "for id_camera in bin_images:\n",
    "    for angle in bin_images[id_camera]:\n",
    "        projection = calibrations[id_camera].get_projection(angle)\n",
    "\n",
    "        image_ref = None\n",
    "        if id_camera == \"side\" and angle in refs_angle_list:\n",
    "            image_ref = bin_images[id_camera][angle]\n",
    "\n",
    "        inclusive = False\n",
    "        if id_camera == \"top\":\n",
    "            inclusive = True\n",
    "\n",
    "        image_views.append(\n",
    "            phm_obj.ImageView(\n",
    "                bin_images[id_camera][angle],\n",
    "                projection,\n",
    "                inclusive=inclusive,\n",
    "                image_ref=image_ref,\n",
    "            )\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.2 Do the mul-view reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_tolerance = 1\n",
    "voxels_size = 4\n",
    "\n",
    "voxel_grid = phm_mvr.reconstruction_3d(\n",
    "    image_views, voxels_size=voxels_size, error_tolerance=error_tolerance\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "phm_display_notebook.show_voxel_grid(voxel_grid, size=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Skeletonization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph = phm_seg.graph_from_voxel_grid(voxel_grid, connect_all_point=True)\n",
    "voxel_skeleton = phm_seg.skeletonize(voxel_grid, graph)\n",
    "\n",
    "# Select images\n",
    "image_projection = list()\n",
    "for angle in bin_images[\"side\"]:\n",
    "    projection = calibrations[\"side\"].get_projection(angle)\n",
    "    image_projection.append((bin_images[\"side\"][angle], projection))\n",
    "\n",
    "voxel_skeleton_denoised = phm_seg.segment_reduction(\n",
    "    voxel_skeleton, image_projection, required_visible=4, nb_min_pixel=100\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "phm_display_notebook.show_skeleton(voxel_skeleton_denoised, with_voxel=True, size=1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Cereals Segmentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vms = phm_seg.maize_segmentation(voxel_skeleton_denoised, graph)\n",
    "vmsi = phm_seg.maize_analysis(vms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Phenomenal measurements of each organs\n",
    "pm_rows = [vo.info for vo in vmsi.voxel_organs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "phm_display_notebook.show_segmentation(vmsi, size=1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Synthetic measurement\n",
    "\n",
    "### 3.1 Diplay leaf meta-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phm_display_notebook.show_synthetic_plant(vertices, faces, meta_data=meta_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Extract measurements from meta-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_duplicate(poly, width):\n",
    "    index, d = list(), dict()\n",
    "    for i, item in enumerate(map(tuple, list(poly))):\n",
    "        if item not in d:\n",
    "            index.append(i)\n",
    "            d[item] = True\n",
    "    return poly[index], width[index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm_rows = list()\n",
    "x, y, z, r = zip(*meta_data[\"stem\"])\n",
    "polyline = numpy.array(list(zip(x, y, z))) * 10 - numpy.array([0, 0, 750])\n",
    "r = numpy.array(r) * 10\n",
    "polyline, r = remove_duplicate(polyline, r)\n",
    "\n",
    "row = dict()\n",
    "row[\"sm_label\"] = \"stem\"\n",
    "row[\"sm_max_width\"] = max(r)\n",
    "row[\"sm_average_width\"] = numpy.mean(r)\n",
    "\n",
    "row[\"sm_surface\"] = scipy.integrate.simpson(\n",
    "    r, x=phm_seg.compute_curvilinear_abscissa(polyline, 1)\n",
    ")\n",
    "row[\"sm_length\"] = phm_seg.compute_length_organ(polyline)\n",
    "angle, stem_vector_mean = phm_seg.compute_azimuth_angle(polyline)\n",
    "sm_rows.append(row)\n",
    "\n",
    "ranks = meta_data[\"leaf_order\"]\n",
    "polylines = {\n",
    "    n: map(numpy.array, list(zip(*meta_data[\"leaf_polylines\"][i])))\n",
    "    for i, n in enumerate(ranks)\n",
    "}\n",
    "\n",
    "for leaf_order in polylines:\n",
    "    x, y, z, r = polylines[leaf_order]\n",
    "    polyline = numpy.array(list(zip(x, y, z))) * 10 - numpy.array([0, 0, 750])\n",
    "    r = numpy.array(r) * 10\n",
    "    polyline, r = remove_duplicate(polyline, r)\n",
    "\n",
    "    row = dict()\n",
    "    row[\"sm_leaf_number\"] = leaf_order\n",
    "    row[\"sm_label\"] = \"leaf\"\n",
    "    row[\"sm_position_tip\"] = tuple(polyline[-1])\n",
    "    row[\"sm_position_base\"] = tuple(polyline[0])\n",
    "    row[\"sm_max_width\"] = max(r)\n",
    "    row[\"sm_average_width\"] = numpy.mean(r)\n",
    "\n",
    "    row[\"sm_surface\"] = scipy.integrate.simpson(\n",
    "        r, x=phm_seg.compute_curvilinear_abscissa(polyline, 1)\n",
    "    )\n",
    "\n",
    "    row[\"sm_length\"] = phm_seg.compute_length_organ(polyline)\n",
    "    angle, vector_mean = phm_seg.compute_azimuth_angle(polyline)\n",
    "    row[\"sm_azimuth_angle\"] = angle\n",
    "    insertion_angle, vector = phm_seg.compute_insertion_angle(\n",
    "        polyline, stem_vector_mean\n",
    "    )\n",
    "    row[\"sm_insertion_angle_vector\"] = vector\n",
    "    row[\"sm_insertion_angle\"] = insertion_angle\n",
    "    row[\"sm_inclination_angle\"] = phm_seg.compute_inclination_angle(polyline)\n",
    "    row[\"sm_full_length\"] = row[\"sm_length\"] + polyline[0][2] + 750\n",
    "\n",
    "    sm_rows.append(row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Measures registration\n",
    "\n",
    "Merge synthetic and phenomenal row measure together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def registration_leafs(pm_rows, sm_rows):\n",
    "    def cannot_registred(pm_row, sm_row):\n",
    "        for label in [\"stem\", \"plant\", \"unknown\"]:\n",
    "            if pm_row[\"pm_label\"] == label or sm_row[\"sm_label\"] == label:\n",
    "                return True\n",
    "        if math.isnan(pm_row[\"pm_length\"]):\n",
    "            return True\n",
    "        return False\n",
    "\n",
    "    def same_label(pm_row, sm_row):\n",
    "        for label in [\"stem\", \"plant\"]:\n",
    "            if pm_row[\"pm_label\"] == label and sm_row[\"sm_label\"] == label:\n",
    "                return True\n",
    "        return False\n",
    "\n",
    "    rows_distance = []\n",
    "    registered_row = []\n",
    "    for pm_row in pm_rows:\n",
    "        for sm_row in sm_rows:\n",
    "            if same_label(pm_row, sm_row):\n",
    "                registered_row.append((pm_row, sm_row, 0))\n",
    "                continue\n",
    "            if cannot_registred(pm_row, sm_row):\n",
    "                continue\n",
    "\n",
    "            pos1 = numpy.array(pm_row[\"pm_position_tip\"])\n",
    "            pos2 = numpy.array(sm_row[\"sm_position_tip\"])\n",
    "            distance = numpy.linalg.norm(pos2 - pos1)\n",
    "\n",
    "            rows_distance.append((pm_row, sm_row, distance))\n",
    "\n",
    "    while rows_distance:\n",
    "        (pm_row, sm_row, d) = min(rows_distance, key=lambda t: t[2])\n",
    "        registered_row.append((pm_row, sm_row, d))\n",
    "        rows_distance = [\n",
    "            v for v in rows_distance if pm_row not in v and sm_row not in v\n",
    "        ]\n",
    "    return registered_row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "registred_rows = list()\n",
    "for pm_row, sm_row, d in registration_leafs(pm_rows, sm_rows):\n",
    "    sm_row.update(pm_row)\n",
    "    sm_row[\"distance_registration\"] = d\n",
    "    registred_rows.append(sm_row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Display result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pandas.DataFrame(registred_rows)\n",
    "# df = df[df[\"distance_registration\"] < 100]\n",
    "# df_filtred_label = df[df[\"pm_label\"] == \"mature_leaf\"]\n",
    "\n",
    "# print df[\"sm_inclination_angle\"], df[\"pm_inclination_angle\"],\n",
    "\n",
    "plt.plot(df[\"sm_full_length\"], df[\"pm_full_length\"], \"k+\")\n",
    "plt.plot(df[\"sm_surface\"] / 100.0, df[\"pm_surface\"] / 100.0, \"c+\")\n",
    "plt.plot(df[\"sm_inclination_angle\"] * 100.0, df[\"pm_inclination_angle\"] * 100.0, \"co\")\n",
    "for label, color in [(\"mature_leaf\", \"b+\"), (\"growing_leaf\", \"r+\"), (\"stem\", \"g+\")]:\n",
    "    df_filtred_label = df[df[\"pm_label\"] == label]\n",
    "    plt.plot(df_filtred_label[\"sm_length\"], df_filtred_label[\"pm_length\"], color)\n",
    "\n",
    "plt.plot(df_filtred_label[\"sm_length\"], df_filtred_label[\"pm_length\"], color)\n",
    "plt.xlabel(\"Phenomenal measurements\")\n",
    "plt.ylabel(\"Synthetic measurements\")\n",
    "\n",
    "m = int(\n",
    "    max(\n",
    "        [\n",
    "            df[k].max()\n",
    "            for k in [\"pm_length\", \"sm_length\", \"pm_full_length\", \"sm_full_length\"]\n",
    "        ]\n",
    "    )\n",
    ")\n",
    "plt.plot(range(m), range(m), \"k-\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.17"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}