Basic ray marching done in python, based of this tutorial: http://jamie-wong.com/2016/07/15/ray-marching-signed-distance-functions/

raymarching.py

import math
from PIL import Image

class Matrix:
    """not really a matrix object"""
    @classmethod
    def fromVector(cls, v):
        m = [[] for _ in range(3)]
        m[0].append(v.x)
        m[1].append(v.y)
        m[2].append(v.z)
        return m
    @classmethod
    def toVector(cls, matrix):
        return Vector(matrix[0][0],
                      matrix[1][0],
                      matrix[2][0] if len(matrix) > 2 else 0)
    @classmethod
    def mult(cls, a, b):
        if isinstance(b, Vector):
            m = cls.fromVector(b)
            r = cls.mult(a, m)
            return cls.toVector(r)

        colsA, rowsA = len(a[0]), len(a)
        colsB, rowsB = len(b[0]), len(b)

        assert colsA == rowsB

        result = []
        for j in range(rowsA):
            result.append([])
            for i in range(colsB):
                s = 0
                for n in range(colsA):
                    s += a[j][n] * b[n][i]
                result[j].append(s)
        return result

class Vector:
    def __init__(self, x=0, y=0, z=0):
        self.x = x
        self.y = y
        self.z = z

    def __getitem__(self, index):
        return (self.x, self.y, self.z)[index]
    def __len__(self): return 3
    def __iter__(self):
        for i in range(len(self)):
            yield self[i]

    def __repr__(self): return f"<Vector({','.join(map(str, self))})>"
    
    def length(self):
        return math.sqrt(self.x**2 + self.y**2 + self.z**2)

    @classmethod
    def _convertIfNotVec(cls, other):
        if not isinstance(other, cls): return cls(other, other, other)
        else: return other

    def dot(self, other):
        other = Vector._convertIfNotVec(other)
        return sum(self * other)

    def norm(self):
        return math.sqrt(self.dot(self))

    def normalize(self):
        n = self.norm()
        p = 1 / n if n != 0 else 0
        return self * p

    def floor(self):
        return Vector(*map(int, self))

    def cross(self, other):
        return Vector(self.y*other.z - self.z*other.y,
                      self.z*other.x - self.x*other.z,
                      self.x*other.y - self.y*other.x)

    def __abs__(self): return Vector(*map(abs, self))

    def __add__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: sum(x), zip(self, other)))
    def __radd__(self, other): return self + other

    def __sub__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: x[0]-x[1], zip(self, other)))
    def __rsub__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: x[1]-x[0], zip(self, other)))

    def __mul__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: x[0]*x[1], zip(self, other)))
    def __rmul__(self, other): return self * other

    def __truediv__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: x[0]/x[1], zip(self, other)))
    def __rtruediv__(self, other):
        other = Vector._convertIfNotVec(other)
        return Vector(*map(lambda x: x[1]/x[0], zip(self, other)))
    
    def __neg__(self):
        return self * -1

    @classmethod
    def min(cls, a, b):
        return cls(*map(lambda x: min(*x), zip(a, b)))
    @classmethod
    def max(cls, a, b):
        return cls(*map(lambda x: max(*x), zip(a, b)))

def linearMap(x,xm,xl,ym,yl): return (x-xm)/(xl-xm)*(yl-ym)+ym
def lerp(a, b, k): return linearMap(k, 0, 1, a, b)
def clamp(x, a, b):
    if x < a: return a
    if x > b: return b
    return x

def smoothMin(a, b, k):
    h = clamp(0.5+0.5*(b-a)/k, 0, 1);
    return lerp(b, a, h) - k*h*(1-h);

def intersectSDF(a, b): return max(a, b)
def unionSDF(a, b):     return min(a, b)
def differenceSDF(a, b):return max(a, -b)

def sphereSDF(p, radius, pos=Vector()):
    p += pos
    return p.length() - radius

def torusSDF(p, t, pos=Vector()):
    p += pos
    q = Vector(Vector(p.x, p.z).length()-t.x, p.y)
    return q.length()-t.y

def boxSDF(p, b, pos=Vector()):
    p += pos
    d = abs(p) - b
    return Vector.max(d, Vector()).length() + min(max(d.x, d.y, d.z), 0)

minDist = 0
maxDist = 100
maxSteps = 255
# Minimum distance to be considered a hit
epsilon = 0.0001
offset = 0

def shortestDist(eye:Vector, dir:Vector, start:float, end:float):
    depth = start
    for _ in range(maxSteps):
        dist = sceneSDF(eye + dir * depth)
        if dist < epsilon: return depth
        depth += dist
        if depth >= end:
            return end
    return end

def rayDir(fov, size, fragCoord):
    xy = fragCoord - size / 2
    z = size.y / math.tan((fov * math.pi / 180) / 2)
    return Vector(xy.x, xy.y, -z).normalize()

def estimateNormal(p):
    return Vector(
        sceneSDF(Vector(p.x + epsilon, p.y, p.z)) - sceneSDF(Vector(p.x - epsilon, p.y, p.z)),
        sceneSDF(Vector(p.x, p.y + epsilon, p.z)) - sceneSDF(Vector(p.x, p.y - epsilon, p.z)),
        sceneSDF(Vector(p.x, p.y, p.z + epsilon)) - sceneSDF(Vector(p.x, p.y, p.z - epsilon)),
        ).normalize()

def rotateX(theta):
    c = math.cos(theta)
    s = math.sin(theta)
    return [
        [1, 0, 0],
        [0, c,-s],
        [0, s, c],
    ]
def rotateZ(theta):
    c = math.cos(theta)
    s = math.sin(theta)
    return [
        [ c, 0, s],
        [ 0, 1, 0],
        [-s, 0, c],
    ]
def lookAt(eye, look, dir):
    delta = eye - look
    pitch = -math.atan2(delta.y, math.sqrt(delta.x * delta.x + delta.z * delta.z))
    yaw = math.pi / 2 - math.atan2(delta.z,delta.x)

    rotPitch = rotateX(pitch)
    rotYaw = rotateZ(yaw)

    dir = Matrix.mult(rotPitch, dir)
    dir = Matrix.mult(rotYaw, dir)
    return dir.normalize()    

def reflect(I, N): return I - 2 * N.dot(I) * N

def phongContribForLight(k_d, k_s, alpha, p, eye, lightPos, lightIntensity):
    N = estimateNormal(p)
    L = (lightPos - p).normalize()
    V = (eye - p).normalize()
    R = reflect(-L, N).normalize()
    
    dotLN = L.dot(N)
    dotRV = R.dot(V)
    
    if dotLN < 0:
        return Vector(0, 0, 0)
    
    if dotRV < 0:
        return lightIntensity * (k_d * dotLN)

    return lightIntensity * (k_d * dotLN + k_s * pow(dotRV, alpha))

def phongIllumination(k_a, k_d, k_s, alpha, p, eye):
    lights = [(Vector(4 * math.sin(offset), -2, 4 * math.cos(offset)), Vector(1, 1, 1)),
              (Vector(2 * math.sin(0.37 * offset), 2 * math.cos(0.37 * offset), 2), Vector(1, 1, 1))]
    ambientLight = Vector(1, 1, 1) * k_a
    color = ambientLight * k_a
    for light in lights:
        color += phongContribForLight(k_d, k_s, alpha, p, eye, light[0], light[1])
    return color

def sceneSDF(p):
    #torus = torusSDF(p, Vector(1, 0.5))
    sphere = sphereSDF(p, 1.2, pos=Vector(math.sin(offset)*5, 0, 0))
    cube = boxSDF(p, Vector()+1, pos=Vector(0, -2.5, 0))
    return smoothMin(cube, sphere, 1)

offset = 0
frames = 30
def render(width, height):
    global offset
    img = Image.new('RGB', (width, height), (100, 100, 100)) 
    data = list(img.getdata())
    for y in range(height):
        for x in range(width):
            i = y * width + x
            dir = rayDir(45, Vector(width, height), Vector(x, y))
            eye = Vector(8, -5, 7)

            dir = lookAt(eye, Vector(0, 0, 0), dir)

            dist = shortestDist(eye, dir, minDist, maxDist)

            if dist > maxDist - epsilon:
                data[i] = (0,0,0)
            else:
                p = eye + dir * dist
                K_a = Vector() + 0.2
                K_d = Vector(0.5) + 0.2
                K_s = Vector() + 1
                shininess = 10
                color = phongIllumination(K_a, K_d, K_s, shininess, p, eye)
                data[i] = tuple(Vector.min((color * 255).floor(),Vector(255, 255, 255)))
    img.putdata(data)
    offset += 2*math.pi/frames
    return img
img = render(160, 90).resize((1280, 720))
img.save('output.png')
img.show()
# def nf(n, digits): return ('0'*digits)[len(str(n)):] + str(n)
# for i in range(frames):
#     render(120, 60).resize((1280, 720)).save(f'{nf(i,len(str(frames)))}.png')