notebooks/chapter05_hpc/05_ray_1.ipynb

{
 "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: pure Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we will render a sphere with a diffuse and specular material. The principle is to model a scene with a light source and a camera, and use the physical properties of light propagation to calculate the light intensity and color of every pixel of the screen."
   ]
  },
  {
   "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": [
    "w, h = 200, 200  # Size of the screen in pixels.\n",
    "\n",
    "def normalize(x):\n",
    "    # This function normalizes a vector.\n",
    "    x /= np.linalg.norm(x)\n",
    "    return x\n",
    "\n",
    "def intersect_sphere(O, D, S, R):\n",
    "    # Return the distance from O to the intersection \n",
    "    # of the ray (O, D) with the sphere (S, R), or \n",
    "    # +inf if there is no intersection.\n",
    "    # O and S are 3D points, D (direction) is a \n",
    "    # normalized vector, R is a scalar.\n",
    "    a = np.dot(D, D)\n",
    "    OS = O - S\n",
    "    b = 2 * np.dot(D, OS)\n",
    "    c = np.dot(OS, OS) - R*R\n",
    "    disc = b*b - 4*a*c\n",
    "    if disc > 0:\n",
    "        distSqrt = np.sqrt(disc)\n",
    "        q = (-b - distSqrt) / 2.0 if b < 0 \\\n",
    "            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 np.inf\n",
    "\n",
    "def trace_ray(O, D):\n",
    "    # Find first point of intersection with the scene.\n",
    "    t = intersect_sphere(O, D, position, radius)\n",
    "    # No intersection?\n",
    "    if t == np.inf:\n",
    "        return\n",
    "    # Find the point of intersection on the object.\n",
    "    M = O + D * t\n",
    "    N = normalize(M - position)\n",
    "    toL = normalize(L - M)\n",
    "    toO = normalize(O - M)\n",
    "    # Ambient light.\n",
    "    col = ambient\n",
    "    # Lambert shading (diffuse).\n",
    "    col += diffuse * max(np.dot(N, toL), 0) * color\n",
    "    # Blinn-Phong shading (specular).\n",
    "    col += specular_c * color_light * \\\n",
    "        max(np.dot(N, normalize(toL + toO)), 0) \\\n",
    "           ** specular_k\n",
    "    return col\n",
    "\n",
    "def run():\n",
    "    img = np.zeros((h, w, 3))\n",
    "    # Loop through all pixels.\n",
    "    for i, x in enumerate(np.linspace(-1., 1., w)):\n",
    "        for j, y in enumerate(np.linspace(-1., 1., h)):\n",
    "            # Position of the pixel.\n",
    "            Q[0], Q[1] = x, y\n",
    "            # Direction of the ray going through the optical center.\n",
    "            D = normalize(Q - O)\n",
    "            # Launch the ray and get the color of the pixel.\n",
    "            col = trace_ray(O, D)\n",
    "            if col is None:\n",
    "                continue\n",
    "            img[h - j - 1, i, :] = np.clip(col, 0, 1)\n",
    "    return img\n",
    "\n",
    "# Sphere properties.\n",
    "position = np.array([0., 0., 1.])\n",
    "radius = 1.\n",
    "color = np.array([0., 0., 1.])\n",
    "diffuse = 1.\n",
    "specular_c = 1.\n",
    "specular_k = 50\n",
    "\n",
    "# Light position and color.\n",
    "L = np.array([5., 5., -10.])\n",
    "color_light = np.ones(3)\n",
    "ambient = .05\n",
    "\n",
    "# Camera.\n",
    "O = np.array([0., 0., -1.])  # Position.\n",
    "Q = np.array([0., 0., 0.])  # Pointing to.\n",
    "\n",
    "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
}