{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5.5. Ray tracing: Cython with tuples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%load_ext cythonmagic\n",
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't use NumPy anymore for computations on 3D vectors, we use tuples which have less overhead for small vectors. We need to reimplement all element-wise computations with tuples, but that's fine since our vectors always have 3 elements only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "DBL = np.double\n",
    "ctypedef np.double_t DBL_C\n",
    "from libc.math cimport sqrt\n",
    "\n",
    "cdef int w, h\n",
    "w, h = 200, 200\n",
    "\n",
    "cdef dot(tuple x, tuple y):\n",
    "    return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]\n",
    "\n",
    "cdef normalize(tuple x):\n",
    "    cdef double n\n",
    "    n = sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])\n",
    "    return (x[0] / n, x[1] / n, x[2] / n)\n",
    "\n",
    "cdef max(double x, double y):\n",
    "    return x if x > y else y\n",
    "\n",
    "cdef min(double x, double y):\n",
    "    return x if x < y else y\n",
    "\n",
    "cdef clip_(double x, double m, double M):\n",
    "    return min(max(x, m), M)\n",
    "\n",
    "cdef clip(tuple x, double m, double M):\n",
    "    return (clip_(x[0], m, M), clip_(x[1], m, M), clip_(x[2], m, M),)\n",
    "\n",
    "cdef add(tuple x, tuple y):\n",
    "    return (x[0] + y[0], x[1] + y[1], x[2] + y[2])\n",
    "\n",
    "cdef subtract(tuple x, tuple y):\n",
    "    return (x[0] - y[0], x[1] - y[1], x[2] - y[2])\n",
    "\n",
    "cdef minus(tuple x):\n",
    "    return (-x[0], -x[1], -x[2])\n",
    "\n",
    "cdef multiply(tuple x, tuple y):\n",
    "    return (x[0] * y[0], x[1] * y[1], x[2] * y[2])\n",
    "    \n",
    "cdef multiply_s(tuple x, double c):\n",
    "    return (x[0] * c, x[1] * c, x[2] * c)\n",
    "    \n",
    "cdef intersect_sphere(tuple O, \n",
    "                      tuple D, \n",
    "                      tuple S, \n",
    "                      double R):\n",
    "    # Return the distance from O to the intersection of the ray (O, D) with the \n",
    "    # sphere (S, R), or +inf if there is no intersection.\n",
    "    # O and S are 3D points, D (direction) is a normalized vector, R is a scalar.\n",
    "    cdef double a, b, c, disc, distSqrt, q, t0, t1\n",
    "    cdef tuple OS\n",
    "    \n",
    "    a = dot(D, D)\n",
    "    OS = subtract(O, S)\n",
    "    b = 2 * dot(D, OS)\n",
    "    c = dot(OS, OS) - R * R\n",
    "    disc = b * b - 4 * a * c\n",
    "    if disc > 0:\n",
    "        distSqrt = sqrt(disc)\n",
    "        q = (-b - distSqrt) / 2.0 if b < 0 else (-b + distSqrt) / 2.0\n",
    "        t0 = q / a\n",
    "        t1 = c / q\n",
    "        t0, t1 = min(t0, t1), max(t0, t1)\n",
    "        if t1 >= 0:\n",
    "            return t1 if t0 < 0 else t0\n",
    "    return float('inf')\n",
    "\n",
    "cdef trace_ray(tuple O, tuple D,):\n",
    "    \n",
    "    cdef double t, radius, diffuse, specular_k, specular_c, DF, SP\n",
    "    cdef tuple M, N, L, toL, toO, col_ray, \\\n",
    "        position, color, color_light, ambient\n",
    "\n",
    "    # Sphere properties.\n",
    "    position = (0., 0., 1.)\n",
    "    radius = 1.\n",
    "    color = (0., 0., 1.)\n",
    "    diffuse = 1.\n",
    "    specular_c = 1.\n",
    "    specular_k = 50.\n",
    "    \n",
    "    # Light position and color.\n",
    "    L = (5., 5., -10.)\n",
    "    color_light = (1., 1., 1.)\n",
    "    ambient = (.05, .05, .05)\n",
    "    \n",
    "    # Find first point of intersection with the scene.\n",
    "    t = intersect_sphere(O, D, position, radius)\n",
    "    # Return None if the ray does not intersect any object.\n",
    "    if t == float('inf'):\n",
    "        return\n",
    "    # Find the point of intersection on the object.\n",
    "    M = (O[0] + D[0] * t, O[1] + D[1] * t, O[2] + D[2] * t)\n",
    "    N = normalize(subtract(M, position))\n",
    "    toL = normalize(subtract(L, M))\n",
    "    toO = normalize(subtract(O, M))\n",
    "    DF = diffuse * max(dot(N, toL), 0)\n",
    "    SP = specular_c * max(dot(N, normalize(add(toL, toO))), 0) ** specular_k\n",
    "    \n",
    "    return add(ambient, add(multiply_s(color, DF), multiply_s(color_light, SP)))\n",
    "\n",
    "def run():\n",
    "    cdef DBL_C[:,:,:] img = np.zeros((h, w, 3))\n",
    "    cdef tuple img_\n",
    "    cdef int i, j\n",
    "    cdef double x, y\n",
    "    cdef tuple O, Q, D, col_ray\n",
    "        \n",
    "    # Camera.\n",
    "    O = (0., 0., -1.)  # Position.\n",
    "        \n",
    "    # Loop through all pixels.\n",
    "    for i in range(w):\n",
    "        for j in range(h):\n",
    "            x = -1. + 2*float(i)/w\n",
    "            y = -1. + 2*float(j)/h\n",
    "            Q = (x, y, 0.)\n",
    "            D = normalize(subtract(Q, O))\n",
    "            col_ray = trace_ray(O, D)\n",
    "            if col_ray is None:\n",
    "                continue\n",
    "            img_ = clip(col_ray, 0., 1.)\n",
    "            img[h - j - 1, i, 0] = img_[0]\n",
    "            img[h - j - 1, i, 1] = img_[1]\n",
    "            img[h - j - 1, i, 2] = img_[2]\n",
    "    return img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img = run()\n",
    "plt.imshow(img);\n",
    "plt.xticks([]); plt.yticks([]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%timeit run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).\n",
    "\n",
    "> [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages)."
   ]
  }
 ],
 "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.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}