Is there a good way to speed up raytracing in Python?

Question

I've been following the Raytracing in one weekend tutorial series (which is coded in C++), but I've run into a problem. I know Python is not really the language for raytracing, and it shows. Rendering 4 spheres with 100 samples per pixel at 600x300px with multithreading took a little over an quarter of an hour.

I figured that multithreading wasn't really doing anything. I needed multiprocessing. I tried to implement it, but it gave me the error that pygame.Surface wasn't serializable. Is there a solution for this, or maybe even a better way of implementing this? Or is the best way forward to recreate this in C#/C++?

The render:

Multithreaded code:

vec3.py

import math
import random

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

    def __add__(self, other):
        return vec3(
            self.x + other.x,
            self.y + other.y,
            self.z + other.z
        )

    def __sub__(self, other):
        return vec3(
            self.x - other.x,
            self.y - other.y,
            self.z - other.z
        )

    def __mul__(self, val):
        if type(val) == vec3:
            return vec3(
                self.x * val.x,
                self.y * val.y,
                self.z * val.z
            )
        else:
            return vec3(
                self.x * val,
                self.y * val,
                self.z * val
            )

    def __truediv__(self, val):
        return vec3(
            self.x / val,
            self.y / val,
            self.z / val
        )
    
    def __pow__(self, pow):
        return vec3(
            self.x ** pow,
            self.y ** pow,
            self.z ** pow
        )
    
    def __rshift__(self, val):
        return vec3(
            self.x >> val,
            self.y >> val,
            self.z >> val
        )
    
    def __repr__(self):
        return "vec3({}, {}, {})".format(self.x, self.y, self.z)

    def __str__(self):
        return "({}, {}, {})".format(self.x, self.y, self.z)

    def dot(self, other):
        return (
            self.x * other.x +
            self.y * other.y +
            self.z * other.z
        )
    
    def cross(self, other):
        return vec3(
            (self.y * other.z) - (self.z * other.y),
           -(self.x * other.z) - (self.z * other.x),
            (self.x * other.y) - (self.y * other.x)
        )

    def length(self):
        return math.sqrt(
            (self.x ** 2) + 
            (self.y ** 2) + 
            (self.z ** 2)
        )

    def squared_length(self):
        return (
            (self.x ** 2) + 
            (self.y ** 2) + 
            (self.z ** 2))

    def make_unit_vector(self):
        return self / (self.length() + 0.0001)
    
    def reflect(self, normal):
        return self - normal * (self.dot(normal) * 2)

Raytracing.py

import uuid
import math
import pygame
from random import random

from vec3 import vec3
from threading import Thread
#import time

imprecision = 0.00000001
mindist = 0.001
maxdist = 999999999

class Material:
    def random_in_unit_sphere():
        p = vec3(random(), random(), random())*2 - vec3(1, 1, 1)
        while p.squared_length() >= 1:
            p = vec3(random(), random(), random())*2 - vec3(1, 1, 1)
        return p
    class Matte:
        def __init__(self, albedo):
            self.albedo = albedo
        
        def scatter(self, ray, rec, scattered):
            target = rec.p + rec.normal + Material.random_in_unit_sphere()
            scattered.origin = rec.p
            scattered.direction = target - rec.p
            return True, self.albedo
    
    class Metal:
        def __init__(self, albedo, fuzz = 0):
            self.albedo = albedo
            self.fuzz = fuzz

        def scatter(self, ray, rec, scattered):
            reflected = ray.direction.make_unit_vector().reflect(rec.normal)
            scattered.origin = rec.p
            scattered.direction = reflected + (Material.random_in_unit_sphere() * self.fuzz)
            return (scattered.direction.dot(rec.normal) > 0), self.albedo

    class Transparent:
        def __init__(self, refractive_index):
            self.ri = refractive_index
        
        def refract(self, v, normal, ni_over_nt):
            uv = v.make_unit_vector()
            dt = uv.dot(normal)
            discriminant = 1 - (ni_over_nt * ni_over_nt * (1-dt*dt))
            if discriminant > 0:
                refracted = (uv - normal*dt) * ni_over_nt - (normal * math.sqrt(discriminant))
                return True, refracted
            else:
                return False, None

        def scatter(self, ray, rec, scattered):
            reflected = ray.direction.reflect(rec.normal)

            if ray.direction.dot(rec.normal) > 0:
                outward_normal = 0 - rec.normal
                ni_over_nt = self.ri
            else:
                outward_normal = rec.normal
                ni_over_nt = 1 / self.ri
            
            refracted = self.refract(ray.direction, outward_normal, ni_over_nt)
            return bool(refracted[0]), vec3(1, 1, 1)

class Object:
    class Sphere:
        def __init__(self, pos, r, material=Material.Matte(0.5)):
            self.center = pos
            self.r = r
            self.material = material

        def hit(self, ray, t_min, t_max, rec):
            oc = ray.origin - self.center
            a = ray.direction.dot(ray.direction)
            b = oc.dot(ray.direction)
            c = oc.dot(oc) - (self.r * self.r)
            discriminant = (b*b) - (a*c)
            if discriminant < 0:
                return False
            else:
                t = (-b - math.sqrt(discriminant)) / (a + imprecision)
                if (t_min < t < t_max):
                    rec.t = t
                    rec.p = ray.point_at_parameter(t)
                    rec.normal = (rec.p - self.center) / self.r
                    rec.material = self.material
                    return True

                t = (-b + math.sqrt(discriminant)) / (a + imprecision)
                if (t_min < t < t_max):
                    rec.t = t
                    rec.p = ray.point_at_parameter(t)
                    rec.normal = (rec.p - self.center) / self.r
                    rec.material = self.material
                    return True
class Scene:
    def genUUID(self):
        id = uuid.uuid4()
        while id in self.objects.keys():
            id = uuid.uuid4()
        return id

    def __init__(self):
        self.objects = {}
    
    def add(self, object):
        id = self.genUUID()
        self.objects[id] = object
        return id

    def hit(self, ray, t_min, t_max, hit_record):
        tmp_hit_record = Camera.HitRec()
        hit_any = False
        closest = maxdist
        for object in self.objects.keys():
            if self.objects[object].hit(ray, t_min, closest, tmp_hit_record):
                hit_any = True
                closest = tmp_hit_record.t
                hit_record.t = tmp_hit_record.t
                hit_record.p = tmp_hit_record.p
                hit_record.normal = tmp_hit_record.normal
                hit_record.material = tmp_hit_record.material
        return hit_any

class Camera:
    class Ray:
        def __init__(self, A, B):
            self.origin = A
            self.direction = B

        def point_at_parameter(self, t):
            return self.origin + (self.direction * t)
    class HitRec:
        def __init__(self):
            self.t = 0
            self.p = vec3(0, 0, 0)
            self.normal = vec3(0, 0, 0)

    def __init__(self, w, screen, scene = Scene(), screen_ratio = 16/9, focal_length = 1):
        self.w = w
        self.h = int(w / screen_ratio)
        self.screen = screen
        self.scene = scene
        self.screen_ratio = screen_ratio
        self.screen_unit_width = 4
        self.focal_length = focal_length
        self.sampels_per_pixel = 1
        self.max_hit_depth = 5
        self.max_block_size = 100

#        self.depthT = pygame.Surface((w, self.h))

        self.origin            = vec3(0,  0,  0)
        self.horizontal        = vec3(self.screen_unit_width, 0,  0)
        self.vertical          = vec3(0, (-self.screen_unit_width / self.screen_ratio),  0)
        self.top_left_corner   = \
            self.origin - (self.horizontal/2) - (self.vertical/2) - vec3(0, 0, self.focal_length)
    
    def color(self, ray, depth=0):

        # Check for sphere hit
        rec = self.HitRec()
        if self.scene.hit(ray, mindist, maxdist, rec):
            scattered = self.Ray(vec3(0, 0, 0), vec3(0, 0, 0))
            hit = rec.material.scatter(ray, rec, scattered) # Scattered is the output ray
            if depth < self.max_hit_depth and hit[0]:
                return self.color(scattered, depth + 1) * hit[1]
            else:
                return hit[1]
#            target = rec.p + rec.normal + self.random_in_unit_sphere()
#            return self.color(self.Ray(rec.p, target - rec.p)) * 0.5
##            return vec3(rec.normal.x + 1, rec.normal.y + 1, rec.normal.z + 1) * 0.5, rec.t

        # Background if nothing's hit
        unit_dir = ray.direction.make_unit_vector()
        t = 0.5 * (unit_dir.y + 1)
        return (vec3(1, 1, 1) * (1 - t)) + (vec3(0.5, 0.7, 1) * t)#, maxdist

    def getRay(self, u, v):
        return self.Ray(
                self.origin,
                self.top_left_corner + (self.horizontal * (u + imprecision)) + (self.vertical * (v + imprecision))
            )

    def renderSquare(self, xoff, yoff, w, h):
        pygame.draw.rect(self.screen, (255, 255, 255), (xoff, yoff, w, h))
        pygame.display.update()
        for y in range(h):
            for x in range(w):
                col = vec3(0, 0, 0)
                for n in range(self.sampels_per_pixel):
                    u = (x + xoff + random()-0.5) / self.w
                    v = (y + yoff + random()-0.5) / self.h
                    r = self.getRay(u, v)
                    col += (self.color(r) / self.sampels_per_pixel)
                col = vec3(math.sqrt(col.x), math.sqrt(col.y), math.sqrt(col.z))
                self.screen.set_at((x + xoff, y + yoff), (col.x * 255, col.y * 255, col.z * 255))
            pygame.display.update()

    def render(self):
        blocks = []
        for i in range(math.ceil(self.w / self.max_block_size)):
            for j in range(math.ceil(self.h / self.max_block_size)):
                blocks.append(Thread( target = self.renderSquare, args=(
                        i * self.max_block_size,
                        j * self.max_block_size,
                        self.max_block_size,
                        self.max_block_size
                    )))
        for job in blocks:
            job.start()
        for job in blocks:
            job.join()

if __name__ == '__main__':
    from Raytracing import Camera, Scene, Object
    import pygame

    pygame.init()

    WIDTH, HEIGHT = 600, 300

    screen = pygame.display.set_mode((WIDTH, HEIGHT))
    pygame.display.set_caption('Multithreaded raytracer')

    scene = Scene()

    scene.add(Object.Sphere(vec3( 0, 0, -1), 0.5, Material.Matte(vec3(0, 0.3, 0.8))))
    scene.add(Object.Sphere(vec3(-1, 0, -1), 0.5, Material.Metal(vec3(0.8, 0.8, 0.8), 0.1)))
    scene.add(Object.Sphere(vec3( 1, 0, -1), 0.5, Material.Metal(vec3(0.8, 0.6, 0.2), 1)))

    scene.add(Object.Sphere(vec3(0, -100.5, -1), 100, Material.Matte(vec3(0.8, 0.8, 0))))

    tracer = Camera(WIDTH, screen, scene, screen_ratio=2 / 1)
#    tracer.render()
    tracer.sampels_per_pixel = 100
    tracer.render()

    loop = True
    while loop:
        for event in pygame.event.get():
            if event.type==pygame.KEYDOWN:
                if event.key == pygame.K_ESCAPE:
                    loop = False
            if event.type==pygame.QUIT:
                loop = False
                

        pygame.display.update()
    
    pygame.quit()