#!/usr/bin/env python
# coding: utf-8

# > 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.
# 

# # 5.5. Ray tracing: naive Cython

# 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.

# In[ ]:


import numpy as np
import matplotlib.pyplot as plt


# In[ ]:


get_ipython().run_line_magic('matplotlib', 'inline')


# In[ ]:


#%load_ext cythonmagic
get_ipython().run_line_magic('load_ext', 'Cython')


# In[ ]:


get_ipython().run_cell_magic('cython', '', 'import numpy as np\ncimport numpy as np\n\nw, h = 200, 200  # Size of the screen in pixels.\n\ndef normalize(x):\n    # This function normalizes a vector.\n    x /= np.linalg.norm(x)\n    return x\n\ndef 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\ndef 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\ndef 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            depth = 0\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.\nposition = np.array([0., 0., 1.])\nradius = 1.\ncolor = np.array([0., 0., 1.])\ndiffuse = 1.\nspecular_c = 1.\nspecular_k = 50\n\n# Light position and color.\nL = np.array([5., 5., -10.])\ncolor_light = np.ones(3)\nambient = .05\n\n# Camera.\nO = np.array([0., 0., -1.])  # Position.\nQ = np.array([0., 0., 0.])  # Pointing to.\n')


# In[ ]:


img = run()
plt.imshow(img);
plt.xticks([]); plt.yticks([]);


# In[ ]:


get_ipython().run_line_magic('timeit', '-n1 -r1 run()')


# > You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).
# 
# > [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages).