rt-original.py

#!/usr/bin/python

from math import sqrt, ceil
import random
import sys

# Define a vector class with constructor and operator: 'v'
class v(object):
    # Vector has three float attributes.
    def __init__(self, x=0.0, y=0.0, z=0.0):
        self.x = float(x)
        self.y = float(y)
        self.z = float(z)

    def __add__(self, r): # +, vector add
        return v(self.x+r.x, self.y+r.y, self.z+r.z)

    def __mul__(self, f): # *, vector scaling
        return v(self.x*f, self.y*f, self.z*f)

    def __mod__(self, r): # %, vector dot product
        return self.x * r.x + self.y * r.y + self.z * r.z

    def __xor__(self, r): # ^, vector cross product
        return v(self.y*r.z-self.z*r.y,
                      self.z*r.x-self.x*r.z,
                      self.x*r.y-self.y*r.x)

    def __neg__(self): # -, Used later for normalizing the vector. 'not' cannot be used
        return self * (1.0/sqrt(self % self))

# The set of sphere positions describing the world.
# Those integers are in fact bit vectors.
G = [247570,280596,280600,249748,18578,18577,231184,16,16]

def R():
    return random.random()

def T(o, d, t, n):
    """
       The intersection test for line [o,v].
       Return 2 if a hit was found (and also return distance t and bouncing ray n).
       Return 0 if no hit was found but ray goes upward
       Return 1 if no hit was found but ray goes downward
    """
    t = 1e9
    m = 0
    p = -o.z/d.z

    if .01 < p:
        t = p
        n = v(0, 0, 1)
        m = 1

    # The world is encoded in G, with 9 lines and 19 columns
    for k in range(18, -1, -1): # For each columns of objects
        for j in range(8, -1, -1): # For each line on that columns
            if G[j] & 1 << k: # For this line j, is there a sphere at column i ?
                # There is a sphere but does the ray hit it ?
                p = o + v(-k, 0, -j-4)
                b = p % d
                c = p % p - 1
                q = b*b-c

                # Does the ray hit the sphere ?
                if q > 0:
                    # It does, compute the distance camera-sphere
                    s = -b-sqrt(q)

                    if 0.01 < s < t:
                        # So far this is the minimum distance, save it. And also compute the bouncing ray vector into 'n'
                        t = s
                        n = -(p+d*t)
                        m = 2

    return m, t, n

def S(o, d):
    """
        (S)ample the world and return the pixel color for
         a ray passing by point o (Origin) and d (Direction)
    """
    t = 0.0
    n = v()

    # Search for an intersection ray Vs World.
    m, t, n = T(o, d, t, n)

    if not m: # m == 0
        # No sphere found and the ray goes upward: Generate a sky color
        return v(.7,.6, 1) * ((1-d.z) ** 4)

    # A sphere was maybe hit

    h = o+d*t # h = intersection coordinate
    l = -(v(9+R(), 9+R(), 16) + h * -1) # 'l' = direction to light (with random delta for soft-shadows).
    r = d + n * (n%d*-2) # r = The half-vector

    # lambertian factor
    b = l % n

    # illumination factor (lambertian coefficient > 0 or in shadow)?
    if b < 0 or T(h, l, t, n)[0]:
        b = 0

    # calculate the color 'p' with diffuse and specular component
    p = (l % r * (1 if b > 0 else 0))**99

    if m %2 == 1: # m -- 1
        h = h * .2 # No sphere was hit and the ray was going downward: Generate a floor color
        return(v(3,1,1)if(ceil(h.x)+ceil(h.y))%2 else v(3,3,3))*(b*.2+.1)

    # m == 2 sphere was hit. Cast an ray bouncing from the sphere surface.
    return v(p, p, p) + S(h, r) * .5 # Attenuate color by 50% since it is bouncing (* .5)


if __name__ == '__main__':
    # The main function. It generates a PPM image to stdout.
    # Usage of the program is hence: pypy rt.py > erk.ppm
    sys.stdout.write("P6 512 512 255 ")

    # the - is for normalizing each vector with '-' operator
    g = -v(-6, -16, 0)  # camera direction
    a = -(v(0,0,1)^g)*.002 # Camera up vector...Seem Z is pointing up :/ WTF !
    b = -(g^a) * 0.002 # The right vector, obtained via traditional cross-product
    c = (a+b)*-256+g # WTF ? See https://news.ycombinator.com/item?id=6425965 for more.

    for y in range(511, -1, -1):
        for x in range(511, -1, -1):
            # reuse the vector class to store not XYZ but a RGB pixel color
            p = v(13, 13, 13)

            # cast 64 rays per pixel (For blur (stochastic sampling) and soft-shadows.
            for r in range(63, -1, -1):
                # The delta to apply to the origin of the view (For Depth of View blur).
                t = a*(R()-.5)*99+b*(R()-.5)*99
                # Set the camera focal point v(17,16,8) and Cast the ray
                # Accumulate the color returned in the p variable
                p = S(v(17, 16, 8)+t,
                           -(t*-1+(a*(R()+x)+b*(y+R())+c)*16) # Ray Direction with random deltas
                                                              # for stochastic sampling
                          )*3.5+p # +p for color accumulation

            sys.stdout.write("%c%c%c" % (int(p.x), int(p.y), int(p.z)))