#!/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: Cython with tuples

# In[ ]:


import numpy as np
import matplotlib.pyplot as plt


# In[ ]:


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


# In[ ]:


import cython


# In[ ]:


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


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

# In[ ]:


get_ipython().run_cell_magic('cython', '', "import numpy as np\ncimport numpy as np\nDBL = np.double\nctypedef np.double_t DBL_C\nfrom libc.math cimport sqrt\n\ncdef int w, h\nw, h = 200, 200\n\ncdef dot(tuple x, tuple y):\n    return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]\n\ncdef 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\ncdef max(double x, double y):\n    return x if x > y else y\n\ncdef min(double x, double y):\n    return x if x < y else y\n\ncdef clip_(double x, double m, double M):\n    return min(max(x, m), M)\n\ncdef 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\ncdef add(tuple x, tuple y):\n    return (x[0] + y[0], x[1] + y[1], x[2] + y[2])\n\ncdef subtract(tuple x, tuple y):\n    return (x[0] - y[0], x[1] - y[1], x[2] - y[2])\n\ncdef minus(tuple x):\n    return (-x[0], -x[1], -x[2])\n\ncdef multiply(tuple x, tuple y):\n    return (x[0] * y[0], x[1] * y[1], x[2] * y[2])\n    \ncdef multiply_s(tuple x, double c):\n    return (x[0] * c, x[1] * c, x[2] * c)\n    \ncdef 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\ncdef 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\ndef 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\n")


# In[ ]:


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


# In[ ]:


get_ipython().run_line_magic('timeit', '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).