A (slow, badly written) Python implementation of https://tmcw.github.io/literate-raytracer/.

trace.py

from PIL import Image
import tqdm
import math
import vector


height = 320
width = 240

scene = {}

scene['camera'] = {
    'point': [0, 1.8, 10],
    'fieldOfView': 45,
    'vector': [0, 3, 0]
}

scene['lights'] = [[-30, -10, 20]]

scene['objects'] = [
    {
        'type': 'sphere',
        'point': [0, 3.5, 0],
        'color': [155, 200, 155],
        'specular': 0.2,
        'lambert': 0.7,
        'ambient': 0.1,
        'radius': 3
    },
    {
        'type': 'sphere',
        'point': [1, 2, 1],
        'color': [155, 155, 155],
        'specular': 0.1,
        'lambert': 0.9,
        'ambient': 0.0,
        'radius': 0.4
    },
    {
        'type': 'sphere',
        'point': [2, 1, 1],
        'color': [255, 255, 255],
        'specular': 0.2,
        'lambert': 0.7,
        'ambient': 0.1,
        'radius': 1
    },



]

img = Image.new('RGB', (width, height))
data = img.load()

def trace(ray, scene, depth):
    if depth > 3:
        return
    distObject = intersectScene(ray, scene)

    if distObject[0] == math.inf:
        return vector.WHITE

    dist = distObject[0]
    obj = distObject[1]

    pointAtTime = vector.add(ray['point'], vector.scale(ray['vector'], dist))

    return surface(ray, scene, obj, pointAtTime, sphereNormal(obj, pointAtTime), depth)

def intersectScene(ray, scene):
    closest = [math.inf, None]

    for obj in scene['objects']:
        dist= sphereIntersection(obj, ray)
        if dist != None and dist < closest[0]:
            closest = [dist, obj]

    return closest

def sphereIntersection(sphere, ray):
    eye_to_center = vector.subtract(sphere['point'], ray['point'])
    v = vector.dotProduct(eye_to_center, ray['vector'])
    eoDot = vector.dotProduct(eye_to_center, eye_to_center)
    discriminant = (sphere['radius'] * sphere['radius']) - eoDot + (v*v)

    if discriminant < 0:
        return
    else:
        return v - math.sqrt(discriminant)

def sphereNormal(sphere, pos):
    return vector.unitVector(vector.subtract(pos, sphere['point']))

def surface(ray, scene, obj, pointAtTime, normal, depth):
    b = obj['color']
    c = vector.ZERO
    lambertAmount = 0

    if 'lambert' in obj:
        for lightPoint in scene['lights']:
            if not isLightVisible(pointAtTime, scene, lightPoint):
                continue

            contribution = vector.dotProduct(vector.unitVector(vector.subtract(lightPoint, pointAtTime)), normal)

            if contribution > 0:
                lambertAmount += contribution

    if 'specular' in obj:
        reflectedRay = {
            'point': pointAtTime,
            'vector': vector.reflectThrough(ray['vector'], normal)
        }

        depth += 1
        reflectedColor = trace(reflectedRay, scene, depth)
        if reflectedColor is not None:
            c = vector.add(c, vector.scale(reflectedColor, obj['specular']))

    lambertAmount = min(1, lambertAmount)

    return vector.add3(c,
                       vector.scale(b, lambertAmount * obj['lambert']),
                       vector.scale(b, obj['ambient']))

def isLightVisible(pt, scene, light):
    distObject = intersectScene({
        'point': pt,
        'vector': vector.unitVector(vector.subtract(pt, light))
    }, scene)

    return distObject[0] > -0.005

def render(scene):
    camera = scene['camera']
    lights = scene['lights']
    objects = scene['objects']

    eyeVector = vector.unitVector(vector.subtract(camera['vector'], camera['point']))
    vpRight = vector.unitVector(vector.crossProduct(eyeVector, vector.UP))
    vpUp = vector.unitVector(vector.crossProduct(vpRight, eyeVector))

    fovRadians = math.pi * (camera['fieldOfView']/2)/180
    heightWidthRatio = height/width
    halfWidth = math.tan(fovRadians)
    halfHeight = heightWidthRatio * halfWidth
    cameraWidth = halfWidth * 2
    cameraHeight = halfHeight * 2
    pixelWidth = cameraWidth / (width-1)
    pixelHeight = cameraHeight / (height - 1)

    ray = {
        'point': camera['point']
    }

    for x in tqdm.tqdm(range(width)):
        for y in range(height):
            xcomp = vector.scale(vpRight, (x * pixelWidth) - halfWidth)
            ycomp = vector.scale(vpUp, (y * pixelHeight) - halfHeight)

            ray['vector'] = vector.unitVector(vector.add3(eyeVector, xcomp, ycomp))

            color = [int(x) for x in trace(ray, scene, 0)]

            data[x, y] = tuple(color)


render(scene)
img.show()

vector.py

'''
Given the use of numpy, probably don't actually need this file at all, but I thought I'd implement it anyway to keep
function names in line with original JS code.
'''


import numpy as np

'''
# Vector Operations

These are general-purpose functions that deal with vectors - in this case,
three-dimensional vectors represented as objects in the form

    { x, y, z }

Since we're not using traditional object oriented techniques, these
functions take and return that sort of logic-less object, so you'll see
`add(a, b)` rather than `a.add(b)`.
'''

# Constants
UP = [0, 1, 0]
ZERO = [0, 0, 0]
WHITE = [255, 255, 255]

# Operations


def dotProduct(a, b):
    '''
    ## [Dot Product](https://en.wikipedia.org/wiki/Dot_product)
    is different than the rest of these since it takes two vectors but
    returns a single number value.
    '''

    return np.dot(a, b)


def crossProduct(a, b):
    '''
    ## [Cross Product](https://en.wikipedia.org/wiki/Cross_product)

    generates a new vector that's perpendicular to both of the vectors
    given.
    '''

    return np.cross(a, b)


def scale(a, t):
    '''
    Enlongate or shrink a vector by a factor of `t`
    '''

    return np.multiply(a, t)


def unitVector(a):
    '''
    # [Unit Vector](http://en.wikipedia.org/wiki/Unit_vector)

    Turn any vector into a vector that has a magnitude of 1.

    If you consider that a [unit sphere](http://en.wikipedia.org/wiki/Unit_sphere)
    is a sphere with a radius of 1, a unit vector is like a vector from the
    center point (0, 0, 0) to any point on its surface.
    '''

    return scale(a, 1/length(a))


def add(a, b):
    '''
    Add two vectors to each other, by simply combining each
    of their components
    '''

    return np.add(a, b)


def add3(a, b, c):
    '''
    A version of `add` that adds three vectors at the same time. While
    it's possible to write a clever version of `Vector.add` that takes
    any number of arguments, it's not fast, so we're keeping it simple and
    just making two versions.
    '''
    return np.add(np.add(a, b), c)


def subtract(a, b):
    '''
    Subtract one vector from another, by subtracting each component
    '''

    return np.subtract(a, b)


def length(a):
    return np.linalg.norm(a)


def reflectThrough(a, normal):
    '''
    Given a vector `a`, which is a point in space, and a `normal`, which is
    the angle the point hits a surface, returna  new vector that is reflect
    off of that surface
    '''
    d = scale(normal, dotProduct(a, normal))
    return subtract(scale(d, 2), a)