{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Le lancer de rayons - 1ère partie\n",
        "\n",
        "# Illumination, ombres, réflexions\n",
        "\n",
        "Marc Lorenzi\n",
        "\n",
        "13 octobre 2018"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.pyplot as plt\n",
        "from math import *"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "plt.rcParams['figure.figsize'] = (15, 10)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Le __lancer de rayons__ (Ray Tracing en anglais) est une technique permettant le _rendu_ sur un écran d'une __scène__ composée d'__objets__.\n",
        "\n",
        "Dans toute notre discussion, nous travaillerons dans l'espace $\\mathbb R^3$ euclidien. Le repère canonique $(O,\\overrightarrow i, \\overrightarrow j, \\overrightarrow k)$ sera orienté un peu différemment de ce que l'on fait d'habitude : nous prendrons $\\overrightarrow i=(0, 1, 0)$, $\\overrightarrow j=(0, 0, 1)$ et $\\overrightarrow k=(0,0,-1)$. Pour faire simple, le plan $(O,\\overrightarrow i,\\overrightarrow j)$ est supposé \"vertical\"."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 1. Opérations sur les vecteurs et les points"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Nous aurons besoin d'effectuer toutes les opérations classiques sur les points et les vecteurs, que nous représenterons par des triplets de flottants : addition, soustraction, produit scalaire, distance, etc. Ces opérations se définissent sans difficulté."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 1.1 L'espace vectoriel $\\mathbb R^3$."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def add(u, v):\n",
        "    return (u[0] + v[0], u[1] + v[1], u[2] + v[2])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "add((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def sub(u, v):\n",
        "    return (u[0] - v[0], u[1] - v[1], u[2] - v[2])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "sub((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def oppose(u):\n",
        "    return (-u[0], -u[1], -u[2])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "oppose((1, 2, 3))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def prod(mu, u):\n",
        "    return (mu * u[0], mu * u[1], mu * u[2])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "prod(2, (1, 2, 3))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 1.3 Vecteur $\\overrightarrow {AB}$"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def vecteur(A, B):\n",
        "    return sub(B, A)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "vecteur((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 1.4 Produit scalaire, grandeurs associées"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def prodscal(u, v):\n",
        "    return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "prodscal((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def norme(u):\n",
        "    return sqrt(prodscal(u, u))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "norme((1, 2, 3))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def normaliser(u):\n",
        "    return prod(1 / norme(u), u)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "normaliser((1, 2, 3))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def norme2(u):\n",
        "    return prodscal(u, u)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "norme2((1, 2, 3))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def cosangle(u, v):\n",
        "    return prodscal(u, v) / (norme(u) * norme(v))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cosangle((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Et le sinus de l'angle, me direz-vous ? Nous en aurons besoin lorsque nous aurons à modéliser les lois de Snell-Descartes pour la réfraction. Nous le ferons le moment venu."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def distance(A, B):\n",
        "    return norme(vecteur(A, B))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "distance((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 1.5 Produit vectoriel"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def prodvect(u, v):\n",
        "    x = u[1] * v[2] - u[2] * v[1]\n",
        "    y = u[2] * v[0] - u[0] * v[2]\n",
        "    z = u[0] * v[1] - u[1] * v[0]\n",
        "    return (x, y, z) "
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "prodvect((1, 2, 3), (4, 5, 6))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 2. Caméras"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Pour observer une scène, il faut une __caméra__. La caméra est placée à une certaine position `pos`. Elle pointe dans une certaine direction `dirVue`. Ensuite, tout bon photographe fait en sorte que son appareil ne soit pas tout de travers. On associe à la caméra une base orthonormée directe $(\\overrightarrow u,\\overrightarrow v,\\overrightarrow w)$ telle que $\\overrightarrow w$ soit opposé à la direction de vue et $\\overrightarrow u$ soit horizontal. Enfin, il est souhaitable que la caméra ne soit pas à l'envers : le vecteur $\\overrightarrow v$ doit pointer vers le haut.\n",
        "\n",
        "Avoir une caméra c'est bien, mais il faut ensuite projeter le film. On projette sur un plan perpendiculaire à la direction de vue, à une distance $d$ de la caméra. L'écran a une certaine résolution (mesurée en pixels) et une certaine étendue mesurée, disons, en mètres : les bornes de l'écran.\n",
        "\n",
        "La classe `Camera` résume tout cela. Elle définit également une méthode `point`. Celle ci prend deux enrtiers $i$ et $j$ en paramètres, qui sont des coordonnées de pixels de l'écran. Elle renvoie le trpilet $(x, y, z)$ des coordonnées du point concerné dans le repère canonique."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Camera:\n",
        "    \n",
        "    def __init__(self, pos, dirVue, d, resolEcran, bornes):\n",
        "        self.nx, self.ny = resolEcran\n",
        "        self.xmin, self.xmax, self.ymin, self.ymax = bornes\n",
        "        self.pos = pos\n",
        "        self.dirVue = dirVue\n",
        "        self.w = normaliser(oppose(dirVue))\n",
        "        self.d = d\n",
        "        x, y, z = dirVue\n",
        "        self.u = normaliser((-self.w[2], 0, self.w[0]))\n",
        "        self.v = prodvect(self.w, self.u)\n",
        "        if self.v[1] < 0:\n",
        "            self.u = oppose(self.u)\n",
        "            self.v = oppose(self.v)\n",
        "        \n",
        "    def point(self, i, j):\n",
        "        x = self.xmin + i * (self.xmax - self.xmin) / self.nx\n",
        "        y = self.ymin + j * (self.ymax - self.ymin) / self.ny\n",
        "        z = -self.d\n",
        "        P = add(self.pos,add(prod(x, self.u),add(prod(y, self.v),prod(z, self.w))))\n",
        "        return P        "
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Voici la caméra que nous utiliserons de façon courante dans le notebook."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cameraSD = Camera(pos=(0, 0, 5), dirVue=(0, 0, -1), d = 5, resolEcran=(300, 200), bornes=(-3, 3, -2, 2))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cameraSD.u, cameraSD.v, cameraSD.w"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cameraSD.point(163, 189)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Je définis également une caméra HD, possédant un écran avec plus de pixels."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cameraHD = Camera(pos=(0, 0, 5), dirVue=(0, 0, -1), d = 5, resolEcran=(600, 400), bornes=(-3, 3, -2, 2))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "__Note sur les temps de calcul__ : au fur et à mesure que nous avancerons dans l'étude du lancer de rayons, nos algorithmes deviendront de plus en plus gourmands en temps de calcul. La caméra HD, par exemple, nécessite un temps de calcul 4 fois supérieur à la caméra classique. Plus la scène que nous voulons photographier possède d'obkets, plus elle possède de lumières, plus les calculs seront longs.\n",
        "\n",
        "Lorsque nous aborderons la réflexion, puis, enfin, la transparence, la complexité des algoriothmes va \"exploser\". La prise en compte de la transparence donne à l'algorithme de lancer de rayon une complexité ... exponentielle.\n",
        "\n",
        "__Moralité__ : restons modestes. Je travaillerai sur des scènes ayant 3 ou 4 objets, et une ou deux lumières. La caméra HD ne sera utilisée que rarement. Mais libre à vous, bien sûr, de faire tous les essais que vous voudrez !"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "cameraUHD = Camera(pos=(0, 0, 5), dirVue=(0, 0, -1), d = 5, resolEcran=(1200, 800), bornes=(-3, 3, -2, 2))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3. Objets et scènes"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Une scène est composée de plusieurs objets. Nous ne considérerons que des objets très simples : les __sphères__ et les __plans__. Rien ne vous empêche de considérer d'autres types d'objets, notre programme va être structuré pour qu'il soit très facile de rajouter de nouveaux types d'objets : cylindres, polygones, etc.\n",
        "\n",
        "Nous allons donc définir une classe __Sphere__ et une classe __Plan__. Une sphère est caractérisée par son centre et son rayon. Un plan, quant à lui, est défini par un point `org` et un vecteur normal `normal`."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.1 Sphères"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Sphere:\n",
        "    \n",
        "    def __init__(self, centre, rayon):\n",
        "        self.centre = centre\n",
        "        self.rayon = rayon "
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "s = Sphere(centre=(1,2,3), rayon=2)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "s.centre, s.rayon"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.2 Plans"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Plan:\n",
        "    \n",
        "    def __init__(self, org, normal):\n",
        "        self.org = org\n",
        "        self.normal = normal"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "p = Plan(org=(1, 1, 1), normal=(2, 0, 1))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "p.org, p.normal"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.3 Rayons"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Si notre algorithme du jour s'appelle le lancer de rayons, il va bien falloir se demander ce qu'est un rayon : c'est une demi-droite. Un rayon est donc défini par une origine `org` et une direction `u`. "
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Rayon:\n",
        "    \n",
        "    def __init__(self, org, u):\n",
        "        self.org = org\n",
        "        self.dir = u"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "r = Rayon(org=(1, 2, 3), u=(1, 1, 1))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "r.org, r.dir"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 3.4 Scène"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Définissions enfin la classe `Scene`. Une scène est composée d'objets ... et de lumières, parce que sinon on n'y voit rien :-). Pour l'instant, ne nous occupons pas des lumières, cela viendra en temps voulu."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Scene:\n",
        "    \n",
        "    def __init__(self, objets, lumieres):\n",
        "        self.objets = objets\n",
        "        self.lumieres = lumieres"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 4. Calculs d'intersections"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Il est temps d'expliquer le principe du lancer de rayon. Nous disposons d'une caméra avec un écran, et d'une scène. Pour chaque pixel de l'écran, la caméra lance un rayon qui passe par ce pixel. Le rayon intersecte la scène en un ou plusieurs objets. On calcule alors la valeur de l'éclairement pour ce pixel en fonction des points des objets frappés par le rayon.\n",
        "\n",
        "Bien entendu, tout est dans le \"On calcule alors la valeur de l'éclairement\". Le reste du notebook et le notebook suivant sont là pour mettre en place des modèles permettant ce calcul. Mais pour l'instant, occupons nous d'intersections !"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 4.1 Intersecter un rayon avec une sphere"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Soit $r$ un rayon d'origine le point $A$, dirigé par un vecteur $\\overrightarrow u$. Un point quelconque du rayon est de la forme $M = A + \\mu \\overrightarrow u$, où $\\mu$ est un réel positif ou nul.\n",
        "\n",
        "Soit $S$ la sphère de centre $B$ et de rayon $R$. Le point $M$ est sur la sphère si et seulement si $||\\overrightarrow{BM}||=R$, ce qui s'écrit encore $||\\overrightarrow{BA}+ \\mu \\overrightarrow u||=R$. Élevons au carré pour obtenir $<\\overrightarrow{BA}+ \\mu \\overrightarrow u,\\overrightarrow{BA}+ \\mu \\overrightarrow u> = R^2$, où les crochets désignent le produit scalaire. On développe :\n",
        "\n",
        "$$||\\overrightarrow u||^2 \\mu^2 + 2<\\overrightarrow{BA},\\overrightarrow u>\\mu+||\\overrightarrow{BA}||^2-R^2=0$$\n",
        "\n",
        "On obtient donc une équation du second degré en $\\mu$ du type $a\\mu^2+b\\mu+c=0$, où $a=||\\overrightarrow u||^2$, $b=2<\\overrightarrow{BA},\\overrightarrow u>$ et $c=||\\overrightarrow{BA}||^2-R^2$. Selon le signe du discriminant de cette équation, le rayon a 0, 1 ou deux points d'intersection avec la sphère.\n",
        "\n",
        "Nous redéfinissons donc la classe `Sphere`pour qu'une sphère soit capable de calculer son intersection avec un rayon. La méthode `intersection`renvoie un couple formé de la sphère et de la liste, éventuellement vide, des points d'intersection. La méthode reprend point par point le cheminement du calcul ci-dessus.\n",
        "\n",
        "J'en profite pour rajouter une méthode `vecteur_normal` qui calcule un vecteur normal à une sphère en un point (ceci est immédiat)."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Sphere:\n",
        "    \n",
        "    def __init__(self, centre, rayon):\n",
        "        self.centre = centre\n",
        "        self.rayon = rayon\n",
        "        \n",
        "    def intersection(self, rayon):\n",
        "        A, u = rayon.org, rayon.dir\n",
        "        B, R = self.centre, self.rayon\n",
        "        v = vecteur(B, A)\n",
        "        a = norme(u) ** 2\n",
        "        b = 2 * prodscal(u, v)\n",
        "        c = norme(v) ** 2 - R ** 2\n",
        "        Delta = b ** 2 - 4 * a * c\n",
        "        if Delta < 0: return []\n",
        "        elif Delta == 0:\n",
        "            mu = - b / (2 * a)\n",
        "            if mu > 1e-5:\n",
        "                return [(self, add(A, prod(mu, u)))]\n",
        "            else: return 0\n",
        "        elif Delta > 0:\n",
        "            s = []\n",
        "            delta = sqrt(Delta)\n",
        "            mu1 = (-b - delta) / (2 * a)\n",
        "            if mu1 > 1e-5:\n",
        "                s.append((self, add(A, prod(mu1, u))))\n",
        "            mu2 = (-b + delta) / (2 * a)\n",
        "            if mu2 > 1e-5:\n",
        "                s.append((self, add(A, prod(mu2, u))))\n",
        "            return s\n",
        "            \n",
        "    def vecteur_normal(self, A):\n",
        "        return vecteur(self.centre, A)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "s = Sphere((1, 0, -3), -2)\n",
        "s.intersection(Rayon((0, 0, 5), (0, 0, -1)))"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 4.2 Intersecter un rayon avec un plan"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Soit $r$ un rayon d'origine le point $A$, dirigé par un vecteur $\\overrightarrow u$. Un point quelconque du rayon est de la forme $M = A + \\mu \\overrightarrow u$, où $\\mu$ est un réel positif ou nul.\n",
        "\n",
        "Soit $P$ le plan passant par $B$ et de vecteur normal $\\overrightarrow n$. Le point $M$ est dans le plan $P$ si et seulement si $<\\overrightarrow{BM}, \\overrightarrow n> = 0$, ou encore $<\\overrightarrow{BA}+ \\mu \\overrightarrow u,\\overrightarrow n>=0$. On développe pour obtenir une équation de degré 1 en $\\mu$ :\n",
        "\n",
        "$$<\\overrightarrow{BA},\\overrightarrow n>+ <\\overrightarrow u,\\overrightarrow n>\\mu=0$$\n",
        "\n",
        "Cette équation possède en général une unique solution sauf, bien évidemment, si les vecteur $\\overrightarrow u$ et $\\overrightarrow n$ sont orthogonaux, puisque dans ce cas le rayon est parallèle au plan.\n",
        "\n",
        "Comme pour la classe `Sphere`, redéfinissons la classe `Plan`."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Plan:\n",
        "    \n",
        "    def __init__(self,org, normal):\n",
        "        self.org = org\n",
        "        self.normal = normal\n",
        "        \n",
        "    def intersection(self, rayon):\n",
        "        A, u = rayon.org, rayon.dir\n",
        "        B, n = self.org, self.normal\n",
        "        v = vecteur(B, A)\n",
        "        a = prodscal(u, n)\n",
        "        b = prodscal(v, n)\n",
        "        if a == 0: return []\n",
        "        else:\n",
        "            mu = - b / a\n",
        "            if mu > 1e-5:\n",
        "                return [(self, add(A, prod(mu, u)))]\n",
        "            else:\n",
        "                return []\n",
        "        \n",
        "    def vecteur_normal(self, A):\n",
        "        return self.normal"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 4.3 Premiers tests"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Créons un écran (une matrice de pixels) et une scène."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def ecran(m, n):\n",
        "    t = m * [None]\n",
        "    for k in range(m):\n",
        "        t[k] = n * [(0.3, 0.3, 1.)]\n",
        "    return t"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "ecran(5, 3)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Notre scène comporte trois sphères posées sur un plan."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "sphere1 = Sphere(centre=(2, 0, -5), rayon=2)\n",
        "sphere2 = Sphere(centre=(-1, -1, -3), rayon=1)\n",
        "sphere3 = Sphere(centre=(0, 1, -20), rayon=3)\n",
        "plan = Plan(org=(0, -2, 0), normal=(0, 1, 0))\n",
        "scene = Scene([sphere1, sphere2, sphere3, plan], [])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "La fonction `lancer0` lance des rayons depuis la caméra vers des pixels de l'écran. Lorsqu'un rayon rencontre un objet de la scène, on allume le pixel correspondant en blanc.\n",
        "\n",
        "Deux mots sur les couleurs : une couleur est représentée par un triplet $(r, g, b)$ où  $r,g,b\\in[0,1]$ représentent des quantités de rouge, de vert et de bleu. Ainsi, $(0,0,0)$ représente la couleur noire, $(1,1,1)$ la couleur blanche, et, par exemple, $(1,0,0)$ la couleur rouge."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancer0(scene, camera):\n",
        "    t = ecran(camera.ny, camera.nx)\n",
        "    for i in range(camera.nx):\n",
        "        for j in range(camera.ny):\n",
        "            M = camera.point(i, j)\n",
        "            r = Rayon(camera.pos, vecteur(camera.pos, M))\n",
        "            s = []\n",
        "            for objet in scene.objets:\n",
        "                s += objet.intersection(r)\n",
        "            if s != []: t[camera.ny - 1 - j][i] = (1., 1., 1.)\n",
        "    plt.axis('off')\n",
        "    plt.imshow(t)\n",
        "    plt.show()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer0(scene, cameraSD)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Évidemment, tout ceci est pour l'instant _très_ décevant. C'est tout à fait normal, puisque nous n'avons pas encore mis en place de modèle d'éclairement. Nous allons faire cela progressivement."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 5. Illumination"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 5.1 La classe `Lumiere`"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Qu'est-ce qu'une source de lumière ? Nous nous considérerons ici que des sources ponctuelles. Une telle source est caractérisée par sa position, et par l'intensité de rouge, de vert, et de bleu qu'elle émet (sa couleur, si vous préférez). D'où la classe `Lumiere`."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Lumiere:\n",
        "    \n",
        "    def __init__(self, org, I):\n",
        "        self.org = org\n",
        "        self.I = I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 5.2 Éclairement diffus"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Lorsque nous regardons un objet, nous percevons la lumière émise par cet objet.\n",
        "\n",
        "Un objet qui n'est pas une source de lumière est visible parce qu'il réfléchit la lumière qui lui est envoyée par d'autres sources. Il existe divers types de réflexion. Le premier type est la __réflexion diffuse__, étudiée par Lambert au dix-huitième siècle : la quantité d'énergie reçue par une surface dépend de l'angle d'arrivée de la lumière sur cette surface. Si la lumière arrive perpendiculairement à la surface, cette énergie est maximale. Si la lumière arrive tangentiellement, l'énergie reçue est nulle. Ceci mène au __modèle d'éclairement de Lambert__ :\n",
        "\n",
        "$$L = k_dI\\max(0, <\\overrightarrow u, \\overrightarrow v>)$$\n",
        "\n",
        "où\n",
        "\n",
        "- $k_d$ est le coefficient diffus, propre à l'objet illuminé. Plus $k_d$ est grand, plus l'objet réfléchit la lumière.\n",
        "- $I$ est l'intensité de la source lumineuse (en fait un triplet $(r, g, b)$ d'intensités)\n",
        "- $\\overrightarrow u$ est le vecteur unitaire normal à la surface\n",
        "- $\\overrightarrow v$ est le vecteur unitaire dirigé vers la source lumineuse.\n",
        "\n",
        "Pourquoi le max ? Un produit scalaire négatif correspond à une lumière arrivant de l'arrière (l'angle est supérieur à $\\frac\\pi 2$). L'objet n'est alors pas éclairé."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Laissons tomber pour l'instant le coefficient $k_d$ (prenons-le égal à 1) et ruons nous sur un test."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def diffuse(A, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        v = vecteur(A, L.org)\n",
        "        u = objet.vecteur_normal(A)\n",
        "        I = add(I, prod(max(0, cosangle(u, v)), L.I))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Prenons une scène composée de trois sphères, d'un plan et de deux sources lumineuses blanches (changez les paramètres, testez, testez !)."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "sphere1 = Sphere(centre=(2, 0, -5), rayon=2)\n",
        "sphere2 = Sphere(centre=(-1, -1, -3), rayon=1)\n",
        "sphere3 = Sphere(centre=(0, 1, -20), rayon=3)\n",
        "plan = Plan(org=(0, -2, 0), normal=(0, 1, 0))\n",
        "lum1 = Lumiere(org=(-100, 20, 100), I=(0.7, 0.7, 0.7))\n",
        "lum2 = Lumiere(org=(100, 100, 100), I=(0.7, 0.7, 0.7))\n",
        "scene = Scene([sphere1, sphere2, sphere3, plan], [lum1, lum2])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Réécrivons la fonction de lancer de rayons. Pour chaque rayon lancé de la caméra, on sélectionne, s'il y en a un, l'objet le plus proche frappé par le rayon. On calcule l'éclairement du point frappé au moyen du modèle de Lambert et on allume le pixel correspondant de l'écran à la couleur adéquate."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancer1(scene, camera):\n",
        "    t = ecran(camera.ny, camera.nx)\n",
        "    for i in range(camera.nx):\n",
        "        for j in range(camera.ny):\n",
        "            M = camera.point(i, j)\n",
        "            r = Rayon(camera.pos, vecteur(camera.pos, M))\n",
        "            s = []\n",
        "            for objet in scene.objets:\n",
        "                s += objet.intersection(r)\n",
        "            if s != []:\n",
        "                (objet, A) = plusProche(camera.pos, s)\n",
        "                Ix, Iy, Iz = diffuse(A, objet, scene)\n",
        "                t[camera.ny - 1 - j][i] = (min(Ix, 1), min(Iy, 1), min(Iz, 1))\n",
        "    plt.axis('off')\n",
        "    plt.imshow(t)\n",
        "    plt.show()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "La fonction `plusProche` renvoie, parmi une liste de couples (objet, point), le couple dont le point est le plus proche d'un point $B$ donné."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def plusProche(B, objets):\n",
        "    obj0, A0 = objets[0]\n",
        "    for obj, A in objets:\n",
        "        if distance(B, A) < distance(B, A0):\n",
        "            obj0 = obj\n",
        "            A0 = A\n",
        "    return obj0, A0"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer1(scene, cameraSD)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Ah oui, c'est quand même bien mieux que notre première fonction de lancer !"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 5.3 Éclairement spéculaire"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "L'éclairement diffus est indépendant de la position de la caméra. Cependant, beaucoup d'objets réels brillent plus ou moins et produisent des __réflexions spéculaires__ qui, elles, dépendent de l'endroit d'où vous regardez l'objet. Alors que le modèle de Lambert est adapté aux surfaces mates (le caoutchouc par exemple), le modèle que nous allons considérer est adapté aux surfaces qui brillent (l'acier poli par exemple). C'est le __modèle de Blinn-Phong__.\n",
        "\n",
        "$$ L = k_s I \\max(0, <\\overrightarrow v, \\overrightarrow h>)^p$$\n",
        "\n",
        "où\n",
        "\n",
        "- $k_s$ est un coefficient propre à l'objet, le coefficient spéculaire\n",
        "- $I$ est l'intensité de la source lumineuse\n",
        "- $\\overrightarrow v$ est le vecteur unitaire dirigé vers la source lumineuse\n",
        "- $p$ est un paramètre de brillance\n",
        "\n",
        "Quant au vecteur $\\overrightarrow h$, il est défini comme suit : soit $\\overrightarrow u$ le vecteur normal à la surface. Soit $\\overrightarrow w$ le vecteur unitaire dirigé vers l'observateur. Alors\n",
        "\n",
        "$$\\overrightarrow h=\\frac{\\overrightarrow u + \\overrightarrow  w}{||\\overrightarrow u + \\overrightarrow  w||}$$"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Quelles valeurs prendre pour $p$ ? Plus $p$ est grand plus ça brille. Voici quelques valeurs typiques de $p$ :\n",
        "\n",
        "- 10 : coquille d'oeuf\n",
        "- 100 : légèrement brillant\n",
        "- 1000 : vraiment brillant\n",
        "- 10000 : quasi-miroir\n",
        "\n",
        "Je prendrai 100 dans ce qui suit."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def phong(A, obs, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        v = vecteur(A, L.org)\n",
        "        v = prod(1 / norme(v), v)\n",
        "        u = objet.vecteur_normal(A)\n",
        "        w = vecteur(A, obs)\n",
        "        w = prod(1 / norme(w), w)\n",
        "        h = add(v, w)\n",
        "        I = add(I, prod(max(0, cosangle(u, h)) ** 100, L.I))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancer2(scene, camera):\n",
        "    t = ecran(camera.ny, camera.nx)\n",
        "    for i in range(camera.nx):\n",
        "        for j in range(camera.ny):\n",
        "            M = camera.point(i, j)\n",
        "            r = Rayon(camera.pos, vecteur(camera.pos, M))\n",
        "            s = []\n",
        "            for objet in scene.objets:\n",
        "                s += objet.intersection(r)\n",
        "            if s != []:\n",
        "                (objet, A) = plusProche(camera.pos, s)\n",
        "                Ix, Iy, Iz = add(diffuse(A, objet, scene), phong(A, camera.pos, objet, scene))\n",
        "                t[camera.ny - 1 - j][i] = (min(Ix, 1), min(Iy, 1), min(Iz, 1))\n",
        "    plt.axis('off')\n",
        "    plt.imshow(t)\n",
        "    plt.show()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer2(scene, cameraSD)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "__Remarque__ : Je profite du dessin qui vient d'apparaître pour parler d'un phénomène dont je ne discuterai plus dans la suite ...\n",
        "\n",
        "Les bords des sphères sont crénelés. C'est le phénomène d'__aliasage__. Ceci peut être résolu en découpant chaque pixel en 4 (par exemple) et en lançant des rayons sur chacun des quarts de pixels puis en faisant les moyennes des éclairements obtenus. Cette méthode d'anti-aliasage multiplie évidemment le temps de calcul par 4."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 5.4 Éclairement ambiant"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Il existe enfin un dernier type d'éclairement : __l'éclairement ambiant__. Une surface qui ne reçoit aucune lumière est normalement complètement noire. Dans la réalité ce n'est pas le cas, les réflexions indirectes d'autres surfaces (entre autres) produisent un léger éclairement. Cette \"lumière ambiante\" peut être ajoutée à notre modèle d'éclairement :\n",
        "\n",
        "$$L = k_aI_a$$\n",
        "\n",
        "où\n",
        "\n",
        "- $k_a$ est un coefficient ambiant propre à l'objet\n",
        "- $I_a$ est l'intensité de la lumière ambiante causée par la source lumineuse.\n",
        "\n",
        "Remarquons que nous avons maintenant trois paramètres, $k_a$, $k_d$ et $k_s$, propres à chaque objet, et deux paramètres, $I_a$ et $I$ propres à chaque source lumineuse. Il est temps de réécrire nos classes `Sphere`, `Plan` et `Lumiere` pour tenir compte de ces nouveaux paramètres."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Sphere:\n",
        "    \n",
        "    def __init__(self, centre, rayon, color, ka, kd, ks):\n",
        "        self.centre = centre\n",
        "        self.rayon = rayon\n",
        "        self.color = color\n",
        "        self.ka = ka\n",
        "        self.kd = kd\n",
        "        self.ks = ks\n",
        "        \n",
        "    def intersection(self, rayon):\n",
        "        A, u = rayon.org, rayon.dir\n",
        "        B, R = self.centre, self.rayon\n",
        "        v = vecteur(B, A)\n",
        "        a = norme(u) ** 2\n",
        "        b = 2 * prodscal(u, v)\n",
        "        c = norme(v) ** 2 - R ** 2\n",
        "        Delta = b ** 2 - 4 * a * c\n",
        "        if Delta < 0: return []\n",
        "        elif Delta == 0:\n",
        "            mu = - b / (2 * a)\n",
        "            if mu > 1e-5:\n",
        "                return [(self, add(A, prod(mu, u)))]\n",
        "            else:\n",
        "                return []\n",
        "        elif Delta > 0:\n",
        "            delta = sqrt(Delta)\n",
        "            mu1 = (-b - delta) / (2 * a)\n",
        "            mu2 = (-b + delta) / (2 * a)\n",
        "            s = []\n",
        "            if mu1 > 1e-5: s.append((self, add(A, prod(mu1, u))))\n",
        "            if mu2 > 1e-5: s.append((self, add(A, prod(mu2, u))))\n",
        "            return s\n",
        "        \n",
        "    def vecteur_normal(self, A):\n",
        "        return vecteur(self.centre, A)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Plan:\n",
        "    \n",
        "    def __init__(self, org, normal, color, ka, kd, ks):\n",
        "        self.org = org\n",
        "        self.normal = normal\n",
        "        self.color = color\n",
        "        self.ka = ka\n",
        "        self.kd = kd\n",
        "        self.ks = ks\n",
        "        \n",
        "    def intersection(self, rayon):\n",
        "        A, u = rayon.org, rayon.dir\n",
        "        B, n = self.org, self.normal\n",
        "        v = vecteur(B, A)\n",
        "        a = prodscal(u, n)\n",
        "        b = prodscal(v, n)\n",
        "        if a == 0: return []\n",
        "        else:\n",
        "            mu = - b / a\n",
        "            if mu > 1e-5:\n",
        "                M = add(A, prod(mu, u))\n",
        "                return [(self, M)]\n",
        "            else:\n",
        "                return []\n",
        "        \n",
        "    def vecteur_normal(self, A):\n",
        "        return self.normal"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "class Lumiere:\n",
        "    \n",
        "    def __init__(self, org, I, Ia):\n",
        "        self.org = org\n",
        "        self.I = I\n",
        "        self.Ia = Ia"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Les modèles déclairement doivent eux-aussi être réécrits pour tenir compte des nouvelles caractéristiques des objets et des lumières."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def diffuse(A, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        v = vecteur(A, L.org)\n",
        "        u = objet.vecteur_normal(A)\n",
        "        I1 = prod(objet.kd, L.I)\n",
        "        I = add(I, prod(max(0, cosangle(u, v)), I1))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def phong(A, obs, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        v = vecteur(A, L.org)\n",
        "        v = prod(1 / norme(v), v)\n",
        "        u = objet.vecteur_normal(A)\n",
        "        w = vecteur(A, obs)\n",
        "        w = prod(1 / norme(w), w)\n",
        "        h = add(v, w)\n",
        "        I1 = prod(objet.ks, L.I)\n",
        "        I = add(I, prod(max(0, cosangle(u, h)) ** 100, I1))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def ambient(A, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        I = add(I, prod(objet.ka, L.Ia))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "La fonction `illumination` additionne les trois modèles.\n",
        "\n",
        "__Remarque__ : si nos sources lumineuses sont trop intenses, des phénomènes de saturation peuvent apparaître dans le rendu. Une intensité lumineuse strictement supérieure à 1 est ramenée à 1 : on a un phénomène de __surexposition__."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def illumination(A, obs, objet, scene):\n",
        "    I = ambient(A, objet, scene)\n",
        "    I = add(I, diffuse(A, objet, scene))\n",
        "    I = add(I, phong(A, obs, objet, scene))\n",
        "    Ix, Iy, Iz = I\n",
        "    cx, cy, cz = objet.color\n",
        "    return (Ix * cx, Iy * cy, Iz * cz)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Maintenant, notre nouvelle fonction de lancer. L'occasion est trop belle, je scinde la fonction en deux parties pour la rendre plus lisible (et préparer le terrain pour la suite :-))."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancer3(scene, camera):\n",
        "    t = ecran(camera.ny, camera.nx)\n",
        "    for i in range(camera.nx):\n",
        "        for j in range(camera.ny):\n",
        "            M = camera.point(i, j)\n",
        "            r = Rayon(camera.pos, vecteur(camera.pos, M))\n",
        "            Ix, Iy, Iz = lancerRayon(r, camera)\n",
        "            t[camera.ny - 1 - j][i] = (min(Ix, 1), min(Iy, 1), min(Iz, 1))\n",
        "    plt.axis('off')\n",
        "    plt.imshow(t)\n",
        "    plt.show()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "La fonction de lancer est donc extrêmenent courte et lisible. Reste à écrire `lancerRayon`."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancerRayon(r, camera):\n",
        "    s = []\n",
        "    for objet in scene.objets:\n",
        "        s += objet.intersection(r)\n",
        "    if s != []:\n",
        "        (objet, A) = plusProche(camera.pos, s)\n",
        "        return illumination(A, camera.pos, objet, scene)\n",
        "    else: return (0.5, 0.5, 1.)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Voici notre scène."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "sphere1 = Sphere(centre=(2, -1, -5), rayon=1, color=(0.8, 0.8, 0.8), ka=0.1, kd=1, ks=1)\n",
        "sphere2 = Sphere(centre=(-1, -1, -3), rayon=1, color=(0.8, 0.3, 0.3), ka=0.1, kd=1, ks=0.)\n",
        "sphere3 = Sphere(centre=(0, 1, -12), rayon=3, color=(1, 0.8, 0.4), ka=0.1, kd=1, ks=1)\n",
        "plan = Plan(org=(0, -2, 0), normal=(0, 1, 0), color=(0.9, 0.6, 0.6), ka=0.1, kd=1, ks=0)\n",
        "lum1 = Lumiere(org=(60, 50, 100), I=(1, 1, 1), Ia=(1, 1, 1))\n",
        "scene = Scene([sphere1, sphere2, sphere3, plan], [lum1])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Lançons des rayons !"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer3(scene, cameraSD)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Il y a un petit souci : les sphères sont soit-disant posées sur le plan mais on a l'impression qu'elles flottent dans les airs. Que manque-t-il ? Il manque les __ombres__. "
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 6. Ombres"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def diffuse(A, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        r = Rayon(A, L.org)\n",
        "        black = False\n",
        "        for obj in scene.objets:\n",
        "            if obj.intersection(r) != []:\n",
        "                black = True\n",
        "                break\n",
        "        if not black:\n",
        "            v = vecteur(A, L.org)\n",
        "            u = objet.vecteur_normal(A)\n",
        "            I1 = prod(objet.kd, L.I)\n",
        "            I = add(I, prod(max(0, cosangle(u, v)), I1))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def phong(A, obs, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        r = Rayon(A, L.org)\n",
        "        black = False\n",
        "        for obj in scene.objets:\n",
        "            if obj.intersection(r) != []:\n",
        "                black = True\n",
        "                break\n",
        "        if not black:\n",
        "            v = vecteur(A, L.org)\n",
        "            v = prod(1 / norme(v), v)\n",
        "            u = objet.vecteur_normal(A)\n",
        "            w = vecteur(A, obs)\n",
        "            w = prod(1 / norme(w), w)\n",
        "            h = add(v, w)\n",
        "            h = prod(1 / norme(h), h)\n",
        "            I1 = prod(objet.ks, L.I)\n",
        "            I = add(I, prod(max(0, cosangle(u, h)) ** 100, I1))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def ambient(A, objet, scene):\n",
        "    I = (0, 0, 0)\n",
        "    for L in scene.lumieres:\n",
        "        r = Rayon(A, L.org)\n",
        "        black = False\n",
        "        for obj in scene.objets:\n",
        "            if obj.intersection(r) != []:\n",
        "                black = True\n",
        "                break\n",
        "        if not black:\n",
        "            I = add(I, prod(objet.ka, L.Ia))\n",
        "    return I"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def illumination(A, obs, objet, scene):\n",
        "    I = ambient(A, objet, scene)\n",
        "    I = add(I, diffuse(A, objet, scene))\n",
        "    I = add(I, phong(A, obs, objet, scene))\n",
        "    Ix, Iy, Iz = I\n",
        "    cx, cy, cz = objet.color\n",
        "    return (Ix * cx, Iy * cy, Iz * cz)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer3(scene, cameraSD)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Et voilà, les sphères ont atterri sur le plan.\n",
        "\n",
        "Quelle est l'étape suivante ? Rappelez-vous le modèle de Blin-Phong pour l'éclairement spéculaire. Une surface brillante est souvent aussi __réfléchissante__ : elle renvoie l'image des objets situés dans son voisinage. Nous allons nous attaquer au problème de la réflexion, qui sera le dernier problème évoqué dans ce notebook."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 7. Réflexion"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Comment prendre en compte les phénomènes de réflexion ? Lançons un rayon : ce rayon frappe un objet, puis se réfléchit selon les lois de Descartes. Le rayon réfléchi peut à nouveau se réfléchir sur un second objet, et ainsi de suite. La fonction `lancer4` est quasiment analogue à la fonction `lancer3`, à ceci près qu'elle prend également en paramètre un entier `prof` permettant d'éviter un nombre infini de réflexions."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancer4(scene, camera, prof, prof0):\n",
        "    t = ecran(camera.ny, camera.nx)\n",
        "    for i in range(camera.nx):\n",
        "        for j in range(camera.ny):\n",
        "            M = camera.point(i, j)\n",
        "            u = vecteur(camera.pos, M)\n",
        "            r = Rayon(camera.pos, u)\n",
        "            Ix, Iy, Iz = lancerRayon(r, 1, prof, prof0)\n",
        "            t[camera.ny - 1 - j][i] = (min(Ix, 1), min(Iy, 1), min(Iz, 1))\n",
        "    plt.axis('off')\n",
        "    plt.imshow(t)\n",
        "    plt.imsave('scene.png', t)\n",
        "    plt.show()"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Tout se passe dans la fonction `lancerRayon`. Celle-ci prend 4 paramètres :\n",
        "\n",
        "- un rayon `r`\n",
        "- un flottant `att` qui est un facteur d'atténuation : chaque réflexion diminue l'intensité du rayon lumineux.\n",
        "- le fameux paramètre `prof` qui contrôle le nombre de réflexions.\n",
        "\n",
        "On lance comme on en a l'habitude le rayon $r$ et on frappe un objet. Si `prof` est non nul, on réfléchit ce rayon et on rappelle récursivement `lancerRayon` avec une prof(ondeur) diminuée de 1 et un facteur d'atténuation modifié en conséquence.\n",
        "\n",
        "Rappelons à toutes fins utiles la loi de la réflexion : si le rayon a pour direction $\\overrightarrow u$ et frappe une surface en un point d'une surface de vecteur normal $\\overrightarrow n$, ce rayon se réfléchit selon une direction $\\overrightarrow v$, où \n",
        "\n",
        "$$\\overrightarrow v = \\overrightarrow u - 2 <\\overrightarrow u, \\overrightarrow n>\\overrightarrow n$$"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "def lancerRayon(r, att, prof, prof0):\n",
        "    s = []\n",
        "    for objet in scene.objets:\n",
        "        s += objet.intersection(r)\n",
        "    if s != []:\n",
        "        (objet, A) = plusProche(r.org, s)\n",
        "        I =  prod(att, illumination(A, r.org, objet, scene))\n",
        "        if prof > 0:\n",
        "            u = r.dir\n",
        "            u = normaliser(u)\n",
        "            n = objet.vecteur_normal(A)\n",
        "            n = normaliser(n)\n",
        "            v = sub(u, prod(2 * prodscal(u, n), n))\n",
        "            r1 = Rayon(A, v)\n",
        "            I1 = lancerRayon(r1, 0.5 * att, prof - 1, prof0)\n",
        "            I = add(I, I1)\n",
        "        return I\n",
        "    elif prof == prof0: return (0.5, 0.5, 1.)\n",
        "    else: return (0., 0., 0.)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "sphere1 = Sphere(centre=(1, -1, -1), rayon=1, color=(0.8, 0.8, 0.3), ka=0.1, kd=1, ks=1)\n",
        "sphere2 = Sphere(centre=(-3, -1, -3), rayon=1, color=(0.8, 0.3, 0.3), ka=0.1, kd=1, ks=0.)\n",
        "sphere3 = Sphere(centre=(0, 0, -5), rayon=2, color=(0.8, 0.8, 0.8), ka=0.1, kd=1, ks=1)\n",
        "plan = Plan(org=(0, -2, 0), normal=(0, 1, 0), color=(0.9, 0.6, 0.6), ka=0.1, kd=0.5, ks=0.6)\n",
        "lum1 = Lumiere(org=(60, 50, 100), I=(1, 1, 1), Ia=(0.5, 0.5, 0.5))\n",
        "scene = Scene([sphere1, sphere2, sphere3, plan], [lum1])"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Ce sera notre dernière expérience pour ce notebook, alors prenons notre caméra Haute Définition. Soyez patients (une vingtaine de secondes sur ma machine)."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "lancer4(scene, cameraHD, 3, 3)"
      ],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Et maintenant ? Le notebook suivant abordera la question de la __transparence__. Certains objets se laissent partiellement traverser par la lumière : un rayon arrivant sur l'objet se subdivise en deux rayons. \n",
        "\n",
        "- L'un des deux est un rayon __réfléchi__.\n",
        "- L'autre est un rayon __réfracté__. \n",
        "\n",
        "Les lois de Snell-Descartes seront nos amies :-)."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [],
      "outputs": [],
      "execution_count": null,
      "metadata": {}
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.4",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "kernel_info": {
      "name": "python3"
    },
    "nteract": {
      "version": "0.14.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}