ray_tracing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Le présent travail a été réalisé par une étudiante dans le cadre du cours \"Agriculture en environnement controlé\" donné en Master 2 à Gembloux Agro-bio Tech (Belgique) par le professeur Fédéric Lebeau. Ce travail a été également encadré par Arnaud Bouvry, assistant-doctorant.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#                                     Ray-tracing dans une chambre de culture avec le logiciel VTK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###                                                               Marie Picron   \n",
    "#### Janvier 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Introduction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La croissance démographique entrainant une pression sur la suffisance alimentaire et l'agriculture menacée par les conditions métérologiques et par la disponibilité des ressources (Tian, 2016), il est nécessaire aujourd'hui de repenser à une nouvelle forme d'agriculture. \n",
    "\n",
    "Les *Plant Factories with Artificial Lamps* (PFALs) connaissent un essor depuis quelques décennies. Ces infrastructures constituent à verticaliser les cultures sous forme de chambres et à suivre et contrôler l'environnement et ses divers paramètres climatiques, les entrées et les sorties.\n",
    "\n",
    "Les éclairages artificiels sont des composants des chambres PFALs. Ils apportent une source lumineuse aux cultures. Ils sont primordiales pour la photosynthèse des plantes. Une bonne configuration et une gestion correcte des paramètres de ces éclairages sont nécessaires pour garantir une croissance optimale des plantes.\n",
    "\n",
    "En plus d'approvisonner en alimentation, les PFALs permettent aussi de faciliter la recherche et l'éducation. En effet, des projets émergent grâce à la capacité d'intéragir et de manipuler les différents paramètres sans dépendre d'autres variables. Ainsi, l'étude et la compréhension des différents phénomènes sont de plus en plus précises et complètes. Parmi les études, l'influence des propriétés lumineuses des lampes sur l'éclairement énergétique et sur la croissance optimale des plantes dans la chambre de culture (Fiorillo et al., 2010; Tian, 2016) sont des exemples.\n",
    " \n",
    "Les différentes expérimentations et manipulations dans les chambres peuvent être chronophages et épuisent des ressources. La tendance de l'industrie 4.0 où on tend vers un monde digitalisé, permet d'apporter des solutions. En effet, les PFALs peuvent être virtualisés pour modèliser l'environnement réelle et les phénomènes accompagnant. Ainsi, la recherche peut être rapide et précise grâce à la virtualisation.\n",
    "\n",
    "Parmi ces phénomènes à modéliser, se trouve le ray-tracing des lampes de la chambre. Le ray-tracing consiste à simuler le parcours de la lumière. Il permet de déterminer l'éclairement lumineux des objets en utilisant les lois optiques. Comme précédemment mentionné, les lampes sont indispensables dans une chambre de culture, la modèlisation du ray-tracing est cruciale pour évaluer la quantité de lumière reçue par les plantes.\n",
    "\n",
    "L'objectif du travail présent est de trouver un algorithme permettant de virtualiser le ray-tracing des lampes horticoles. Cet algorithme permettra d'étudier, à long terme, la composition du spectre lumineux reçu par une surface photosynthètiquement active en faisant intervenir des paramètres variables des lampes. \n",
    "\n",
    "Divers logiciels existent pour permettre le ray-tracing comme le *GPUFlux model in GroIMP* (Henke & Buck-Sorlin, 2017) ou *PHOTONSIM* (Kostro et al., 2007). Le travail présent utilise le logiciel _Visualization ToolKit_ (VTK). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ***Le logiciel VTK***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce logiciel open source et puissant, permet de réaliser différentes opérations dans un espace en 3 dimensions comme la visualisation de données volumiques ou le traitement d'image (VTK,2020). \n",
    "\n",
    "VTK permet de manipuler et d'intéragir facilement avec les objets sans se soucier des coordonnées spatiales et des différentes propriétés géométriques. En effet, cette librairie VTK possède un nombre important de packages et d'outils permettant de résoudre différentes opérations. \n",
    "\n",
    "Bien que ce logiciel présente de nombreux avantages, il contient, cependant, aussi de nombreux inconvénients. En effet, pour réaliser une simple opération, il est parfois nécessaire de coder plusieurs lignes et d'utiliser plusieurs fonctions afin d'obtenir l'objectif. Il existe des architectures de hiérarchie de classes à procéder et à respecter avec éventuellement des options (sous-classes) à ajouter. Il faut les associer entre elles ces différentes classes pour obtenir le résultat souhaité. Par exemple, pour visualiser une forme géométrique, la procédure à suivre est la suivante (Schroede et al., 2000) : \n",
    "\n",
    "$$\n",
    "Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \\xrightarrow{\\text{}} Renderer\n",
    "$$\n",
    "\n",
    "Ces différentes classes et sous-classes seront expliquées ultérieurement.\n",
    "                                                  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Rappels Théoriques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est important de mentionner quelques définitions utilisées dans le présent cas.\n",
    "\n",
    "- Un photon est le quantum d'énergie associé aux ondes lumineuses.\n",
    "- Un rayon est défini comme la trajectoire droite d'un paquet contenant une infinité de photons.\n",
    "- La luminance est définie comme le flux lumineux provenant d'une surface éclairée réfléchi dans une direction donnée."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 Lois de réflexion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les lois de la réflexion seront utilisées dans le cas du présent travail et sont rappelées dans ce paragraphe.\n",
    "Lorsqu'un rayon lumineux (rayon incident) arrive sur une surface (point d'incidence), ce rayon peut subir une réflexion de différentes manières: soit de façon spéculaire, soit de façon diffuse (fig. 1). Ces réflexions dépendent de la rugosité de la surface (Lopez, 2007).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='lo_speculardiffuse.jpg' width=\"350\">\n",
    "Fig 1 : Réflexion spéculaire et Réflexion diffuse (Lopez, 2007)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***Réflexion Spéculaire***\n",
    "\n",
    "La réflexion spéculaire est un phénomène où le rayon suit les lois de Snell-Descartes. Le rayon incident forme un angle avec la normale du point incident, appelé angle d'incidence $\\theta_i$. De la même manière, l'angle formé entre le rayon réfléchi et la normale est l'angle réfléchi $\\theta_r$. La loi de Snell-Descartes démontre que les valeurs d'angles incident et réfléchi doivent être identiques comme présenté sur la figure 2. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='snell_descartes.PNG' width=\"250\">\n",
    "Fig 2 : Schéma de la diffusion spéculaire selon la loi de Snell-Descartes(Imen,2010)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Réflexion Diffuse***\n",
    "\n",
    "La réflexion diffuse, quant à elle, signifie que la surface réfléchi les photons du rayon incident aléatoirement dans toutes les directions du point d'incidence. Le résultat est une diffusion de la lumière. A l'inverse de la surface spéculaire où elle est parfaitement polie comme par exemple un miroir, la surface diffuse est rugueuse. \n",
    "\n",
    "A l'échelle microscopique, chaque photon contenu dans le rayon est intersecté sur un point de la surface et est réfléchi selon la loi de Snell-Descartes. En effet, la normale de ce point local (en vert) peut être d'orientation différente de la normale générale (en rouge) de la surface comme démontre la figure 3. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='diffuse.PNG' width=\"250\">\n",
    "Fig 3 : réflexions spéculaires des photons sur une surface rugueuse (Chaire CRSNG,2004)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Réflexion mixte***\n",
    "\n",
    "Ces deux modes de réflexion évoqués sont deux cas extrèmes et irréalistes. En effet, une surface n'est pas parfaitement diffuse, ni spéculaire, mais un compromis entre les deux (Lopez, 2007). \n",
    "\n",
    "A l'échelle mascropique, la lumière est un peu dispersée autour du vecteur de la réflexion spéculaire. Sa distribution correspond à celle du lobe bleu comme démontre la figure 4. \n",
    "La forme du lobe est fonction de la rugosité de la surface (Lopez,2007). Plus une surface est rugueuse, plus le lobe sera dispersé autour de l'axe de réflexion spéculaire. A l'inverse, plus une surface sera polie, plus le lobe sera fin et ressemblera à la réflexion spéculaire parfaite. \n",
    "\n",
    "Si nous émettons l'hypothèse que la surface est de nature isotrope, la forme du lobe sera symétrique autour de l'axe de réflexion spéculaire (Lopez, 2007). En effet, tout déplacement d'un certain angle $\\Omega$ de l'axe entraine une diminution de la luminance. Si, en gardant ce même angle d'écart $\\Omega$, on se déplace autour de l'axe, la luminance sera identique.\n",
    "\n",
    "La diffusion peut être modélisée comme la diffusion lambertienne (Lopez, 2007) où la lumière est éparpillée également dans toutes les directions de la surface. La luminance est identique dans toutes les directions. Cette diffusion est représentée par le demi-cercle vert de la figure 4. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='reflection.PNG' width=\"300\">\n",
    "Fig 4. Schéma de la réflexion mixte (Lopez, 2007)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2 Coefficients de Fresnel et Loi de Snell"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'avenir d'un photon lorsqu'il est à un point d'incidence entre deux milieux distincts est soit d'être réfléchi, soit d'être transmis dans le nouveau milieu  (Greve, 2004). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='reflection_transmission.PNG' width=\"300\">\n",
    "Fig 5. Trajectoire des photons: réflexion et transmission(Greve, 2004)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La loi de Snell permet de déterminer l'angle entre la normale et le photon transmis, appelé l'angle de transmission ($\\theta_t$) selon l'équation suivante:\n",
    "$$\n",
    "\\eta_2sin\\theta_t = \\eta_1sin\\theta_i   \\iff  \\theta_t=arcsin\\left(\\frac {\\eta_1} {\\eta_2}sin(\\theta_i)\\right)     \n",
    "$$\n",
    "Cet angle varie en fonction de l'angle d'incidence et des indices de réfraction des milieux.\n",
    "Il est intéressant de souligner que l'angle de transmission est limité à 90° de la normale du point incident. Au-delà de 90°, le vecteur de transmission du photon est considéré comme un vecteur de réflexion. Dans ce cas, la réflexion est totale , il n'y a pas de réfraction possible des photons. Ce phénomène se produit lorque l'angle incident est supérieur à une certaine valeur d'angle, appelée angle critique $\\theta_c$ (Perrote,2009)\n",
    "$$\n",
    "\\theta_c=arcsin\\left(\\frac {\\eta_2} {\\eta_1}sin(90)\\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les probabilités qu'un photon soit réfléchi ou transmis se calculent selon les équations de Fresnel (Greve, 2004). Ces équations dépendent de l'indice de réfraction de chacun des deux milieux ($\\eta_1$et $\\eta_2$) ainsi de l'angle d'incidence ($\\theta_i$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La probabilité qu'un photon soit réfléchi est formulée comme la moyenne de la composante horizontale et de la composante verticale du vecteur réfléchi.\n",
    "Dans le cas où il y a une réflexion totale, la probabilité vaut naturellement 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "R(\\theta_i) = \\begin{cases}\\frac{R_\\perp(\\theta_i) + R_\\lVert(\\theta_i)}{2} & \\text{if not total reflection} \\\\1 & \\text{if total reflection}\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La composante verticale et la composante horizontale du vecteur réfléchi sont calculées selon les équations suivantes:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "R_\\perp(\\theta_i) = \\left(\\frac {\\eta_1cos(\\theta_i) - \\eta_2cos(\\theta_t) } {\\eta_1cos(\\theta_i) + \\eta_2cos(\\theta_t)}\\right)^2\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "R_\\lVert(\\theta_i) = \\left(\\frac {\\eta_2cos(\\theta_i) - \\eta_1cos(\\theta_t) } {\\eta_1cos(\\theta_i) + \\eta_2cos(\\theta_t)}\\right)^2\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "avec $$ cos(\\theta_t) = \\sqrt{1-\\left(\\frac {\\eta_1} {\\eta_2}\\right)^2 sin(\\theta_i) ^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La probabilité qu'un photon soit transmis est déterminée par la différence entre l'unité et la probailité de réflexion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T(\\theta_i) = 1 - R(\\theta_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Script "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce chapitre est consacré à la méthodologie de la modélisation du ray-tracing d'un photon en un instant $t$ dans une chambre de culture avec trois lampes grâce au logiciel VTK couplé avec Python. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le code a été inspiré de *From Ray Casting to Ray Tracing with Python and VTK* sous le pseudo *somada141* (2014) représentant le ray tracing du soleil jusqu'à la terre et la réflexion spéculaire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le ray-tracing est le **processus** du programme. Un **pré-processus** est effectué pour mettre en place l'environnement virtuelle et ses composantes. Le paramètrage des propriétés de la lampe peut également être réalisé. Une ébauche est proposée dans le script *pre_processing.ipynb* rassemblant quelques notions de photométrie.  \n",
    "Un **post-processus** récolte les informations des photons intersectés par la plante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Données d'entrée\n",
    "\n",
    "Les données d'entrée nécessaires pour le cas présent sont les dimensions de la chambre, les indices de réfraction des différents milieux et des informations sur les lampes.\n",
    "\n",
    "Parmi les informations des lampes, se trouve l'angle d'ouverture (par défaut 180 pour chaque lampe), la localisation de chaque lampe, ainsi que les longueurs d'onde (choix arbitraire).\n",
    "\n",
    "Les indices de réfraction des parois et de l'air sont choisis suivant les données wikipédia (2020) et on considére que la paroi est de nature Polytéréphtalate d'éthylène (Mylar). \n",
    "En ce qui concerne la plante, elle absorbe ou réfléchi la lumière selon la longueur d'onde. Pour simplifier, nous considérons que les photons intersectés sur la plante suivront les probabilités selon les coefficients de Fresnel comme pour les parois. L'indice de réfraction pour la plante a été choisi aléatoirement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Importation des différentes librairies et du script \"définition\" contenant tous les fonctions utilisées. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le script *definition.py* est importé pour l'exécution du code. Le script *definition.pynb* contient les commentaires associés aux fonctions pour faciliter la compréhension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 427,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import sys\n",
    "import vtk \n",
    "import numpy\n",
    "import random\n",
    "from definition import *\n",
    "from datetime import datetime "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Fonctions qui permettent de convertir list/tuple en numpy et vice-versa. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les deux fonctions *lambda* permettent de convertir une liste ou un tuple en un tableau numérique et vice-versa. \n",
    "En effet, VTK utilise et génére des données telles que les coordonnées ou les vecteurs en object de type 'liste' ou 'tuple' qu'il faut convertir en objet de type 'numpy' pour effectuer des opérations vectorielles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2n = lambda l: numpy.array(l)\n",
    "n2l = lambda n: list(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Création du dictionnaire contenant les résultats & début du chrono de l'exécution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un dictionnaire _'plant_result'_ est créé. Il stockera les informations pour le post-processing. On enregistre aussi la date et l'heure actuelle pour calculer la durée de l'exécution du programme. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-01-02 19:10:58.676277\n"
     ]
    }
   ],
   "source": [
    "plant_result = {}\n",
    "startTime = datetime.now()\n",
    "print(startTime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Lecture fichier d'entrée .json"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le fichier d'entrée est de type json. Il fonctionne comme un dictionnaire mais permet d'être séparé du code. Cela est utile pour utiliser n'importe quel fichier de même structure dans un code. Les données des lampes et les informations de la chambre de culture sont inscrites dans ce fichier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 430,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "   with open('donnee_entree.json') as f:\n",
    "     data = json.load(f)\n",
    "except:\n",
    "   sys.exit(\"Error: enable to read file\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1 Pré-processus : Création de l'environnement virtuel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.1.1 Création de l'environnement visuel 3D**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction _vtkRenderer_ permet de créer l'environnement visuel où les différents objets seront ajoutés et le processus sera représenté. La couleur blanche pour le fond est choisie. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 431,
   "metadata": {},
   "outputs": [],
   "source": [
    "ren = vtk.vtkRenderer()\n",
    "ren.SetBackground(1, 1, 1) #(R,G,B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.1.2 La chambre de culture**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Création de la chambre"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La création et la visualisation d'un objet géometrique sont définies par la méthodologie suivante:\n",
    "$\n",
    "Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \\xrightarrow{\\text{}} Renderer\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Première étape :*\n",
    "\n",
    "$\n",
    "Source\n",
    "$\n",
    "\n",
    "$\n",
    "cube \n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La chambre de culture est représentée par un cube grâce à la fonction *vtkCubeSource* de la classe *Source*. La chambre de culture est conçue avec les dimensions indiquées dans le dictionnaire d'entrée Json grâce aux fonctions *SetXLength*, *SetYLength* et *SetZlength*. Ces longueurs sont référencées par rapport au centre de gravité du polygèdre. Ainsi, par exemple pour une hauteur de 30 cm, nous aurons -15 cm et +15 cm dans l'axe Y. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='axes_chambre.PNG' width=\"300\">\n",
    "Fig 6. Axes de repères avec la chambre de culture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 432,
   "metadata": {},
   "outputs": [],
   "source": [
    "xwidth = data[\"room\"][\"width\"]\n",
    "ywidth = data[\"room\"][\"height\"]\n",
    "zwidth = data[\"room\"][\"ceiling-height\"]\n",
    "cube = vtk.vtkCubeSource()\n",
    "cube.SetXLength(xwidth)\n",
    "cube.SetYLength(ywidth)\n",
    "cube.SetZLength(zwidth)\n",
    "cube.Update()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Deuxième étape:*\n",
    "\n",
    "$\n",
    "Source \\xrightarrow{\\text{}} Mapper \n",
    "$\n",
    "\n",
    "$\n",
    "cube \\xrightarrow{\\text{}} cubeMapper \n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La variable cube créée auparavant dans la classe *Source* est connectée dans la fonction *vtkPolyDataMapper* de la classe *Mapper* pour transformer la géométrie en primitives graphiques."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 433,
   "metadata": {},
   "outputs": [],
   "source": [
    "cubeMapper = vtk.vtkPolyDataMapper()\n",
    "cubeMapper.SetInputConnection(cube.GetOutputPort())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Troisième étape:*\n",
    "\n",
    "$\n",
    "Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \n",
    "$\n",
    "\n",
    "$\n",
    "cube \\xrightarrow{\\text{}} cubeMapper \\xrightarrow{\\text{}} cubeActor \n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les primitives graphiques sont transformées en entités géométriques grâce à la classe _actor_. La forme cube est officiellement créée et est prête à être visualisée grâce à la fonction _vtkActor_.\n",
    "La sous-classe *GetProperty* permet d'ajouter des options comme définir la position du centre de gravité (*SetOrigin*), la couleur (*setColor*) et l'opacité (*setOpacity*). Le centre de gravité est défini comme l'origine du repère de l'environnement virtuel. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "metadata": {},
   "outputs": [],
   "source": [
    "cubeActor = vtk.vtkActor()\n",
    "cubeActor.SetMapper(cubeMapper)\n",
    "cubeActor.SetOrigin(0, 0, 0) \n",
    "cubeActor.GetProperty().SetOpacity(0.1)\n",
    "cubeActor.GetProperty().SetColor(0, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Quatrième étape:*\n",
    "\n",
    "$\n",
    "Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \\xrightarrow{\\text{}} Renderer\n",
    "$\n",
    "\n",
    "$\n",
    "cube \\xrightarrow{\\text{}} cubeMapper \\xrightarrow{\\text{}} cubeActor \\xrightarrow{\\text{}} ren\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'acteur est ajouté à l'environnement virtuel _ren_ prédéfini auparavant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "metadata": {},
   "outputs": [],
   "source": [
    "ren.AddActor(cubeActor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- OBBTree "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction *oriented bounding box (OBB) Tree* est une des plus importantes du code. Elle va emboiter l'objet *cube*. Cet emboitement permet d'indiquer que les surfaces d'un objet sont des surfaces d'intersection. Autrement dit, cette fonction spécifie que ces surfaces sont les obstacles du ray-tracing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 436,
   "metadata": {},
   "outputs": [],
   "source": [
    "obbSurface = vtk.vtkOBBTree()\n",
    "obbSurface.SetDataSet(cube.GetOutput())\n",
    "obbSurface.BuildLocator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Coordonnées & Normales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On extrait les coordonnées des points et les normales des surfaces grâce à la fonction *cellcenter_normal* prédéfinie dans le script *definition*. Ces normales seront utilisées pour les lois des réflexions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 437,
   "metadata": {},
   "outputs": [],
   "source": [
    "pointsCellCentersSurface, normalsSurface = cellcenter_normal(cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.1.3 La plante**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Création de la plante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour simplifier dans un premier temps, une surface plane mime une plante à laquelle certains photons vont l'intersecter. La méthode est similaire à celle de la création de la chambre de culture $\n",
    "Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \\xrightarrow{\\text{}} Renderer\n",
    "$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction *Set Origin* permet de définir un premier sommet. Les fonctions *SetPoint1* et *SetPoint2* définissent deux autres sommets adjacents. La forme est le résultat d'un reliage entre ces trois points et d'un quatrième sommet créé automatiquement suite aux lois géométriques.\n",
    "Les fonctions *SetXResolution* et *SetYResolution* détermine l'échelle du maillage. La surface sera divisée en 15 x 15 cellules. Ces cellules permettront de classer les informations sur les intersections des photons avec la plante selon la celulle intersectée.\n",
    "Le maillage et la division de la surface en cellule est visible grâce à la fonction *EdgeVisibilityOn* de la sous classe *GetProperty* de la classe *Actor*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 438,
   "metadata": {},
   "outputs": [],
   "source": [
    "plant = vtk.vtkPlaneSource()\n",
    "plant.SetOrigin(-30, -75, -10)  #premier point\n",
    "plant.SetPoint1(-30, -75, 30)   #deuxième point\n",
    "plant.SetPoint2(40, -75, -10)   #troisième point\n",
    "plant.SetXResolution(15)\n",
    "plant.SetYResolution(15)\n",
    "plant.Update()\n",
    "\n",
    "plantMapper = vtk.vtkPolyDataMapper()\n",
    "plantMapper.SetInputConnection(plant.GetOutputPort())\n",
    "\n",
    "plantActor = vtk.vtkActor()\n",
    "plantActor.SetMapper(plantMapper)\n",
    "plantActor.GetProperty().SetColor(0, 1, 0)\n",
    "plantActor.GetProperty().EdgeVisibilityOn()  # show edges/wireframe\n",
    "plantActor.GetProperty().SetEdgeColor(0, 1, 0)\n",
    "\n",
    "ren.AddActor(plantActor)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- OBBTree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On indique que cette surface plane est une surface d'intersection grâce à la fonction *ObbTree*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 439,
   "metadata": {},
   "outputs": [],
   "source": [
    "obbPlant = vtk.vtkOBBTree()\n",
    "obbPlant.SetDataSet(plant.GetOutput())\n",
    "obbPlant.BuildLocator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Coordonnées & Normales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On extrait les coordonnées des points et les normales des surfaces de la plante grâce à la fonction *cellcenter_normal* prédéfinie dans le script *definition*. Ces normales seront utilisées pour les lois des réflexions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'cellcenter_normal' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-9265219edb6a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mpointsCellCentersplante\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnormalsplante\u001b[0m \u001b[1;33m=\u001b[0m\u001b[0mcellcenter_normal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mplant\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'cellcenter_normal' is not defined"
     ]
    }
   ],
   "source": [
    "pointsCellCentersplante,normalsplante =cellcenter_normal(plant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.1.4 Les lampes**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Création des lampes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La suite du code est une boucle itérative où chaque lampe est calculée individuellement. En effet, pour chaque lampe, la forme sera créée et s'en suit le processus du ray-tracing de cette lampe.\n",
    "Pour des raisons de clareté, la boucle est scindée en plusieurs paragraphes pour pouvoir de manière plus explicite, commenter les étapes. Le code de la boucle complet se trouve à la fin de ces explications et peut être exécuté."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En fonction de son angle d'ouverture, la lampe est représentée comme une portion de la sphére. Dans le cas présent, les lampes ont un angle d'ouverture de 180°, la forme sera donc une demi-sphére. Les fonctions *SetStartTheta* et *SetEndTheta* détermine l'angle de début (180°) et de fin pour la demi-sphère (360°).\n",
    "\n",
    "La méthode est similaire aux deux objets précédents ( $Source \\xrightarrow{\\text{}} Mapper \\xrightarrow{\\text{}} Actor \\xrightarrow{\\text{}} Renderer$). Il est nécessaire de rappeler que l'origine est située au centre du cube. La localisation de la lampe est référencée à cette origine, comme la localisation de la plante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce qui est vraiment important à noter ici, ce sont les résolutions thêta et phi de la demi-sphère réglée par les fonctions *SetThetaResolution* et *SetPhiResolution* respectivement. VTK crée la demi-sphère en utilisant des coordonnées sphériques et les paramètres de résolution définissent le nombre de points le long des directions de longitude (azimut) et de latitude (élevation) respectivement. Une résolution plus faible se traduira par moins de points et donc moins de cellules du maillage définissant la demi-sphère, ce qui lui donnera l'aspect d'un polygèdre grossier.\n",
    "\n",
    "De plus, la résolution de la lampe est intéressante car, en plus d'arrondir au mieux la demi-sphére, le maillage de la lampe va permettre le lancer des photons. C'est à partir des points centraux des cellules du maillage que les photons démarrent leurs trajectoires du ray-tracing. Plus la résolution est grande, plus le nombre de cellules du maillage est élevé, plus le nombre de photons lancés est important.\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "v = 0\n",
    "i = 1\n",
    "while v < 3:\n",
    " x = data[\"array\"][str(v)][\"x\"] #par rapport à l'origine qui est au centre du cube\n",
    " y = data[\"array\"][str(v)][\"y\"]\n",
    " z = data[\"array\"][str(v)][\"z\"]\n",
    "\n",
    " Lamp = vtk.vtkSphereSource()\n",
    " Lamp.SetCenter(x, y, z)\n",
    " Lamp.SetRadius(5)\n",
    " angle_ouverture = data[\"array\"][str(v)][\"angle\"]\n",
    " angle = 180 - angle_ouverture\n",
    " Lamp.SetStartTheta(180 + angle/2)\n",
    " Lamp.SetEndTheta(360 - angle/2)\n",
    " Lamp.SetPhiResolution(30)\n",
    " Lamp.SetThetaResolution(30)\n",
    "\n",
    " mapper = vtk.vtkPolyDataMapper()\n",
    " mapper.SetInputConnection(Lamp.GetOutputPort())\n",
    "\n",
    " pointActor = vtk.vtkActor()\n",
    " pointActor.SetMapper(mapper)\n",
    " pointActor.SetOrigin(0, 0, 0) \n",
    " pointActor.GetProperty().SetColor(color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    " pointActor.GetProperty().EdgeVisibilityOn()  # show edges/wireframe\n",
    " pointActor.GetProperty().SetEdgeColor(color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "\n",
    " ren.AddActor(pointActor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La figure 6 montre un aperçu de l'environnement créé avec la plante et les trois lampes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='axes.png' width=\"300\">\n",
    "Fig 6. Représentation de l'environnement virtuel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Coordonnées & Normales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les coordonnées des points pour chaque cellule du maillage de la lampe sont calculées et le point zéro qui correspond au centre de la cellule est extrait. On extrait également les normales de la lampe grâce à la fonction *cellcenter_normal* définie dans le script *definiton*."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    " pointsCellCentersLamp,normalsLamp =cellcenter_normal(Lamp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2 Processus : Ray-tracing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'étape suivante du code démarre le processus du ray-tracing et est une boucle. En effet, chaque photon et sa trajectoire, un à la fois, sont créés pour chaque cellule de la lampe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour chaque point central de la cellule du maillage de la lampe, ce point et la normale de la cellule sont extraits en variable temporaires. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    " for idx in range(pointsCellCentersLamp.GetNumberOfPoints()):\n",
    "   pointLamp = pointsCellCentersLamp.GetPoint(idx)\n",
    "   normalLamp = normalsLamp.GetTuple(idx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les coordonnées du point terminal de la trajectoire du photon sont calculées en prolongeant la normale d'un facteur de 500 à partir du point central de la cellule étudiée. Le choix du facteur 500 est justifié par la certitude que le photon atteint les parois ou la plante. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "   pointRayTarget = n2l(l2n(pointLamp) + 500*l2n(normalLamp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Grâce aux fonctions préalablement définies dans le fichier *définition*, la fonction *IntersectWithLine* qui est couplée avec la fonction *obbTree* permet de détecter les intersections des photons aux surfaces de la parois et à la plante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.2.1 Intersection avec la plante**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si une intersection du photon avec la plante est trouvée, on extrait l'indice de la cellule intersectée (point incident) et ses coordonnées. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    " if isHit(obbPlant, pointLamp, pointRayTarget):  # intersections were found\n",
    "        ptsInter, cIdsInter = GetIntersect(obbPlant, pointLamp, pointRayTarget) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La normale de la cellule est prélevée et le vecteur incident est mesuré. La variable X permet de déduire si le photon est réfléchi ou transmis suite à un tirage aléatoire selon les probabilités détermininées par les équations Fresnel. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "       normalplante = normalsplante.GetTuple(cIdsInter[0])\n",
    "       vecInc = n2l(l2n(ptsInter[0]) - l2n(pointLamp))\n",
    "       X = Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_plante\"],vecInc, normalplante)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si le photon est transmis dans la plante, par conséquent l'énergie lumineuse est absorbée dans la matière, il est intéressant de récolter diverses données pour le post-processing. La trajectoire du photon est également dessinée pour la visualisation grâce aux fonctions *addLine* et *addPoint*, définies dans le script *definition*. Dans le cas contraire, le photon est réfléchi et donc aucune énergie lumineuse n'est absorbée dans la plante.\n",
    "\n",
    "Des informations comme les coordonnées *ptsInter* du point incident et le numéro de la lampe originaire du photon *v* sont enregistrés dans le dictionnaire *plant_result* créé au préalable. Ces informations sont regroupées par cellule du maillage de la plante. Pour cela, une clé (*key*) avec l'indice de la cellule intersectée est créée. Si ce numéro a déjà été enregistré dans le dictionnaire, on encode les informations du photon intersecté *i* dans cette clé. Si ce numéro n'existait pas encore, on l'ajoute dans le dictionnaire."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "        if X[0] == \"T\" :\n",
    "          addLine(ren, pointLamp, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "          addPoint(ren, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "          \n",
    "          key = \"surface_\"+str(cIdsInter[0])\n",
    "          \n",
    "          if key in plant_result: \n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i] = ptsInter[0]\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "              i+=1\n",
    "              \n",
    "             \n",
    "          else:\n",
    "             plant_result.update({\"surface_\"+str(cIdsInter[0]):{\"coord\":{},\"lampe\":{}}})\n",
    "             plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i]= ptsInter[0]\n",
    "             plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "             i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.2.2 Intersection avec la paroi**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si une intersection du photon avec la paroi est trouvée, on extrait l'indice de la cellule et les coordonnées du point incident. La trajectoire du photon est dessinée pour la visualisation. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "   if isHit(obbSurface, pointLamp, pointRayTarget): \n",
    "        pointsInter, cellIdsInter = GetIntersect(obbSurface, pointLamp, pointRayTarget)\n",
    "        addLine(ren, pointLamp, pointsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]), opacity = 0.1)\n",
    "        addPoint(ren, pointsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]), opacity = 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comme dans le cas de l'intersection avec la plante, la normale de la cellule est prélevée et le vecteur incident de la trajectoire du photon est mesuré. La variable X permet de déduire si le photon est réfléchi ou transmis suite à un tirage aléatoire selon les probablités détermininées selon les équations Fresnel.  "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "        normalSurface = normalsSurface.GetTuple(cellIdsInter[0])\n",
    "        vecInc = n2l(l2n(pointsInter[0]) - l2n(pointLamp))\n",
    "        #vecRef = calcVecR(vecInc, normalSurface)\n",
    " \n",
    "        X =Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_paroi\"],vecInc, normalSurface)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si le photon est réfléchi, la direction de celui-ci est choisi aléatoirement puisque la rugosité de la surface est inconnue et donc la normale locale du point incident aussi.\n",
    "\n",
    "Le point d'incidence est réprésenté par une sphère. Un point central d'une cellule du maillage de cette sphère est tiré aléatoirement et la normale de cette cellule sera la trajectoire du photon réfléchi. Un nouveau point terminal de la trajectoire du photon réfléchi est déterminé. \n",
    "\n",
    "Dans le cas contraire, si le photon est transmis dans la paroi, par conséquent l'énergie lumineuse est absorbée dans la matière, il n'y a pas de réflexion de la lumière dans la chambre."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "        if X[0] == \"R\":\n",
    "            pointsCellCenterspoint, normalspoint =cellcenter_normal(point)\n",
    "            m = random.choice(range(pointsCellCenterspoint.GetNumberOfPoints()))\n",
    "            vecRef = normalspoint.GetTuple(m)\n",
    "            \n",
    "               \n",
    "            pointRayReflectedTarget = n2l(l2n(pointsInter[0]) + 500*l2n(vecRef))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mettons en évidence la faiblesse de cette partie du code. Le point incident est réprésenté par une sphère dont la surface de la paroi la traverse. Autrement dit, une demi-sphère est orientée à l'intérieur de la chambre et l'autre moitié est orientée à l'extérieur de la chambre. La normale choisie aléatoirement peut donc être une qui est dans le milieu extérieur. Dans ce cas, nous n'avons pas un photon réfléchi. Le problème peut être résolu en schématisant le point incident comme uniquement une demi-sphère orientée vers l'intérieure de la chambre. Cependant, l'orientation de la demi-sphère varie selon la paroi de la chambre de culture. En effet, pour la paroi de gauche, la demi-sphère devrait être orientée vers la gauche mais pour la paroi de droite, elle devrait être orientée vers la droite. A ce jour, une méthode pour éviter de produire un photon réfléchi dans le milieu extérieur n'a pas encore été trouvée. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.2.3 Deuxième réflexion depuis la paroi**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est possible de détecter l'intersection du photon réfléchi d'une paroi avec soit la plante soit une autre paroi. La méthode est similaire que le premier lancer d'un photon de la lampe. "
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "    try :          \n",
    "     if isHit(obbSurface,  pointsInter[0], pointRayReflectedTarget): \n",
    "          \n",
    "           pointsInter, cellIdsInter = GetIntersect(obbSurface,  pointsInter[0], pointRayReflectedTarget)\n",
    "           addLine(ren, pointsInter[0], pointsInter[1], color(data[\"array\"][str(v)][\"wavelength\"]),opacity=1)\n",
    "           addPoint(ren, pointsInter[1], color(data[\"array\"][str(v)][\"wavelength\"]), opacity =0.1)\n",
    "\n",
    "     if isHit(obbPlant, pointsInter[0], pointRayReflectedTarget):  \n",
    "           ptsInter, cIdsInter = GetIntersect(obbPlant,  pointsInter[0], pointRayReflectedTarget)\n",
    "           \n",
    "           normalplante = normalsplante.GetTuple(cIdsInter[0])\n",
    "           vecInc = n2l(l2n(ptsInter[0]) - l2n(pointsInter[0]))\n",
    "           \n",
    "           X =Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_plante\"],vecInc, normalSurface)\n",
    "        \n",
    "           if X[0] ==\"T\": \n",
    "               addLine(ren, pointsInter[0], ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "               addPoint(ren, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "       \n",
    "               key = \"surface_\"+str(cIdsInter[0])\n",
    "        \n",
    "               if key in plant_result: \n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i] = ptsInter[0]\n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"])\n",
    "                   i+=1\n",
    "            \n",
    "             \n",
    "               else:\n",
    "                   plant_result.update({\"surface_\"+str(cIdsInter[0]):{\"coord\":{},\"longueur d'onde\":{},\"lampe\":{}}})\n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i]= ptsInter[0]\n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "                   plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"])\n",
    "                   i+= 1\n",
    "     \n",
    "    except:\n",
    "       a = 0 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est intéressant de réfléchir à une boucle pour effectuer plusieurs réflexions succesives.\n",
    "Il est à nuancer aussi avec la limite de réflexion où la quantité de lumière devient négligeable. En effet, l'énergie lumineuse diminue avec la distance parcourue du photon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.2.4 Le code complet de la boucle**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 441,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 0\n",
    "i = 1\n",
    "while v < 3: # la boucle se termine à la fin du code!!\n",
    "\n",
    "  x = data[\"array\"][str(v)][\"x\"] #par rapport à l'origine qui est au centre du cube\n",
    "  y = data[\"array\"][str(v)][\"y\"]\n",
    "  z = data[\"array\"][str(v)][\"z\"]\n",
    "\n",
    "  Lamp = vtk.vtkSphereSource()\n",
    "  Lamp.SetCenter(x, y, z)\n",
    "  Lamp.SetRadius(5)\n",
    "  angle_ouverture = data[\"array\"][str(v)][\"angle\"]\n",
    "  angle = 180 - angle_ouverture\n",
    "  Lamp.SetStartTheta(180 + angle/2)\n",
    "  Lamp.SetEndTheta(360 - angle/2)\n",
    "  Lamp.SetPhiResolution(30)\n",
    "  Lamp.SetThetaResolution(30)\n",
    "\n",
    "  mapper = vtk.vtkPolyDataMapper()\n",
    "  mapper.SetInputConnection(Lamp.GetOutputPort())\n",
    "\n",
    "  pointActor = vtk.vtkActor()\n",
    "  pointActor.SetMapper(mapper)\n",
    "  pointActor.SetOrigin(0, 0, 0) \n",
    "  pointActor.GetProperty().SetColor(color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "  pointActor.GetProperty().EdgeVisibilityOn()  # show edges/wireframe\n",
    "  pointActor.GetProperty().SetEdgeColor(color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "\n",
    "  ren.AddActor(pointActor)\n",
    "\n",
    "  pointsCellCentersLamp,normalsLamp =cellcenter_normal(Lamp)\n",
    "\n",
    "\n",
    "# # #-----------------------------------------------------------------------------\n",
    "# # # RAY-TRACING\n",
    "# # #----------------------------------------------------------------------------\n",
    "\n",
    "\n",
    "  for idx in range(pointsCellCentersLamp.GetNumberOfPoints()):\n",
    "    \n",
    "    pointLamp = pointsCellCentersLamp.GetPoint(idx)\n",
    "    normalLamp = normalsLamp.GetTuple(idx)\n",
    "     \n",
    "    pointRayTarget = n2l(l2n(pointLamp) + 500*l2n(normalLamp))\n",
    "    \n",
    "    if isHit(obbSurface, pointLamp, pointRayTarget):  \n",
    " \n",
    "    \n",
    "        pointsInter, cellIdsInter = GetIntersect(obbSurface, pointLamp, pointRayTarget)\n",
    "        \n",
    "        addLine(ren, pointLamp, pointsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]), opacity = 0.1)\n",
    "        point = addPoint(ren, pointsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]), opacity = 0.1)\n",
    "\n",
    "        normalSurface = normalsSurface.GetTuple(cellIdsInter[0])\n",
    "        vecInc = n2l(l2n(pointRayTarget) - l2n(pointLamp))\n",
    "        \n",
    "        #vecRef = calcVecR(vecInc, normalSurface)\n",
    "         \n",
    "        X =Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_paroi\"],vecInc, normalSurface)\n",
    "        \n",
    "        if X[0] == \"R\":\n",
    "            pointsCellCenterspoint, normalspoint =cellcenter_normal(point)\n",
    "            m = random.choice(range(pointsCellCenterspoint.GetNumberOfPoints()))\n",
    "            vecRef = normalspoint.GetTuple(m)\n",
    "            while calcAngle(vecRef,normalSurface)> numpy.pi :\n",
    "                m = random.choice(range(pointsCellCenterspoint.GetNumberOfPoints()))\n",
    "                vecRef = normalspoint.GetTuple(m)\n",
    "            \n",
    "            pointRayReflectedTarget = n2l(l2n(pointsInter[0]) + 500*l2n(vecRef))\n",
    "           \n",
    "           \n",
    "           \n",
    "    \n",
    "    if isHit(obbPlant, pointLamp, pointRayTarget): \n",
    "        ptsInter, cIdsInter = GetIntersect(obbPlant, pointLamp, pointRayTarget)\n",
    "        \n",
    "        normalplante = normalsplante.GetTuple(cIdsInter[0])\n",
    "        vecInc = n2l(l2n(pointRayTarget) - l2n(pointLamp))\n",
    "        \n",
    "        X = Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_plante\"],vecInc, normalplante)\n",
    "        if X[0] == \"T\" :\n",
    "          addLine(ren, pointLamp, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "          addPoint(ren, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "          \n",
    "          key = \"surface_\"+str(cIdsInter[0])\n",
    "          \n",
    "          if key in plant_result: \n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i] = ptsInter[0]\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"]\n",
    "              i+=1\n",
    "              \n",
    "             \n",
    "          else:\n",
    "              plant_result.update({\"surface_\"+str(cIdsInter[0]):{\"coord\":{},\"longueur d'onde\":{},\"lampe\":{}}})\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i]= ptsInter[0]\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "              plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"]\n",
    "              i+=1\n",
    "             \n",
    "             \n",
    "    try :          \n",
    "      if isHit(obbSurface,  pointsInter[0], pointRayReflectedTarget):  \n",
    "          \n",
    "              pointsInter, cellIdsInter = GetIntersect(obbSurface,  pointsInter[0], pointRayReflectedTarget)\n",
    "     \n",
    "              addLine(ren, pointsInter[0], pointsInter[1], color(data[\"array\"][str(v)][\"wavelength\"]),opacity=0.1)\n",
    "              addPoint(ren, pointsInter[1], color(data[\"array\"][str(v)][\"wavelength\"]), opacity =0.1)\n",
    "         \n",
    "    \n",
    "      if isHit(obbPlant, pointsInter[0], pointRayReflectedTarget):  \n",
    "\n",
    "                ptsInter, cIdsInter = GetIntersect(obbPlant,  pointsInter[0], pointRayReflectedTarget)\n",
    "    \n",
    "                normalplante = normalsplante.GetTuple(cIdsInter[0])\n",
    "                vecInc = n2l(l2n(ptsInter[0]) - l2n(pointsInter[0]))\n",
    "        \n",
    "                X =Fresnel(data[\"room\"][\"n_air\"],data[\"room\"][\"n_plante\"],vecInc, normalSurface)\n",
    "                if X[0] ==\"T\": # s'il y a transmission des photons\n",
    "                    addLine(ren, pointsInter[0], ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "                    addPoint(ren, ptsInter[0], color(data[\"array\"][str(v)][\"wavelength\"]))\n",
    "                    \n",
    "                    key = \"surface_\"+str(cIdsInter[0])\n",
    "        \n",
    "                    if key in plant_result: \n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i] = ptsInter[0]\n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"]\n",
    "                      i+=1\n",
    "            \n",
    "             \n",
    "                    else:\n",
    "                      plant_result.update({\"surface_\"+str(cIdsInter[0]):{\"coord\":{},\"longueur d'onde\":{},\"lampe\":{}}})\n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])]['coord'][i]= ptsInter[0]\n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])]['lampe'][i] = v\n",
    "                      plant_result[\"surface_\"+str(cIdsInter[0])][\"longueur d'onde\"][i] = data[\"array\"][str(v)][\"wavelength\"]\n",
    "                      i+= 1\n",
    "      \n",
    "    except:\n",
    "        a = 0\n",
    "\n",
    "  v+=1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3 Post-Processus : Résultat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.3.1 Fichier de sortie**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les résultats du ray-tracing sont obtenus dans le dictionnaire _plant_result_. Ce dictionnaire est converti en un fichier externe .json afin de pouvoir l'utiliser pour un éventuel post-processing plus avancé. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 442,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"out.json\", \"w\") as f:\n",
    "        json.dump(plant_result, f,indent= 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le fichier est organisé par cellule de la plante. Dans chaque cellule, les coordonnées d'intersection, l'origine de la lampe et la longueur d'onde de la lampe pour chaque photon d'indice $i$ sont indiqués."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.3.2 Nombres de rayons**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 443,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chaque lampe a émis 1679 photons\n",
      "26 photons ont atteint la plante, soit 1.5485407980941037% des photons émis\n"
     ]
    }
   ],
   "source": [
    "print(\"Chaque lampe a émis \"+str(idx)+\" photons\")\n",
    "pourcentage = ((i-1)/idx)*100\n",
    "print(str(i-1)+\" photons ont atteint la plante, soit \"+str(pourcentage)+\"% des photons émis\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Résolution de la lampe|Nombre de photons émis| Nombre de photons transmis dans la plante| \n",
    "|----------------------|----------------------|------------------------------------------|\n",
    "|         10           |     477 =(3x159)     |               entre 0 et 10              |   \n",
    "|         20           |     2157 =(3x719)    |               entre 10 et 25             | \n",
    "|         30           |     5037 =(3x1679)   |               entre 25 et 35             |   \n",
    "|         50           |     14397 =(3x4799)  |               entre 70 et 110            |      \n",
    "|         100          |     58797 =(3x19599) |               quelques centaines         |    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.3.3 Durée de l'exécution**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Grâce aux dates enregistrées au début du code et à la fin, la durée de l'exécution est obtenue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 444,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed in 0:00:44.708494\n"
     ]
    }
   ],
   "source": [
    "endTime = datetime.now()\n",
    "deltaTime = endTime - startTime\n",
    "print(\"Completed in \"+str(deltaTime))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Comparaison des temps d'exécution en fonction du nombre de photons émis par une lampe._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|Résolution de la lampe|Nombre de photons|Temps d'exécution sur JupyterLab| Temps d'exécution sur Python Spyder |\n",
    "|----------------------|-----------------|--------------------------------|-------------------------------------|\n",
    "|         10           |  477 =(3x159)   |         0:00:02.352695         |            0:00:00.532576           |\n",
    "|         20           |  2157 =(3x719)  |         0:00:03.318503         |            0:00:03.057796           |\n",
    "|         30           |  5037 =(3x1679) |         0:00:07.688443         |            0:00:07.215807           |\n",
    "|         50           | 14397 =(3x4799) |         0:01:37.694579         |            0:01:05.575608           |\n",
    "|         100          | 58797 =(3x19599)|         0:02:44.450176         |            0:18:12.718293           |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sachant que le temps d'exécution varie en fonction du tirage sur base des probabilités de réflexion ou de transmission, le programme sur JupyterLab est capable d'exécuter avec un nombre de photon conséquent en un temps très court. Le temps d'exécution sur Python Spyder est supérieur lorsque le nombre de photon devient grand. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3.3.4 Visualisation**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est possible de visualiser la scène et le processus du ray-tracing. Cependant, il est possible d'exécuter ces lignes de code sur un ordinateur local (et non sur Jupyter NoteBook). Une configuration sur Jupyter NoteBook permettant cette action n'a pas encore été trouvée. Néanmoins, la figure 7 montre un aperçu du résultat."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src ='exemple_ray_tracing.PNG' width=\"300\">\n",
    "Fig 7. Visualisation du processus du ray-tracing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 445,
   "metadata": {},
   "outputs": [],
   "source": [
    "#renWin = vtk.vtkRenderWindow()\n",
    "#renWin.AddRenderer(ren)\n",
    "\n",
    "#iren = vtk.vtkRenderWindowInteractor()\n",
    "#iren.SetRenderWindow(renWin)\n",
    "\n",
    "#iren.Initialize()\n",
    "#renWin.Render()\n",
    "#iren.Start()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Conclusions et perspectives\n",
    "\n",
    "Le logiciel _Visualization ToolKit_ (VTK) permet une simulation du ray-tracing dans un espace de trois dimensions avec un haut niveau de précision. De plus, il a été démontré une vitesse d'exécution relativement rapide. \n",
    "\n",
    "Le code permet donc de modéliser l'intersection des photons venant d'une lampe caractérisée par sa longueur d'onde avec une surface photosynthétiquement active. Rappelons que les plantes absorbent la lumière selon la longueur d'onde. L'information de la longueur d'onde du photon intersecté mérite d'être exploitée en profondeur dans un autre processus afin de caractériser la réponse des plantes comme sa densité de flux de photons photosynthétiques (ppfd).\n",
    "\n",
    "Le fichier de sortie regroupe, par cellule du maillage de la plante, les photons transmis dans la plante et des informations pertinentes comme les coordonnées du point incident,la lampe origine du photon et la longueur d'onde de la lampe. Ce fichier peut être additionné avec les phénomènes et caractéristiques lumineux (intensité, luminance,...). Il peut être ensuite utilisé dans un autre processus pour caractériser la réponse des plantes. \n",
    "\n",
    "Le travail présenté est un point de départ vers d'autres voies d'amélioration et d'optimisation. Ci-dessous, une liste (non-exhaustive) de quelques perspectives concernant le ray-tracing: \n",
    "\n",
    "- Il est possible de texturer un objet dans VTK grâce à l'algorithme *vtkTexture*. Dès lors, une simplification de la réflexion en utilisant exclusivement la réflexion spéculaire est probable. Le temps de calcule se raccourcira et la précision des réflexions sera importante!\n",
    "- Les indices de réfractions $\\eta$ des milieux dépendent des longueurs d'ondes et de la température (Caye et al., 1968). Intégrer une fonction déterminant avec précision l'indice du milieu est indispensable pour déterminer au plus proche de la réalité, la proportion des photons réfléchis ou transmis.\n",
    "- Le code permet de procéder à un ray-tracing à un temps $t$ avec un seul photon lancé pour chaque cellule du maillage du code. Il est intéressant de trouver une méthode pouvant faire exécuter le code plusieurs fois, en fonction de la vitesse des ondes lumineuses, pour simuler une durée d'illumination $\\Delta t$ avec une infinité de lancer de photons.\n",
    "- Le code présent effectue une seule réflexion secondaire de la paroi. Il est envisageable d'effectuer une succession de réflexion du photon sur les parois. Effectuer des réflexions du photon sur la plante est aussi nécessaire pour compléter le processus.\n",
    "- Une coloration des cellules de la plante selon les longueurs d'onde des photons intersectés permettrait de préciser la visualisation et de présenter la situation réelle.\n",
    "- ... \n",
    "\n",
    "Les possibilités pour rendre le processus du ray-tracing sur VTK plus ressemblant de la réalité sont multiples et les études utilisant ce processus seront prometteuses!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5. Sources\n",
    "\n",
    "- Caye R. & Cervelle B., 1968. Détermination de l’indice de réfraction et du coefficient d’absorption des minéraux non transparents. *Bull. la Société française Minéralogie Cristallogr*. 91(3), 284–288.\n",
    "\n",
    "- Chaire CRSNG, 2004. Différente sortes de réflexion. *OPUS*,http://lamh.gmc.ulaval.ca/opus/physique534/resumes/12a.shtml, (02/01/2021).\n",
    "\n",
    "- Fiorillo D., Decaux B., Depta F. & Hu J., 2010. Outils de prototypage pour le dimensionnement de l’éclairement énergétique en milieux confinés. *Cah. Tech. Inra* 71, 5–40.\n",
    "\n",
    "- Greve B., 2004. Reflections and Refractions in Ray Tracing.\n",
    "\n",
    "- Henke M. & Buck-Sorlin G.H., 2017. Using a full spectral raytracer for calculating light microclimate in functional-structural plant modelling. *Comput. Informatics* 36(6), 1492–1522.\n",
    "\n",
    "- Imen E.,2010. Instrumentation optique pour la mesure des périkymaties de la couronne dentaire. Thèse de doctorat : Université de Franche-comté (France).\n",
    "\n",
    "- Kostro A., Huriet B. & Schüler A., 2007. PHOTONSIM: developing and testing a monte Carlo ray tracing software for the simulation of planar luminescent solar concentrators.\n",
    "\n",
    "- Lopez T., 2007. Rapport de Stage:Caractérisation de matériaux. Département informatique: Université de Rennes (France).\n",
    "\n",
    "- Perrote P.,2009. l'angle critique. http://siteheritage.csdecou.qc.ca/heritage2011/perrottep/anglecritique.htm, (31/12/2020).\n",
    "\n",
    "- Schroeder W.J., Avila L.S. & Hoffman W., 2000. Visualizing with VTK: A tutorial. *IEEE Comput. Graph. Appl.* 20(5), 20–27.\n",
    "\n",
    "- Somada141, 2014. From Ray Casting to Ray Tracing with Python and VTK. *PyScience*. https://pyscience.wordpress.com/2014/10/05/from-ray-casting-to-ray-tracing-with-python-and-vtk/, (31/12/20).\n",
    "\n",
    "- Tian F., 2016. Étude Et Optimisation Des Systèmes D’Éclairage Pour La Croissance Des Plantes En Milieu Contrôlé. Thèse de doctorat:Génie Electrique, Electronique et Télécommunications (GEET), Université de Toulouse (France).\n",
    "\n",
    "- VTK,2020.https://vtk.org/doc/nightly/html,(31/12/20).\n",
    "\n",
    "- Wikepedia,2020. Polytéréphtalate d'éthylène.https://fr.wikipedia.org/wiki/Polyt%C3%A9r%C3%A9phtalate_d%27%C3%A9thyl%C3%A8ne, (03/01/21)\n",
    "\n",
    "- Wikepedia,2020. Liste des indices de réfraction.https://fr.wikipedia.org/wiki/Liste_d%27indices_de_r%C3%A9fraction, (03/01/21)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}