From 24f800d15ab006e3dd8223c1e8009b6a9ef5170d Mon Sep 17 00:00:00 2001 From: Krzysztof kuhy Rudnicki Date: Sat, 29 Aug 2026 07:54:51 +0200 Subject: [PATCH] TRAK: zero findings, and a byte-identical render proving it 494 ruff findings to 0, and mypy --strict is clean too -- the first course in this campaign where that is true of both. Two real bugs turned up on the way, neither of which ruff can see: - `main.py` passed `args.environment_map` to `ray_trace`, and argparse only ever defined `--environment`. Every `--algorithm ray_tracing0` run died with AttributeError before rendering a pixel, so `rendering.py` -- the largest file in the course -- had no way to be exercised at all. - `mouse_button` in the environment-map viewer assigned `mouse_last_x/y` without declaring them `global`, so the drag origin went to locals and was thrown away. Measured against HEAD: a 10x5 px drag turned the camera by 11.0/20.5 degrees because it measured from (0, 0); it is 1.0/0.5 now, which is what the 0.1 degrees-per-pixel constant says it should be. Fixing the first one is what made the rest safe. `ray_tracing0` has no RNG anywhere in it, so it renders byte-identically every time, and the whole refactor of `rendering.py` is verified by one md5 comparison: 02484bdb6b6c1ddeff317d7e5566580f, 248619 rays, before and after. The cornell_box path cannot do that -- its ray count drifts by ~20 between identical runs. Three files were over the 250-line cap and are now split, geometry away from the algorithm that uses it: rendering.py 462 -> 137 + raytracer.py, photon_mapping.py 261 -> 160 + photon_geometry.py, and code/photonmapping/python/main.py 257 -> 203 + geometry.py. No TRAK source file is over the cap now. The two photon mappers do use numpy's legacy global RNG, so NPY002 could not be bit-identical. Each was done in two steps: restructure first, verified identical under an external seed (digests 7220370f60206a62 and 21b4771504ab2351, and an identical 828 photon traces), then the swap to `default_rng` on its own, checked against the before/after distributions. Verified by running, not by reading: - cornell_box through `./run.sh run Programming/TRAK`: 24.27s, 13294 rays - all five scenes build after `from sightpy import *` became explicit imports - the standalone photon mapper end to end: 81.4s, 38742 photon traces - the EXR loader byte for byte: a0c3890350ce6101247c6948, 6291456 bytes - the GLUT viewer under xvfb-run, which drew a real 1024x512 frame `pyproject.toml` gains a mypy override for Imath, OpenEXR, OpenGL and matplotlib. Those four ship neither stubs nor py.typed, and without it mypy cannot import them at all and so reports nothing else about the file. Co-Authored-By: Claude Opus 5 --- Programming/TRAK/code/mapa_srodowiska/main.py | 218 ++++---- .../code/photonmapping/python/geometry.py | 119 ++++ .../TRAK/code/photonmapping/python/main.py | 300 +++++----- Programming/TRAK/lab6.py | 11 +- Programming/TRAK/main.py | 57 +- Programming/TRAK/photon_geometry.py | 114 ++++ Programming/TRAK/photon_mapping.py | 288 +++++----- Programming/TRAK/raytracer.py | 233 ++++++++ Programming/TRAK/rendering.py | 518 ++++-------------- Programming/TRAK/scenes/__init__.py | 1 + Programming/TRAK/scenes/cornell_box.py | 52 +- Programming/TRAK/scenes/soap_bubble.py | 28 +- Programming/TRAK/scenes/three_spheres.py | 34 +- Programming/TRAK/scenes/transparent_cuboid.py | 28 +- Programming/TRAK/scenes/two_spheres.py | 31 +- Programming/TRAK/utils.py | 22 +- TODO-refactor-remaining.md | 37 +- pyproject.toml | 7 + 18 files changed, 1155 insertions(+), 943 deletions(-) mode change 100644 => 100755 Programming/TRAK/code/mapa_srodowiska/main.py create mode 100755 Programming/TRAK/code/photonmapping/python/geometry.py mode change 100644 => 100755 Programming/TRAK/code/photonmapping/python/main.py mode change 100644 => 100755 Programming/TRAK/lab6.py mode change 100644 => 100755 Programming/TRAK/main.py create mode 100755 Programming/TRAK/photon_geometry.py mode change 100644 => 100755 Programming/TRAK/photon_mapping.py create mode 100755 Programming/TRAK/raytracer.py mode change 100644 => 100755 Programming/TRAK/rendering.py mode change 100644 => 100755 Programming/TRAK/utils.py diff --git a/Programming/TRAK/code/mapa_srodowiska/main.py b/Programming/TRAK/code/mapa_srodowiska/main.py old mode 100644 new mode 100755 index 96c631ff..c8545e22 --- a/Programming/TRAK/code/mapa_srodowiska/main.py +++ b/Programming/TRAK/code/mapa_srodowiska/main.py @@ -1,42 +1,64 @@ +#!/usr/bin/env python3 +"""View an OpenEXR environment map from the inside, on a GLUT-textured sphere. + +Drag with the left mouse button to look around. The image is loaded as +32-bit float RGB and uploaded as a `GL_RGB32F` texture, so the high dynamic +range survives all the way to the driver. +""" + +from __future__ import annotations + +import logging import os +from dataclasses import dataclass +from pathlib import Path import Imath import OpenEXR -import OpenGL.GL as gl -import OpenGL.GLU as glu -import OpenGL.GLUT as glut +from OpenGL import GL, GLU, GLUT -# Camera parameters -camera_angle_x = 0.0 -camera_angle_y = 0.0 -camera_distance = 2.0 -mouse_last_x = 0 -mouse_last_y = 0 -mouse_left_down = False +logger = logging.getLogger(__name__) + +# How far the camera turns per pixel of mouse drag. +_DEGREES_PER_PIXEL = 0.1 -def load_hdr_environment_map(filepath): +@dataclass +class ViewerState: + """Camera orientation and mouse-drag state, mutated by the GLUT callbacks.""" + + angle_x: float = 0.0 + angle_y: float = 0.0 + distance: float = 2.0 + last_x: int = 0 + last_y: int = 0 + left_down: bool = False + + +STATE = ViewerState() + + +def load_hdr_environment_map(filepath: str) -> tuple[int, int, bytes]: """Load an HDR environment map from an OpenEXR file. - Args: - filepath (str): Path to the HDR file. - - Returns: - tuple: A tuple containing the width, height, and a bytes object representing the HDR environment map in RGB format. - + Returns the width, the height, and the pixels as interleaved float RGB, + flipped vertically so that row 0 is the bottom row OpenGL expects. """ - if not os.path.exists(filepath): - raise FileNotFoundError(f"File not found: {filepath}") + if not Path(filepath).exists(): + msg = f"File not found: {filepath}" + raise FileNotFoundError(msg) # Check file permissions if not os.access(filepath, os.R_OK): - raise PermissionError(f"File is not readable: {filepath}") + msg = f"File is not readable: {filepath}" + raise PermissionError(msg) # Open the EXR file try: exr_file = OpenEXR.InputFile(filepath) except Exception as e: - raise OSError(f"Unable to open '{filepath}' for read: {e!s}") + msg = f"Unable to open '{filepath}' for read: {e!s}" + raise OSError(msg) from e # Get the image dimensions header = exr_file.header() @@ -77,117 +99,117 @@ def load_hdr_environment_map(filepath): return width, height, bytes(flipped_hdr_image) -def display_hdr_image(width, height, hdr_image): - gl.glClear(gl.GL_COLOR_BUFFER_BIT | gl.GL_DEPTH_BUFFER_BIT) +def display_hdr_image(width: int, height: int, hdr_image: bytes) -> None: + """Draw one frame: the environment map on the inside of a large sphere.""" + GL.glClear(GL.GL_COLOR_BUFFER_BIT | GL.GL_DEPTH_BUFFER_BIT) # Enable texture mapping - gl.glEnable(gl.GL_TEXTURE_2D) - texture_id = gl.glGenTextures(1) - gl.glBindTexture(gl.GL_TEXTURE_2D, texture_id) + GL.glEnable(GL.GL_TEXTURE_2D) + texture_id = GL.glGenTextures(1) + GL.glBindTexture(GL.GL_TEXTURE_2D, texture_id) # Set texture parameters - gl.glTexParameteri(gl.GL_TEXTURE_2D, gl.GL_TEXTURE_WRAP_S, gl.GL_CLAMP_TO_EDGE) - gl.glTexParameteri(gl.GL_TEXTURE_2D, gl.GL_TEXTURE_WRAP_T, gl.GL_CLAMP_TO_EDGE) - gl.glTexParameteri(gl.GL_TEXTURE_2D, gl.GL_TEXTURE_MIN_FILTER, gl.GL_LINEAR) - gl.glTexParameteri(gl.GL_TEXTURE_2D, gl.GL_TEXTURE_MAG_FILTER, gl.GL_LINEAR) + GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_S, GL.GL_CLAMP_TO_EDGE) + GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_WRAP_T, GL.GL_CLAMP_TO_EDGE) + GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MIN_FILTER, GL.GL_LINEAR) + GL.glTexParameteri(GL.GL_TEXTURE_2D, GL.GL_TEXTURE_MAG_FILTER, GL.GL_LINEAR) # Load the texture - gl.glTexImage2D( - gl.GL_TEXTURE_2D, + GL.glTexImage2D( + GL.GL_TEXTURE_2D, 0, - gl.GL_RGB32F, + GL.GL_RGB32F, width, height, 0, - gl.GL_RGB, - gl.GL_FLOAT, + GL.GL_RGB, + GL.GL_FLOAT, hdr_image, ) # Set up the viewport and projection - gl.glViewport( - 0, - 0, - glut.glutGet(glut.GLUT_WINDOW_WIDTH), - glut.glutGet(glut.GLUT_WINDOW_HEIGHT), - ) - gl.glMatrixMode(gl.GL_PROJECTION) - gl.glLoadIdentity() - glu.gluPerspective( - 45.0, - glut.glutGet(glut.GLUT_WINDOW_WIDTH) - / float(glut.glutGet(glut.GLUT_WINDOW_HEIGHT)), - 0.1, - 100.0, - ) - gl.glMatrixMode(gl.GL_MODELVIEW) - gl.glLoadIdentity() + window_width = GLUT.glutGet(GLUT.GLUT_WINDOW_WIDTH) + window_height = GLUT.glutGet(GLUT.GLUT_WINDOW_HEIGHT) + GL.glViewport(0, 0, window_width, window_height) + GL.glMatrixMode(GL.GL_PROJECTION) + GL.glLoadIdentity() + GLU.gluPerspective(45.0, window_width / float(window_height), 0.1, 100.0) + GL.glMatrixMode(GL.GL_MODELVIEW) + GL.glLoadIdentity() # Apply camera transformations - gl.glTranslatef(0.0, 0.0, -camera_distance) - gl.glRotatef(camera_angle_y, 1.0, 0.0, 0.0) - gl.glRotatef(camera_angle_x, 0.0, 1.0, 0.0) + GL.glTranslatef(0.0, 0.0, -STATE.distance) + GL.glRotatef(STATE.angle_y, 1.0, 0.0, 0.0) + GL.glRotatef(STATE.angle_x, 0.0, 1.0, 0.0) # Draw a textured sphere to simulate being inside the HDR environment - quadric = glu.gluNewQuadric() - glu.gluQuadricTexture(quadric, gl.GL_TRUE) - glu.gluSphere(quadric, 50.0, 50, 50) - glu.gluDeleteQuadric(quadric) + quadric = GLU.gluNewQuadric() + GLU.gluQuadricTexture(quadric, GL.GL_TRUE) + GLU.gluSphere(quadric, 50.0, 50, 50) + GLU.gluDeleteQuadric(quadric) # Disable texture mapping - gl.glDisable(gl.GL_TEXTURE_2D) + GL.glDisable(GL.GL_TEXTURE_2D) - glut.glutSwapBuffers() + GLUT.glutSwapBuffers() -def mouse_motion(x, y): - global mouse_last_x, mouse_last_y, camera_angle_x, camera_angle_y - if mouse_left_down: - dx = x - mouse_last_x - dy = y - mouse_last_y - camera_angle_x += dx * 0.1 - camera_angle_y += dy * 0.1 - mouse_last_x = x - mouse_last_y = y - glut.glutPostRedisplay() +def mouse_motion(x: int, y: int) -> None: + """Turn the camera by however far the pointer moved since the last event.""" + if STATE.left_down: + dx = x - STATE.last_x + dy = y - STATE.last_y + STATE.angle_x += dx * _DEGREES_PER_PIXEL + STATE.angle_y += dy * _DEGREES_PER_PIXEL + STATE.last_x = x + STATE.last_y = y + GLUT.glutPostRedisplay() -def mouse_button(button, state, x, y): - global mouse_left_down - if button == glut.GLUT_LEFT_BUTTON: - if state == glut.GLUT_DOWN: - mouse_left_down = True - mouse_last_x = x - mouse_last_y = y - elif state == glut.GLUT_UP: - mouse_left_down = False +def mouse_button(button: int, state: int, x: int, y: int) -> None: + """Start a drag on left-button press and end it on release.""" + if button == GLUT.GLUT_LEFT_BUTTON: + if state == GLUT.GLUT_DOWN: + STATE.left_down = True + # These two were assigned without a `global` declaration before, so + # they went to locals and were thrown away: the first motion event + # of every drag measured from wherever the previous drag ended. + STATE.last_x = x + STATE.last_y = y + elif state == GLUT.GLUT_UP: + STATE.left_down = False -def main(filepath): +def main(filepath: str) -> None: + """Load `filepath` and hand control to the GLUT main loop.""" try: width, height, hdr_map = load_hdr_environment_map(filepath) - print("HDR Map Loaded. Dimensions:", width, "x", height) + except OSError: + # Covers the not-found, not-readable and unreadable-EXR cases above. + # Anything else is a real fault and should surface as a traceback. + logger.exception("Could not load '%s'", filepath) + return - # Initialize GLUT and create window - glut.glutInit() - glut.glutInitDisplayMode(glut.GLUT_DOUBLE | glut.GLUT_RGB | glut.GLUT_DEPTH) - glut.glutInitWindowSize(800, 600) # Set initial window size - glut.glutCreateWindow(b"HDR Environment Map") + logger.info("HDR Map Loaded. Dimensions: %s x %s", width, height) - # Set display callback - glut.glutDisplayFunc(lambda: display_hdr_image(width, height, hdr_map)) + # Initialize GLUT and create window + GLUT.glutInit() + GLUT.glutInitDisplayMode(GLUT.GLUT_DOUBLE | GLUT.GLUT_RGB | GLUT.GLUT_DEPTH) + GLUT.glutInitWindowSize(800, 600) # Set initial window size + GLUT.glutCreateWindow(b"HDR Environment Map") - # Set mouse callbacks - glut.glutMotionFunc(mouse_motion) - glut.glutMouseFunc(mouse_button) + # Set display callback + GLUT.glutDisplayFunc(lambda: display_hdr_image(width, height, hdr_map)) - # Start the GLUT main loop - glut.glutMainLoop() - except Exception as e: - print(e) + # Set mouse callbacks + GLUT.glutMotionFunc(mouse_motion) + GLUT.glutMouseFunc(mouse_button) + + # Start the GLUT main loop + GLUT.glutMainLoop() # Example usage if __name__ == "__main__": - filepath = "lilienstein_4k.exr" - main(filepath) + logging.basicConfig(format="%(message)s", level=logging.INFO) + main("lilienstein_4k.exr") diff --git a/Programming/TRAK/code/photonmapping/python/geometry.py b/Programming/TRAK/code/photonmapping/python/geometry.py new file mode 100755 index 00000000..eca39f0f --- /dev/null +++ b/Programming/TRAK/code/photonmapping/python/geometry.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +"""Vectors, photons and the two surfaces the standalone photon mapper traces. + +Split out of `main.py` so each file stays under this repo's 250-line cap; the +classes are unchanged. +""" + +from __future__ import annotations + +import numpy as np + +# A ray whose direction is this close to parallel with a plane never hits it. +_PARALLEL_EPSILON = 1e-6 + +# Define basic vector operations +class Vector3: + """A 3D vector, with just the operations this renderer needs.""" + + def __init__(self, x: float, y: float, z: float) -> None: + """Store the three components.""" + self.x = x + self.y = y + self.z = z + + def __add__(self, other: Vector3) -> Vector3: + """Return the component-wise sum.""" + return Vector3(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other: Vector3) -> Vector3: + """Return the component-wise difference.""" + return Vector3(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar: float) -> Vector3: + """Return this vector scaled by `scalar`.""" + return Vector3(self.x * scalar, self.y * scalar, self.z * scalar) + + def dot(self, other: Vector3) -> float: + """Return the dot product with `other`.""" + return self.x * other.x + self.y * other.y + self.z * other.z + + def norm(self) -> float: + """Return the length of this vector.""" + return float(np.sqrt(self.dot(self))) + + def normalize(self) -> Vector3: + """Return a new vector of unit length pointing the same way.""" + n = self.norm() + return Vector3(self.x / n, self.y / n, self.z / n) + + +# Define the photon +class Photon: + """A packet of light energy travelling in a straight line.""" + + def __init__( + self, position: Vector3, direction: Vector3, power: np.ndarray + ) -> None: + """Store where the photon is, where it is going, and what it carries.""" + self.position = position + self.direction = direction + self.power = power + + +# Define a simple sphere +class Sphere: + """A sphere of uniform colour.""" + + def __init__(self, center: Vector3, radius: float, color: np.ndarray) -> None: + """Store the sphere's geometry and colour.""" + self.center = center + self.radius = radius + self.color = color + + def intersect( + self, ray_origin: Vector3, ray_direction: Vector3 + ) -> tuple[float, Vector3, Vector3] | None: + """Return the nearest hit as (distance, point, normal), or None.""" + # Solve quadratic equation for intersection + oc = ray_origin - self.center + a = ray_direction.dot(ray_direction) + b = 2.0 * oc.dot(ray_direction) + c = oc.dot(oc) - self.radius * self.radius + discriminant = b * b - 4 * a * c + if discriminant < 0: + return None # No intersection + t = (-b - np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + t = (-b + np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + return None + hit_point = ray_origin + ray_direction * t + normal = (hit_point - self.center).normalize() + return (t, hit_point, normal) + + +# Define a simple plane +class Plane: + """An infinite plane of uniform colour.""" + + def __init__(self, point: Vector3, normal: Vector3, color: np.ndarray) -> None: + """Store a point on the plane, its (normalised) normal, and its colour.""" + self.point = point + self.normal = normal.normalize() + self.color = color + + def intersect( + self, ray_origin: Vector3, ray_direction: Vector3 + ) -> tuple[float, Vector3, Vector3] | None: + """Return the nearest hit as (distance, point, normal), or None.""" + denom = self.normal.dot(ray_direction) + if abs(denom) > _PARALLEL_EPSILON: + t = (self.point - ray_origin).dot(self.normal) / denom + if t >= 0: + hit_point = ray_origin + ray_direction * t + return (t, hit_point, self.normal) + return None + + +Surface = Sphere | Plane diff --git a/Programming/TRAK/code/photonmapping/python/main.py b/Programming/TRAK/code/photonmapping/python/main.py old mode 100644 new mode 100755 index 99b99523..7c8dc610 --- a/Programming/TRAK/code/photonmapping/python/main.py +++ b/Programming/TRAK/code/photonmapping/python/main.py @@ -1,88 +1,32 @@ +#!/usr/bin/env python3 +"""Standalone photon mapper: scatter photons from the light, then gather them. + +The lab version, kept whole: the scene, the light and the photon map are +module-level state, and running the file emits photons, renders a 100x100 +image and shows it. `Programming/TRAK/photon_mapping.py` is the later +parameterised rewrite of the same algorithm. +""" + +from __future__ import annotations + +import logging import time +from typing import cast import matplotlib.pyplot as plt import numpy as np +from geometry import Photon, Plane, Sphere, Surface, Vector3 -PHOTONS = 0 +logger = logging.getLogger(__name__) +# Shadow rays start this far off the surface they were cast from. +_SHADOW_OFFSET = 1e-5 -# Define basic vector operations -class Vector3: - def __init__(self, x, y, z): - self.x = x - self.y = y - self.z = z - - def __add__(self, other): - return Vector3(self.x + other.x, self.y + other.y, self.z + other.z) - - def __sub__(self, other): - return Vector3(self.x - other.x, self.y - other.y, self.z - other.z) - - def __mul__(self, scalar): - return Vector3(self.x * scalar, self.y * scalar, self.z * scalar) - - def dot(self, other): - return self.x * other.x + self.y * other.y + self.z * other.z - - def norm(self): - return np.sqrt(self.dot(self)) - - def normalize(self): - n = self.norm() - return Vector3(self.x / n, self.y / n, self.z / n) - - -# Define the photon -class Photon: - def __init__(self, position, direction, power): - self.position = position - self.direction = direction - self.power = power - - -# Define a simple sphere -class Sphere: - def __init__(self, center, radius, color): - self.center = center - self.radius = radius - self.color = color - - def intersect(self, ray_origin, ray_direction): - # Solve quadratic equation for intersection - oc = ray_origin - self.center - a = ray_direction.dot(ray_direction) - b = 2.0 * oc.dot(ray_direction) - c = oc.dot(oc) - self.radius * self.radius - discriminant = b * b - 4 * a * c - if discriminant < 0: - return None # No intersection - t = (-b - np.sqrt(discriminant)) / (2.0 * a) - if t < 0: - t = (-b + np.sqrt(discriminant)) / (2.0 * a) - if t < 0: - return None - hit_point = ray_origin + ray_direction * t - normal = (hit_point - self.center).normalize() - return (t, hit_point, normal) - - -# Define a simple plane -class Plane: - def __init__(self, point, normal, color): - self.point = point - self.normal = normal.normalize() - self.color = color - - def intersect(self, ray_origin, ray_direction): - denom = self.normal.dot(ray_direction) - if abs(denom) > 1e-6: - t = (self.point - ray_origin).dot(self.normal) / denom - if t >= 0: - hit_point = ray_origin + ray_direction * t - return (t, hit_point, self.normal) - return None +# Fraction of its power a photon keeps across one diffuse bounce. +_BOUNCE_ABSORPTION = 0.8 +# One generator for the whole module, matching the module-level scene below. +_rng = np.random.default_rng() # Scene setup sphere = Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0])) # Red sphere @@ -90,14 +34,14 @@ plane = Plane( Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5]) ) # Gray plane -objects = [sphere, plane] +objects: list[Surface] = [sphere, plane] # Light source light_position = Vector3(-5, 5, -5) light_power = np.array([1, 1, 1]) * 1000 # Intense white light # Photon map -photon_map = [] +photon_map: list[tuple[Vector3, np.ndarray]] = [] # Parameters num_photons = 10000 # Number of photons to emit @@ -105,75 +49,26 @@ max_depth = 5 # Maximum number of bounces gather_radius = 0.5 # Radius for radiance estimation -def emit_photons(): - for _ in range(num_photons): - # Emit photons in random directions from the light source - direction = random_unit_vector() - power = light_power / num_photons - photon = Photon(light_position, direction, power) - trace_photon(photon, 0) - - -def trace_photon(photon, depth): - global PHOTONS - PHOTONS += 1 - if depth > max_depth: - return - closest_t = np.inf - hit_object = None - hit_info = None - # Find the nearest intersection - for obj in objects: - result = obj.intersect(photon.position, photon.direction) - if result: - t, hit_point, normal = result - if t < closest_t: - closest_t = t - hit_object = obj - hit_info = (hit_point, normal) - if hit_object: - hit_point, normal = hit_info - photon_map.append((hit_point, photon.power)) - # Diffuse reflection - new_direction = random_hemisphere_direction(normal) - photon.position = hit_point - photon.direction = new_direction - # Absorb some power - photon.power = photon.power * 0.8 # Simple absorption - trace_photon(photon, depth + 1) - - -def random_unit_vector(): - theta = np.random.uniform(0, 2 * np.pi) - z = np.random.uniform(-1, 1) +def random_unit_vector() -> Vector3: + """Draw a direction uniformly over the sphere.""" + theta = _rng.uniform(0, 2 * np.pi) + z = _rng.uniform(-1, 1) r = np.sqrt(1 - z * z) return Vector3(r * np.cos(theta), r * np.sin(theta), z) -def random_hemisphere_direction(normal): - dir = random_unit_vector() - if dir.dot(normal) < 0: - dir = Vector3(-dir.x, -dir.y, -dir.z) - return dir +def random_hemisphere_direction(normal: Vector3) -> Vector3: + """Draw a direction uniformly over the hemisphere around `normal`.""" + direction = random_unit_vector() + if direction.dot(normal) < 0: + direction = Vector3(-direction.x, -direction.y, -direction.z) + return direction -def render_image(width, height): - aspect_ratio = width / height - fov = np.pi / 3 # 60 degrees field of view - image = np.zeros((height, width, 3)) - for y in range(height): - for x in range(width): - # Convert pixel coordinate to camera ray - px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio - py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2) - ray_origin = Vector3(0, 0, 0) - ray_direction = Vector3(px, py, -1).normalize() - color = trace_ray(ray_origin, ray_direction) - image[y, x, :] = np.clip(color, 0, 1) - return image - - -def trace_ray(ray_origin, ray_direction): +def nearest_hit( + ray_origin: Vector3, ray_direction: Vector3 +) -> tuple[Surface, Vector3, Vector3] | None: + """Return the closest surface the ray meets, with its point and normal.""" closest_t = np.inf hit_object = None hit_info = None @@ -185,35 +80,59 @@ def trace_ray(ray_origin, ray_direction): if t < closest_t: closest_t = t hit_object = obj - hit_info = (hit_point, normal, obj.color) - if hit_object: - hit_point, normal, color = hit_info - direct_light = compute_direct_light(hit_point, normal) - indirect_light = estimate_radiance(hit_point, normal) - return color * (direct_light + indirect_light) - return np.array([0, 0, 0]) # Background color + hit_info = (hit_point, normal) + if hit_object is None or hit_info is None: + return None + return (hit_object, hit_info[0], hit_info[1]) -def compute_direct_light(point, normal): +def trace_photon(photon: Photon, depth: int) -> int: + """Bounce one photon, recording each landing, and return the traces made. + + The count is what the original `global PHOTONS` accumulated: one per call, + bounces included. + """ + if depth > max_depth: + return 1 + hit = nearest_hit(photon.position, photon.direction) + if hit is None: + return 1 + _, hit_point, normal = hit + photon_map.append((hit_point, photon.power)) + # Diffuse reflection + photon.position = hit_point + photon.direction = random_hemisphere_direction(normal) + # Absorb some power + photon.power = photon.power * _BOUNCE_ABSORPTION # Simple absorption + return 1 + trace_photon(photon, depth + 1) + + +def emit_photons() -> int: + """Emit `num_photons` from the light and return the total traces made.""" + traces = 0 + for _ in range(num_photons): + # Emit photons in random directions from the light source + direction = random_unit_vector() + power = light_power / num_photons + traces += trace_photon(Photon(light_position, direction, power), 0) + return traces + + +def compute_direct_light(point: Vector3, normal: Vector3) -> np.ndarray: + """Return the light arriving straight from the source, or black if shadowed.""" # Simple Lambertian reflection from light source direction_to_light = (light_position - point).normalize() # Shadow ray - shadow_origin = point + normal * 1e-5 - shadow_ray = direction_to_light - in_shadow = False + shadow_origin = point + normal * _SHADOW_OFFSET for obj in objects: - result = obj.intersect(shadow_origin, shadow_ray) - if result: - in_shadow = True - break - if in_shadow: - return np.array([0, 0, 0]) + if obj.intersect(shadow_origin, direction_to_light): + return np.array([0, 0, 0]) intensity = max(0, normal.dot(direction_to_light)) return intensity * light_power / (4 * np.pi * (light_position - point).norm() ** 2) -def estimate_radiance(point, normal): - # Gather photons within the gather_radius +def estimate_radiance(point: Vector3, normal: Vector3) -> np.ndarray: + """Gather the photons within `gather_radius` of `point` into a radiance.""" accumulated_power = np.array([0.0, 0.0, 0.0]) for photon_pos, photon_power in photon_map: distance = (photon_pos - point).norm() @@ -224,27 +143,62 @@ def estimate_radiance(point, normal): return accumulated_power / (area * num_photons) -# Main execution -if __name__ == "__main__": - print("Emitting photons...") - emit_photons() +def trace_ray(ray_origin: Vector3, ray_direction: Vector3) -> np.ndarray: + """Shade one camera ray with direct light plus gathered photons.""" + hit = nearest_hit(ray_origin, ray_direction) + if hit is None: + return np.array([0, 0, 0]) # Background color + hit_object, hit_point, normal = hit + direct_light = compute_direct_light(hit_point, normal) + indirect_light = estimate_radiance(hit_point, normal) + return cast("np.ndarray", hit_object.color * (direct_light + indirect_light)) - print("Rendering image...") + +def render_image(width: int, height: int) -> np.ndarray: + """Shoot one camera ray per pixel and return the `height` x `width` image.""" + aspect_ratio = width / height + fov = np.pi / 3 # 60 degrees field of view + image = np.zeros((height, width, 3)) + for y in range(height): + for x in range(width): + # Convert pixel coordinate to camera ray + px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio + py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2) + ray_direction = Vector3(px, py, -1).normalize() + image[y, x, :] = np.clip(trace_ray(Vector3(0, 0, 0), ray_direction), 0, 1) + return image + + +def main() -> None: + """Emit photons, render the image and show it.""" + logger.info("Emitting photons...") + photon_traces = emit_photons() + + logger.info("Rendering image...") width = 100 height = 100 t0 = time.time() image = render_image(width, height) - print( - f"Render Took: {round(time.time() - t0, 2)}s\n" - f"resolution: {width}x{height}\n" - f"photons (emitted): {num_photons}\n" - f"photons (reflected): {PHOTONS - num_photons}\n" - f"photons (total): {PHOTONS}\n" - f"rays: {width * height}" + logger.info( + "Render Took: %ss\nresolution: %sx%s\nphotons (emitted): %s\n" + "photons (reflected): %s\nphotons (total): %s\nrays: %s", + round(time.time() - t0, 2), + width, + height, + num_photons, + photon_traces - num_photons, + photon_traces, + width * height, ) # Display the image plt.imshow(image) plt.axis("off") plt.show() + + +# Main execution +if __name__ == "__main__": + logging.basicConfig(format="%(message)s", level=logging.INFO) + main() diff --git a/Programming/TRAK/lab6.py b/Programming/TRAK/lab6.py old mode 100644 new mode 100755 index 0fe9a6a8..b3b07c57 --- a/Programming/TRAK/lab6.py +++ b/Programming/TRAK/lab6.py @@ -1,2 +1,11 @@ +#!/usr/bin/env python3 +"""Lab 6 driver: run another script in this directory by path. + +`script_path` was filled in by hand during the lab, which is why it is empty +here -- the file is kept as submitted. +""" + +import runpy + script_path = "" # path to script -exec(compile(open(script_path).read(), script_path, "exec")) +runpy.run_path(script_path, run_name="__main__") diff --git a/Programming/TRAK/main.py b/Programming/TRAK/main.py old mode 100644 new mode 100755 index 9aa9c4bc..d46375d1 --- a/Programming/TRAK/main.py +++ b/Programming/TRAK/main.py @@ -1,14 +1,25 @@ +#!/usr/bin/env python3 +"""Entry point: render a scene with one of the three algorithms in this course. + +The algorithm, scene, resolution and output path all default to `config.ini` +and are overridable on the command line. +""" + import argparse import importlib -import os +import logging +from pathlib import Path import matplotlib.pyplot as plt from photon_mapping import render_photon_mapping from rendering import ray_trace from utils import load_config, parse_resolution +logger = logging.getLogger(__name__) -def main(): + +def main() -> None: + """Parse the command line over `config.ini` defaults and run the renderer.""" # default config config = load_config("config.ini") @@ -98,58 +109,64 @@ def main(): # Run the selected algorithm if args.algorithm == "ray_tracing0": - print("Starting ray tracing zero...") - # ray_trace(args.scene, args.environment_map, image_width=width, image_height=height, output_file="output_ray_traced.png") + logger.info("Starting ray tracing zero...") ray_trace( args.num_spheres, - args.environment_map, + args.environment, image_width=width, - image_height=height, # na razie generujemy w kodzie, ale potem trzeba będzie obj wczytywać - output_file="output_ray_traced.png", + # na razie generujemy w kodzie, ale potem trzeba będzie obj wczytywać + image_height=height, + output_file=Path("outputs") / args.output, ) elif args.algorithm == "ray_tracing": - print("Starting ray tracing...") + logger.info("Starting ray tracing...") try: - print(args.scene) + logger.info("%s", args.scene) scene_module = importlib.import_module(f"scenes.{args.scene}") except ModuleNotFoundError: - print(f"Error: Scene '{args.scene}' not found in the 'scenes' directory.") + logger.info( + "Error: Scene '%s' not found in the 'scenes' directory.", args.scene + ) return try: scene = scene_module.setup_scene( width=width, height=height, environment=f"{args.environment}" ) except AttributeError: - print( - f"Error: Scene '{args.scene}' does not define a `setup_scene` function." + logger.info( + "Error: Scene '%s' does not define a `setup_scene` function.", + args.scene, ) return # Renderowanie - print( - f"Rendering scene '{args.scene}' with {args.samples_per_pixel} samples per pixel..." + logger.info( + "Rendering scene '%s' with %s samples per pixel...", + args.scene, + args.samples_per_pixel, ) img = scene.render(samples_per_pixel=args.samples_per_pixel) - output_path = os.path.join("outputs", args.output) + output_path = Path("outputs") / args.output img.save(output_path) - print(f"Image saved to {output_path}") + logger.info("Image saved to %s", output_path) if args.show: img.show() elif args.algorithm == "photon_mapping": - print("Starting photon mapping...") + logger.info("Starting photon mapping...") image = render_photon_mapping( width, height, args.num_photons, args.max_depth, args.gather_radius ) plt.imshow(image) plt.axis("off") - output_path = os.path.join("outputs", args.output) + output_path = Path("outputs") / args.output plt.savefig(output_path) - print(f"Image saved to {output_path}") + logger.info("Image saved to %s", output_path) if args.show: plt.show() else: - print(f"Unknown algorithm: {args.algorithm}") + logger.info("Unknown algorithm: %s", args.algorithm) return if __name__ == "__main__": + logging.basicConfig(format="%(message)s", level=logging.INFO) main() diff --git a/Programming/TRAK/photon_geometry.py b/Programming/TRAK/photon_geometry.py new file mode 100755 index 00000000..9ab48b8c --- /dev/null +++ b/Programming/TRAK/photon_geometry.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python3 +"""Vectors, photons and the two surfaces the photon mapper traces. + +Split out of `photon_mapping.py` so each file stays under this repo's +250-line cap; the classes are unchanged. +""" + +from __future__ import annotations + +import numpy as np + +# A ray whose direction is this close to parallel with a plane never hits it. +_PARALLEL_EPSILON = 1e-6 + +class Vector3: + """A 3D vector, with just the operations this renderer needs.""" + + def __init__(self, x: float, y: float, z: float) -> None: + """Store the three components.""" + self.x = x + self.y = y + self.z = z + + def __add__(self, other: Vector3) -> Vector3: + """Return the component-wise sum.""" + return Vector3(self.x + other.x, self.y + other.y, self.z + other.z) + + def __sub__(self, other: Vector3) -> Vector3: + """Return the component-wise difference.""" + return Vector3(self.x - other.x, self.y - other.y, self.z - other.z) + + def __mul__(self, scalar: float) -> Vector3: + """Return this vector scaled by `scalar`.""" + return Vector3(self.x * scalar, self.y * scalar, self.z * scalar) + + def dot(self, other: Vector3) -> float: + """Return the dot product with `other`.""" + return self.x * other.x + self.y * other.y + self.z * other.z + + def norm(self) -> float: + """Return the length of this vector.""" + return float(np.sqrt(self.dot(self))) + + def normalize(self) -> Vector3: + """Return a new vector of unit length pointing the same way.""" + n = self.norm() + return Vector3(self.x / n, self.y / n, self.z / n) + + +class Photon: + """A packet of light energy travelling in a straight line.""" + + def __init__( + self, position: Vector3, direction: Vector3, power: np.ndarray + ) -> None: + """Store where the photon is, where it is going, and what it carries.""" + self.position = position + self.direction = direction + self.power = power + + +class Sphere: + """A sphere of uniform colour.""" + + def __init__(self, center: Vector3, radius: float, color: np.ndarray) -> None: + """Store the sphere's geometry and colour.""" + self.center = center + self.radius = radius + self.color = color + + def intersect( + self, ray_origin: Vector3, ray_direction: Vector3 + ) -> tuple[float, Vector3, Vector3] | None: + """Return the nearest hit as (distance, point, normal), or None.""" + oc = ray_origin - self.center + a = ray_direction.dot(ray_direction) + b = 2.0 * oc.dot(ray_direction) + c = oc.dot(oc) - self.radius * self.radius + discriminant = b * b - 4 * a * c + if discriminant < 0: + return None + t = (-b - np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + t = (-b + np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + return None + hit_point = ray_origin + ray_direction * t + normal = (hit_point - self.center).normalize() + return (t, hit_point, normal) + + +class Plane: + """An infinite plane of uniform colour.""" + + def __init__(self, point: Vector3, normal: Vector3, color: np.ndarray) -> None: + """Store a point on the plane, its (normalised) normal, and its colour.""" + self.point = point + self.normal = normal.normalize() + self.color = color + + def intersect( + self, ray_origin: Vector3, ray_direction: Vector3 + ) -> tuple[float, Vector3, Vector3] | None: + """Return the nearest hit as (distance, point, normal), or None.""" + denom = self.normal.dot(ray_direction) + if abs(denom) > _PARALLEL_EPSILON: + t = (self.point - ray_origin).dot(self.normal) / denom + if t >= 0: + hit_point = ray_origin + ray_direction * t + return (t, hit_point, self.normal) + return None + + +Surface = Sphere | Plane diff --git a/Programming/TRAK/photon_mapping.py b/Programming/TRAK/photon_mapping.py old mode 100644 new mode 100755 index 1719b025..59b1f6c9 --- a/Programming/TRAK/photon_mapping.py +++ b/Programming/TRAK/photon_mapping.py @@ -1,181 +1,150 @@ +#!/usr/bin/env python3 +"""A two-pass photon mapper: scatter photons from the light, then gather them. + +The first pass emits `num_photons` from the light and records every surface +they land on. The second pass shoots one camera ray per pixel and shades each +hit with direct light plus the photons gathered within `gather_radius`. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import cast + import numpy as np +from photon_geometry import Photon, Plane, Sphere, Surface, Vector3 + +# Shadow rays start this far off the surface they were cast from. +_SHADOW_OFFSET = 1e-5 + +# Fraction of its power a photon keeps across one diffuse bounce. +_BOUNCE_ABSORPTION = 0.8 -# Define basic vector operations -class Vector3: - def __init__(self, x, y, z): - self.x = x - self.y = y - self.z = z - - def __add__(self, other): - return Vector3(self.x + other.x, self.y + other.y, self.z + other.z) - - def __sub__(self, other): - return Vector3(self.x - other.x, self.y - other.y, self.z - other.z) - - def __mul__(self, scalar): - return Vector3(self.x * scalar, self.y * scalar, self.z * scalar) - - def dot(self, other): - return self.x * other.x + self.y * other.y + self.z * other.z - - def norm(self): - return np.sqrt(self.dot(self)) - - def normalize(self): - n = self.norm() - return Vector3(self.x / n, self.y / n, self.z / n) -class Photon: - def __init__(self, position, direction, power): - self.position = position - self.direction = direction - self.power = power +@dataclass +class Scene: + """Everything the two passes share: the geometry, the light, the photons.""" + + objects: list[Surface] + rng: np.random.Generator + light_position: Vector3 + light_power: np.ndarray + gather_radius: float + photon_map: list[tuple[Vector3, np.ndarray]] -class Sphere: - def __init__(self, center, radius, color): - self.center = center - self.radius = radius - self.color = color - - def intersect(self, ray_origin, ray_direction): - oc = ray_origin - self.center - a = ray_direction.dot(ray_direction) - b = 2.0 * oc.dot(ray_direction) - c = oc.dot(oc) - self.radius * self.radius - discriminant = b * b - 4 * a * c - if discriminant < 0: - return None - t = (-b - np.sqrt(discriminant)) / (2.0 * a) - if t < 0: - t = (-b + np.sqrt(discriminant)) / (2.0 * a) - if t < 0: - return None - hit_point = ray_origin + ray_direction * t - normal = (hit_point - self.center).normalize() - return (t, hit_point, normal) +def random_unit_vector(rng: np.random.Generator) -> Vector3: + """Draw a direction uniformly over the sphere.""" + theta = rng.uniform(0, 2 * np.pi) + z = rng.uniform(-1, 1) + r = np.sqrt(1 - z * z) + return Vector3(r * np.cos(theta), r * np.sin(theta), z) -class Plane: - def __init__(self, point, normal, color): - self.point = point - self.normal = normal.normalize() - self.color = color +def random_hemisphere_direction(rng: np.random.Generator, normal: Vector3) -> Vector3: + """Draw a direction uniformly over the hemisphere around `normal`.""" + direction = random_unit_vector(rng) + if direction.dot(normal) < 0: + direction = Vector3(-direction.x, -direction.y, -direction.z) + return direction - def intersect(self, ray_origin, ray_direction): - denom = self.normal.dot(ray_direction) - if abs(denom) > 1e-6: - t = (self.point - ray_origin).dot(self.normal) / denom - if t >= 0: - hit_point = ray_origin + ray_direction * t - return (t, hit_point, self.normal) + +def nearest_hit( + objects: list[Surface], ray_origin: Vector3, ray_direction: Vector3 +) -> tuple[Surface, Vector3, Vector3] | None: + """Return the closest surface the ray meets, with its point and normal.""" + closest_t = np.inf + hit_object = None + hit_info = None + for obj in objects: + result = obj.intersect(ray_origin, ray_direction) + if result: + t, hit_point, normal = result + if t < closest_t: + closest_t = t + hit_object = obj + hit_info = (hit_point, normal) + if hit_object is None or hit_info is None: return None + return (hit_object, hit_info[0], hit_info[1]) -def render_photon_mapping(width, height, num_photons, max_depth, gather_radius): - # Photon mapping logic - photon_map = [] - sphere = Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0])) # Red sphere - plane = Plane( - Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5]) - ) # Gray plane - objects = [sphere, plane] +def trace_photon(scene: Scene, photon: Photon, depth: int, max_depth: int) -> None: + """Bounce one photon through the scene, recording every surface it lands on.""" + if depth > max_depth: + return + hit = nearest_hit(scene.objects, photon.position, photon.direction) + if hit is not None: + _, hit_point, normal = hit + scene.photon_map.append((hit_point, photon.power)) + photon.position = hit_point + photon.direction = random_hemisphere_direction(scene.rng, normal) + photon.power = photon.power * _BOUNCE_ABSORPTION + trace_photon(scene, photon, depth + 1, max_depth) - light_position = Vector3(-5, 5, -5) - light_power = np.array([1, 1, 1]) * 1000 - def emit_photons(): - for _ in range(num_photons): - direction = random_unit_vector() - power = light_power / num_photons - photon = Photon(light_position, direction, power) - trace_photon(photon, 0) +def emit_photons(scene: Scene, num_photons: int, max_depth: int) -> None: + """Run the scatter pass: emit `num_photons` from the light and trace each.""" + for _ in range(num_photons): + direction = random_unit_vector(scene.rng) + power = scene.light_power / num_photons + photon = Photon(scene.light_position, direction, power) + trace_photon(scene, photon, 0, max_depth) - def trace_photon(photon, depth): - if depth > max_depth: - return - closest_t = np.inf - hit_object = None - hit_info = None - for obj in objects: - result = obj.intersect(photon.position, photon.direction) - if result: - t, hit_point, normal = result - if t < closest_t: - closest_t = t - hit_object = obj - hit_info = (hit_point, normal) - if hit_object: - hit_point, normal = hit_info - photon_map.append((hit_point, photon.power)) - new_direction = random_hemisphere_direction(normal) - photon.position = hit_point - photon.direction = new_direction - photon.power = photon.power * 0.8 - trace_photon(photon, depth + 1) - def random_unit_vector(): - theta = np.random.uniform(0, 2 * np.pi) - z = np.random.uniform(-1, 1) - r = np.sqrt(1 - z * z) - return Vector3(r * np.cos(theta), r * np.sin(theta), z) - - def random_hemisphere_direction(normal): - dir = random_unit_vector() - if dir.dot(normal) < 0: - dir = Vector3(-dir.x, -dir.y, -dir.z) - return dir - - def trace_ray(ray_origin, ray_direction): - closest_t = np.inf - hit_object = None - hit_info = None - for obj in objects: - result = obj.intersect(ray_origin, ray_direction) - if result: - t, hit_point, normal = result - if t < closest_t: - closest_t = t - hit_object = obj - hit_info = (hit_point, normal, obj.color) - if hit_object: - hit_point, normal, color = hit_info - direct_light = compute_direct_light(hit_point, normal) - indirect_light = estimate_radiance(hit_point, normal) - return color * (direct_light + indirect_light) - return np.array([0, 0, 0]) - - def compute_direct_light(point, normal): - direction_to_light = (light_position - point).normalize() - shadow_origin = point + normal * 1e-5 - shadow_ray = direction_to_light - in_shadow = False - for obj in objects: - result = obj.intersect(shadow_origin, shadow_ray) - if result: - in_shadow = True - break - if in_shadow: +def compute_direct_light(scene: Scene, point: Vector3, normal: Vector3) -> np.ndarray: + """Return the light arriving straight from the source, or black if shadowed.""" + direction_to_light = (scene.light_position - point).normalize() + shadow_origin = point + normal * _SHADOW_OFFSET + for obj in scene.objects: + if obj.intersect(shadow_origin, direction_to_light): return np.array([0, 0, 0]) - intensity = max(0, normal.dot(direction_to_light)) - return ( - intensity * light_power / (4 * np.pi * (light_position - point).norm() ** 2) - ) + intensity = max(0, normal.dot(direction_to_light)) + distance_squared = (scene.light_position - point).norm() ** 2 + return intensity * scene.light_power / (4 * np.pi * distance_squared) - def estimate_radiance(point, normal): - accumulated_power = np.array([0.0, 0.0, 0.0]) - for photon_pos, photon_power in photon_map: - distance = (photon_pos - point).norm() - if distance < gather_radius: - weight = max(0, normal.dot((photon_pos - point).normalize())) - accumulated_power += photon_power * weight - area = np.pi * gather_radius**2 - return accumulated_power / area - emit_photons() +def estimate_radiance(scene: Scene, point: Vector3, normal: Vector3) -> np.ndarray: + """Gather the photons within `gather_radius` of `point` into a radiance.""" + accumulated_power = np.array([0.0, 0.0, 0.0]) + for photon_pos, photon_power in scene.photon_map: + distance = (photon_pos - point).norm() + if distance < scene.gather_radius: + weight = max(0, normal.dot((photon_pos - point).normalize())) + accumulated_power += photon_power * weight + area = np.pi * scene.gather_radius**2 + return accumulated_power / area + + +def trace_ray(scene: Scene, ray_origin: Vector3, ray_direction: Vector3) -> np.ndarray: + """Shade one camera ray with direct light plus gathered photons.""" + hit = nearest_hit(scene.objects, ray_origin, ray_direction) + if hit is None: + return np.array([0, 0, 0]) + hit_object, hit_point, normal = hit + direct_light = compute_direct_light(scene, hit_point, normal) + indirect_light = estimate_radiance(scene, hit_point, normal) + return cast("np.ndarray", hit_object.color * (direct_light + indirect_light)) + + +def render_photon_mapping( + width: int, height: int, num_photons: int, max_depth: int, gather_radius: float +) -> np.ndarray: + """Render the fixed red-sphere-over-grey-plane scene and return the image.""" + scene = Scene( + objects=[ + Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0])), # Red sphere + Plane(Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5])), + ], + rng=np.random.default_rng(), + light_position=Vector3(-5, 5, -5), + light_power=np.array([1, 1, 1]) * 1000, + gather_radius=gather_radius, + photon_map=[], + ) + emit_photons(scene, num_photons, max_depth) aspect_ratio = width / height fov = np.pi / 3 @@ -184,9 +153,8 @@ def render_photon_mapping(width, height, num_photons, max_depth, gather_radius): for x in range(width): px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2) - ray_origin = Vector3(0, 0, 0) ray_direction = Vector3(px, py, -1).normalize() - color = trace_ray(ray_origin, ray_direction) + color = trace_ray(scene, Vector3(0, 0, 0), ray_direction) image[y, x, :] = np.clip(color, 0, 1) return image diff --git a/Programming/TRAK/raytracer.py b/Programming/TRAK/raytracer.py new file mode 100755 index 00000000..26d99bcc --- /dev/null +++ b/Programming/TRAK/raytracer.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python3 +"""Ray/geometry intersection and Phong shading for the `ray_tracing0` renderer. + +These were nested closures inside `rendering.ray_trace`, which is what made +that one function 460 lines long. They are pure now: the scene and the light +constants are passed in rather than captured, so each can be read on its own. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any, cast + +import numpy as np + +# A scene object is a plain dict tagged with "type": "sphere" or "plane". +SceneObject = dict[str, Any] + +# A ray whose direction is this close to parallel with a plane never hits it. +_PARALLEL_EPSILON = 1e-6 + +# Shadow and reflection rays start this far off the surface, so that the +# surface they just left does not intersect them at t=0. +SURFACE_OFFSET = 0.0001 + + +@dataclass(frozen=True) +class Lighting: + """The single point light and the Phong coefficients every surface uses.""" + + position: np.ndarray + color: np.ndarray + ambient: float + diffuse_c: float + specular_c: float + specular_k: int + + +def normalize(vector: np.ndarray) -> np.ndarray: + """Scale `vector` to unit length, in place, and return it.""" + vector /= np.linalg.norm(vector) + return vector + + +def intersect_plane( + ray_origin: np.ndarray, + ray_direction: np.ndarray, + plane_point: np.ndarray, + plane_normal: np.ndarray, +) -> float: + """Distance along the ray to the plane, or +inf if it never reaches it.""" + denom = np.dot(ray_direction, plane_normal) + if np.abs(denom) < _PARALLEL_EPSILON: + return np.inf + d = np.dot(plane_point - ray_origin, plane_normal) / denom + if d < 0: + return np.inf + return float(d) + + +def calculate_sphere_intersection(a: float, b: float, c: float, disc: float) -> float: + """Solve `a*t^2 + b*t + c` for the nearest root ahead of the ray origin. + + `disc` is the discriminant, passed in because the caller already has it. + The root is taken via the numerically stable `q` form rather than the + textbook quadratic formula, which loses precision when `b` dominates. + """ + if disc > 0: + distance_squared = np.sqrt(disc) + q = (-b - distance_squared) / 2.0 if b < 0 else (-b + distance_squared) / 2.0 + t0 = q / a + t1 = c / q + t0, t1 = min(t0, t1), max(t0, t1) + if t1 >= 0: + return float(t1) if t0 < 0 else float(t0) + return np.inf + + +def intersect_sphere( + ray_origin: np.ndarray, + ray_direction: np.ndarray, + sphere_center: np.ndarray, + sphere_radius: float, +) -> float: + """Distance along the ray to the sphere, or +inf if it misses.""" + a = np.dot(ray_direction, ray_direction) + origin_to_center = ray_origin - sphere_center + b = 2 * np.dot(ray_direction, origin_to_center) + radius_squared = sphere_radius * sphere_radius + c = np.dot(origin_to_center, origin_to_center) - radius_squared + disc = b * b - 4 * a * c + return calculate_sphere_intersection(a, b, c, disc) + + +def intersect( + ray_origin: np.ndarray, ray_direction: np.ndarray, object_: SceneObject +) -> float: + """Dispatch to the plane or sphere intersection for `object_`.""" + if object_["type"] == "plane": + return intersect_plane( + ray_origin, ray_direction, object_["position"], object_["normal"] + ) + return intersect_sphere( + ray_origin, ray_direction, object_["position"], object_["radius"] + ) + + +def get_normal(object_: SceneObject, intersection_point: np.ndarray) -> np.ndarray: + """Surface normal of `object_` at `intersection_point`.""" + if object_["type"] == "sphere": + return normalize(intersection_point - object_["position"]) + if object_["type"] == "plane": + return cast("np.ndarray", object_["normal"]) + msg = f"Unknown object type: {object_['type']}" + raise ValueError(msg) + + +def get_color(object_: SceneObject, intersection_point: np.ndarray) -> np.ndarray: + """Return the colour of `object_`: either a constant or a function of position.""" + color = object_["color"] + if not hasattr(color, "__len__"): + color = color(intersection_point) + return cast("np.ndarray", color) + + +def find_intersection( + scene: list[SceneObject], ray_origin: np.ndarray, ray_direction: np.ndarray +) -> tuple[float, int]: + """Return the distance to the nearest hit and that object's scene index. + + The index is -1 when the ray hits nothing, in which case the distance is + +inf and callers stop. + """ + t = np.inf + obj_idx = -1 + for index, object_ in enumerate(scene): + t_obj = intersect(ray_origin, ray_direction, object_) + if t_obj < t: + t, obj_idx = t_obj, index + return t, obj_idx + + +def is_shadowed( + scene: list[SceneObject], + lighting: Lighting, + intersection_point: np.ndarray, + normal: np.ndarray, + obj_idx: int, +) -> bool: + """Report whether anything but object `obj_idx` blocks the light here.""" + to_light = normalize(lighting.position - intersection_point) + shadow_intersections = [ + intersect(intersection_point + normal * SURFACE_OFFSET, to_light, obj_sh) + for k, obj_sh in enumerate(scene) + if k != obj_idx + ] + return bool(shadow_intersections) and min(shadow_intersections) < np.inf + + +def compute_color( + lighting: Lighting, + object_: SceneObject, + intersection_point: np.ndarray, + normal: np.ndarray, + ray_origin: np.ndarray, +) -> tuple[SceneObject, np.ndarray, np.ndarray, np.ndarray]: + """Shade the hit with Blinn-Phong and return it alongside the hit geometry.""" + color = get_color(object_, intersection_point) + to_light = normalize(lighting.position - intersection_point) + to_origin = normalize(ray_origin - intersection_point) + color_ray = lighting.ambient + diffuse_intensity = object_.get("diffuse_c", lighting.diffuse_c) * max( + np.dot(normal, to_light), 0 + ) + color_ray += diffuse_intensity * color + half_vector = normalize(to_light + to_origin) + specular_intensity = ( + object_.get("specular_c", lighting.specular_c) + * max(np.dot(normal, half_vector), 0) ** lighting.specular_k + ) + color_ray += specular_intensity * lighting.color + return object_, intersection_point, normal, color_ray + + +def trace_ray( + scene: list[SceneObject], + lighting: Lighting, + ray_origin: np.ndarray, + ray_direction: np.ndarray, +) -> tuple[SceneObject, np.ndarray, np.ndarray, np.ndarray] | None: + """Trace one ray and return its hit, or None if it misses or lands in shadow.""" + t, obj_idx = find_intersection(scene, ray_origin, ray_direction) + if t == np.inf: + return None + object_ = scene[obj_idx] + intersection_point = ray_origin + ray_direction * t + normal = get_normal(object_, intersection_point) + if is_shadowed(scene, lighting, intersection_point, normal, obj_idx): + return None + return compute_color(lighting, object_, intersection_point, normal, ray_origin) + + +def make_sphere( + position: list[float], radius: float, color: np.ndarray +) -> SceneObject: + """Build a sphere scene object.""" + return { + "type": "sphere", + "position": np.array(position), + "radius": np.array(radius), + "color": np.array(color), + "reflection": 0.5, + } + + +def make_plane( + position: list[float], + normal: list[float], + color_a: np.ndarray, + color_b: np.ndarray, +) -> SceneObject: + """Build a plane scene object with a checkerboard of `color_a`/`color_b`.""" + return { + "type": "plane", + "position": np.array(position), + "normal": np.array(normal), + "color": lambda point: ( + color_a if (int(point[0] * 2) % 2) == (int(point[2] * 2) % 2) else color_b + ), + "diffuse_c": 0.75, + "specular_c": 0.5, + "reflection": 0.25, + } diff --git a/Programming/TRAK/rendering.py b/Programming/TRAK/rendering.py old mode 100644 new mode 100755 index f9b23230..aa12a272 --- a/Programming/TRAK/rendering.py +++ b/Programming/TRAK/rendering.py @@ -1,382 +1,52 @@ #!/usr/bin/env python3 -"""Renders an image using raytracing.""" +"""Renders an image using raytracing. +The geometry and shading live in `raytracer`; what is left here is the scene +this lab builds (a checkerboard plane and a row of spheres) and the pixel loop +that shoots one primary ray per pixel and reflects it up to `DEPTH_MAX` times. +""" + +from __future__ import annotations + +import logging import time +from typing import TYPE_CHECKING import matplotlib.pyplot as plt import numpy as np +from raytracer import ( + SURFACE_OFFSET, + Lighting, + SceneObject, + make_plane, + make_sphere, + normalize, + trace_ray, +) + +if TYPE_CHECKING: + from pathlib import Path + +logger = logging.getLogger(__name__) + +# Default light and material parameters. +AMBIENT = 0.05 +DIFFUSE_C = 1.0 +SPECULAR_C = 1.0 +SPECULAR_K = 50 + +DEPTH_MAX = 5 # Maximum number of light reflections. + +# Progress is logged every this many pixel columns. +_PROGRESS_EVERY = 10 -def ray_trace( - num_spheres, environment, image_width=400, image_height=300, output_file="fig.png" -): - IMAGE_WIDTH = image_width - IMAGE_HEIGHT = image_height - - def normalize(vector): - """Normalize a vector. - - Parameters - ---------- - vector (numpy.ndarray): The input vector to be normalized. - - Returns - ------- - numpy.ndarray: The normalized vector. - - """ - vector /= np.linalg.norm(vector) - return vector - - def intersect_plane(ray_origin, ray_direction, plane_point, plane_normal): - """Calculate the intersection of a ray with a plane. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): A normalized 3D vector representing the - direction of the ray. - plane_point (numpy.ndarray): A 3D point representing a point on the plane. - plane_normal (numpy.ndarray): A normalized 3D vector representing - the normal of the plane. - - Returns - ------- - float: The distance from the origin ray_origin to the intersection - point with the plane. - Returns +inf if there is no intersection or if the intersection is - behind the origin. - - """ - denom = np.dot(ray_direction, plane_normal) - if np.abs(denom) < 1e-6: - return np.inf - d = np.dot(plane_point - ray_origin, plane_normal) / denom - if d < 0: - return np.inf - return d - - def intersect_sphere(ray_origin, ray_direction, sphere_center, sphere_radius): - """Calculate the intersection of a ray with a sphere. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): A normalized 3D vector representing the - direction of the ray. - sphere_center (numpy.ndarray): A 3D point representing - the center of the sphere. - sphere_radius (float): The radius of the sphere. - - Returns - ------- - float: The distance from the origin ray_origin to the intersection - point with the sphere. - Returns +inf if there is no intersection or if the intersection is - behind the origin. - - """ - a = np.dot(ray_direction, ray_direction) - origin_to_center = ray_origin - sphere_center - b = 2 * np.dot(ray_direction, origin_to_center) - radius_squared = sphere_radius * sphere_radius - c = np.dot(origin_to_center, origin_to_center) - radius_squared - disc = b * b - 4 * a * c - return calculate_sphere_intersection(a, b, c, disc) - - def calculate_sphere_intersection(a, b, c, disc): - """Calculate the - intersection distance of a ray with a sphere using the quadratic formula. - - Parameters - ---------- - a (float): Coefficient of t^2 in the quadratic equation. - b (float): Coefficient of t in the quadratic equation. - c (float): Constant term in the quadratic equation. - disc (float): Discriminant of the quadratic equation. - - Returns - ------- - float: - The distance from the origin to the intersection point with the sphere. - Returns +inf if there is no intersection - or if the intersection is behind the origin. - - """ - if disc > 0: - distance_squared = np.sqrt(disc) - # q is used to find the roots of the quadratic equation - if b < 0: - q = (-b - distance_squared) / 2.0 - else: - q = (-b + distance_squared) / 2.0 - t0 = q / a - t1 = c / q - t0, t1 = min(t0, t1), max(t0, t1) - if t1 >= 0: - return t1 if t0 < 0 else t0 - return np.inf - - def intersect(ray_origin, ray_direction, object_): - """Calculate the intersection of a ray with an object. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): A normalized 3D vector representing the - direction of the ray. - obj (dict): A dictionary representing the object with keys - 'type', 'position', 'normal' (for planes), and 'radius' (for spheres). - - Returns - ------- - float: The distance from the origin ray_origin to the intersection - point with the object. - Returns +inf if there is no intersection or if the intersection is - behind the origin. - - """ - if object_["type"] == "plane": - return intersect_plane( - ray_origin, ray_direction, object_["position"], object_["normal"] - ) - # object_['type'] == 'sphere': - return intersect_sphere( - ray_origin, ray_direction, object_["position"], object_["radius"] - ) - - def get_normal(object_, intersection_point): - """Calculate the normal at the intersection point on the object. - - Parameters - ---------- - obj (dict): A dictionary representing the object with keys - 'type' and 'position'. - intersection_point (numpy.ndarray): A 3D point representing the - intersection point on the object. - - Returns - ------- - numpy.ndarray: The normal vector at the intersection point. - - """ - if object_["type"] == "sphere": - normal = normalize(intersection_point - object_["position"]) - elif object_["type"] == "plane": - normal = object_["normal"] - else: - raise ValueError(f"Unknown object type: {object_['type']}") - return normal - - def get_color(object_, intersection_point): - """Get the color of the object at the intersection point. - - Parameters - ---------- - object_ (dict): A dictionary representing the object with a key 'color'. - intersection_point (numpy.ndarray): A 3D point representing the - intersection point on the object. - - Returns - ------- - numpy.ndarray: The color of the object at the intersection point. - - """ - color = object_["color"] - if not hasattr(color, "__len__"): - color = color(intersection_point) - return color - - def trace_ray(ray_origin, ray_direction): - """Trace a ray and find the color at the intersection point. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): - A normalized 3D vector representing the direction of the ray. - - Returns - ------- - tuple: A tuple containing the object, - intersection point, normal at the intersection, - and the color at the intersection point. - Returns None if there is no intersection. - - """ - t, obj_idx = find_intersection(ray_origin, ray_direction) - if t == np.inf: - return None - object_, intersection_point = get_intersection_details( - ray_origin, ray_direction, t, obj_idx - ) - normal, color = ( - get_normal(object_, intersection_point), - get_color(object_, intersection_point), - ) - if is_shadowed(intersection_point, normal, obj_idx): - return None - return compute_color(object_, intersection_point, normal, color, ray_origin) - - def find_intersection(ray_origin, ray_direction): - """Find the intersection of a ray with the objects in the scene. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): - A normalized 3D vector representing the direction of the ray. - - Returns - ------- - tuple: A tuple containing the distance to the intersection point - and the index of the intersected object. - - """ - t = np.inf - obj_idx = -1 - for index, object_ in enumerate(scene): - t_obj = intersect(ray_origin, ray_direction, object_) - if t_obj < t: - t, obj_idx = t_obj, index - return t, obj_idx - - def get_intersection_details(ray_origin, ray_direction, t, obj_idx): - """Get the details of the intersection point on the object. - - Parameters - ---------- - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - ray_direction (numpy.ndarray): - A normalized 3D vector representing the direction of the ray. - t (float): The distance to the intersection point. - obj_idx (int): The index of the intersected object in the scene. - - Returns - ------- - tuple: A tuple containing the intersected object - and the intersection point. - - """ - object_ = scene[obj_idx] - intersection_point = ray_origin + ray_direction * t - return object_, intersection_point - - def is_shadowed(intersection_point, normal, obj_idx): - """Determine if the intersection point is in shadow. - - Parameters - ---------- - intersection_point (numpy.ndarray): - A 3D point representing the intersection point on the object. - normal (numpy.ndarray): The normal vector at the intersection point. - obj_idx (int): The index of the intersected object in the scene. - - Returns - ------- - bool: True if the intersection point is in shadow, False otherwise. - - """ - to_light = normalize(L - intersection_point) - shadow_intersections = [ - intersect(intersection_point + normal * 0.0001, to_light, obj_sh) - for k, obj_sh in enumerate(scene) - if k != obj_idx - ] - return shadow_intersections and min(shadow_intersections) < np.inf - - def compute_color(object_, intersection_point, normal, color, ray_origin): - """Compute the color at the intersection point using shading techniques. - - Parameters - ---------- - object_ (dict): A dictionary representing the intersected object. - intersection_point (numpy.ndarray): - A 3D point representing the intersection point on the object. - normal (numpy.ndarray): The normal vector at the intersection point. - color (numpy.ndarray): The base color of the object. - ray_origin (numpy.ndarray): A 3D point representing the origin of the ray. - - Returns - ------- - tuple: - A tuple containing the intersected object, intersection point, normal, - and the computed color. - - """ - to_light = normalize(L - intersection_point) - to_origin = normalize(ray_origin - intersection_point) - color_ray = AMBIENT - diffuse_intensity = object_.get("diffuse_c", DIFFUSE_C) * max( - np.dot(normal, to_light), 0 - ) - color_ray += diffuse_intensity * color - half_vector = normalize(to_light + to_origin) - specular_intensity = ( - object_.get("specular_c", SPECULAR_C) - * max(np.dot(normal, half_vector), 0) ** SPECULAR_K - ) - color_ray += specular_intensity * color_light - return object_, intersection_point, normal, color_ray - - def add_sphere(position, radius, color): - """Create a dictionary representing a sphere object. - - Parameters - ---------- - position (list or numpy.ndarray): - A 3D point representing the position of the sphere. - radius (float): The radius of the sphere. - color (list or numpy.ndarray): The color of the sphere. - - Returns - ------- - dict: A dictionary representing the sphere object. - - """ - return { - "type": "sphere", - "position": np.array(position), - "radius": np.array(radius), - "color": np.array(color), - "reflection": 0.5, - } - - def add_plane(position, normal): - """Create a dictionary representing a plane object. - - Parameters - ---------- - position (list or numpy.ndarray): - A 3D point representing a point on the plane. - normal (list or numpy.ndarray): - A normalized 3D vector representing the normal of the plane. - - Returns - ------- - dict: A dictionary representing the plane object. - - """ - return { - "type": "plane", - "position": np.array(position), - "normal": np.array(normal), - "color": lambda M: ( - color_plane0 - if (int(M[0] * 2) % 2) == (int(M[2] * 2) % 2) - else color_plane1 - ), - "diffuse_c": 0.75, - "specular_c": 0.5, - "reflection": 0.25, - } - - scene = [] - - # List of objects. +def build_scene(num_spheres: int) -> list[SceneObject]: + """Lay out a checkerboard plane and `num_spheres` spheres in a receding row.""" color_plane0 = 1.0 * np.ones(3) color_plane1 = 0.0 * np.ones(3) - scene.append(add_plane([0.0, -0.5, 0.0], [0.0, 1.0, 0.0])) + scene = [make_plane([0.0, -0.5, 0.0], [0.0, 1.0, 0.0], color_plane0, color_plane1)] + base_radius = 1 / np.sqrt(num_spheres) # Im więcej kul, tym mniejsze base_distance = 4.5 / num_spheres @@ -390,67 +60,87 @@ def ray_trace( color = np.array([i / num_spheres, (num_spheres - i) / num_spheres, 0.5]) # Dodanie kuli do sceny - scene.append(add_sphere([x, y, z], base_radius, color)) + scene.append(make_sphere([x, y, z], base_radius, color)) - # Light position and color. - L = np.array([5.0, 5.0, -10.0]) - color_light = np.ones(3) + return scene - # Default light and material parameters. - AMBIENT = 0.05 - DIFFUSE_C = 1.0 - SPECULAR_C = 1.0 - SPECULAR_K = 50 - DEPTH_MAX = 5 # Maximum number of light reflections. - col = np.zeros(3) # Current color. +def ray_trace( + num_spheres: int, + _environment: str, + image_width: int = 400, + image_height: int = 300, + output_file: str | Path = "fig.png", +) -> None: + """Render `num_spheres` spheres over a checkerboard and write `output_file`. + + `_environment` is accepted for symmetry with the other two algorithms in + `main.py` and ignored: this renderer has no environment map. + """ + scene = build_scene(num_spheres) + lighting = Lighting( + position=np.array([5.0, 5.0, -10.0]), + color=np.ones(3), + ambient=AMBIENT, + diffuse_c=DIFFUSE_C, + specular_c=SPECULAR_C, + specular_k=SPECULAR_K, + ) + + color = np.zeros(3) # Current color. camera_origin = np.array([0.0, 0.35, -1.0]) # Camera. - Q = np.array([0.0, 0.0, 0.0]) # Camera pointing to. - img = np.zeros((IMAGE_HEIGHT, IMAGE_WIDTH, 3)) + camera_target = np.array([0.0, 0.0, 0.0]) # Camera pointing to. + img = np.zeros((image_height, image_width, 3)) - r = float(IMAGE_WIDTH) / IMAGE_HEIGHT + aspect = float(image_width) / image_height # Screen coordinates: x0, y0, x1, y1. - S = (-1.0, -1.0 / r + 0.25, 1.0, 1.0 / r + 0.25) + screen = (-1.0, -1.0 / aspect + 0.25, 1.0, 1.0 / aspect + 0.25) - renderTime = time.time() + render_time = time.time() reflections = 0 rays = 0 - initialRays = 0 + primary_rays = 0 # Loop through all pixels. - for i, x in enumerate(np.linspace(S[0], S[2], IMAGE_WIDTH)): - if i % 10 == 0: - print(round(i / float(IMAGE_WIDTH) * 100, 2), "%") - for j, y in enumerate(np.linspace(S[1], S[3], IMAGE_HEIGHT)): - col[:] = 0 - Q[:2] = (x, y) - D = normalize(Q - camera_origin) - DEPTH = 0 - rayO, rayD = camera_origin, D - REFLECTION = 1.0 - initialRays += 1 + for i, x in enumerate(np.linspace(screen[0], screen[2], image_width)): + if i % _PROGRESS_EVERY == 0: + logger.info("%s %%", round(i / float(image_width) * 100, 2)) + for j, y in enumerate(np.linspace(screen[1], screen[3], image_height)): + color[:] = 0 + camera_target[:2] = (x, y) + ray_origin = camera_origin + ray_direction = normalize(camera_target - camera_origin) + depth = 0 + reflection = 1.0 + primary_rays += 1 # Loop through initial and secondary rays. - while DEPTH < DEPTH_MAX: - traced = trace_ray(rayO, rayD) + while depth < DEPTH_MAX: + traced = trace_ray(scene, lighting, ray_origin, ray_direction) rays += 1 if not traced: break reflections += 1 - obj, M, N, col_ray = traced + obj, hit_point, hit_normal, color_ray = traced # Reflection: create a new ray. - rayO, rayD = M + N * 0.0001, normalize(rayD - 2 * np.dot(rayD, N) * N) - DEPTH += 1 - col += REFLECTION * col_ray - REFLECTION *= obj.get("reflection", 1.0) - img[IMAGE_HEIGHT - j - 1, i, :] = np.clip(col, 0, 1) - renderTime = time.time() - renderTime + ray_origin = hit_point + hit_normal * SURFACE_OFFSET + ray_direction = normalize( + ray_direction - 2 * np.dot(ray_direction, hit_normal) * hit_normal + ) + depth += 1 + color += reflection * color_ray + reflection *= obj.get("reflection", 1.0) + img[image_height - j - 1, i, :] = np.clip(color, 0, 1) + render_time = time.time() - render_time plt.imsave(output_file, img) - print( - f"Image saved as {output_file}\n" - f"resolution: {IMAGE_WIDTH}x{IMAGE_HEIGHT}\n" - f"render time: {round(renderTime, 2)} s\n" - f"reflections: {reflections}\n" - f"rays (initial): {initialRays}\n" - f"rays (secondary): {rays - initialRays}\n" - f"rays (total): {rays}" + logger.info( + "Image saved as %s\nresolution: %sx%s\nrender time: %s s\nreflections: %s\n" + "rays (initial): %s\nrays (secondary): %s\nrays (total): %s", + output_file, + image_width, + image_height, + round(render_time, 2), + reflections, + primary_rays, + rays - primary_rays, + rays, ) diff --git a/Programming/TRAK/scenes/__init__.py b/Programming/TRAK/scenes/__init__.py index e69de29b..eaacfcd0 100644 --- a/Programming/TRAK/scenes/__init__.py +++ b/Programming/TRAK/scenes/__init__.py @@ -0,0 +1 @@ +"""Scene definitions, each exposing a `setup_scene(width, height, environment)`.""" diff --git a/Programming/TRAK/scenes/cornell_box.py b/Programming/TRAK/scenes/cornell_box.py index b4254025..b4156cf2 100644 --- a/Programming/TRAK/scenes/cornell_box.py +++ b/Programming/TRAK/scenes/cornell_box.py @@ -1,14 +1,35 @@ -from sightpy import * +"""The Cornell box: a closed room lit by a ceiling panel. + +The reference scene for global illumination -- the colour bleeding from the +red and green walls onto the white ones is the whole point of it. +""" + +from sightpy import ( + Cuboid, + Diffuse, + Emissive, + Plane, + Refractive, + Scene, + Sphere, + rgb, + vec3, +) -# define materials to use -def setup_scene(width=400, height=300, environment=None): +def setup_scene( + width: int = 400, height: int = 300, environment: str | None = None +) -> Scene: + """Build the Cornell box at `width` x `height`. - Sc = Scene(ambient_color=rgb(0.00, 0.00, 0.00)) + `environment` is accepted so every scene has one interface, and ignored: + the box is closed, so no background is ever visible. + """ + del environment - angle = -0 + scene = Scene(ambient_color=rgb(0.00, 0.00, 0.00)) - Sc.add_Camera( + scene.add_Camera( screen_width=width, screen_height=height, look_from=vec3(278, 278, 800), @@ -23,11 +44,10 @@ def setup_scene(width=400, height=300, environment=None): red_diffuse = Diffuse(diff_color=rgb(0.65, 0.05, 0.05)) white_diffuse = Diffuse(diff_color=rgb(0.73, 0.73, 0.73)) emissive_white = Emissive(color=rgb(15.0, 15.0, 15.0)) - emissive_blue = Emissive(color=rgb(2.0, 2.0, 3.5)) blue_glass = Refractive(n=vec3(1.5 + 0.05e-8j, 1.5 + 0.02e-8j, 1.5 + 0.0j)) # this is the light - Sc.add( + scene.add( Plane( material=emissive_white, center=vec3(213 + 130 / 2, 554, -227.0 - 105 / 2), @@ -39,7 +59,7 @@ def setup_scene(width=400, height=300, environment=None): importance_sampled=True, ) - Sc.add( + scene.add( Plane( material=white_diffuse, center=vec3(555 / 2, 555 / 2, -555.0), @@ -50,7 +70,7 @@ def setup_scene(width=400, height=300, environment=None): ) ) - Sc.add( + scene.add( Plane( material=green_diffuse, center=vec3(-0.0, 555 / 2, -555 / 2), @@ -61,7 +81,7 @@ def setup_scene(width=400, height=300, environment=None): ) ) - Sc.add( + scene.add( Plane( material=red_diffuse, center=vec3(555.0, 555 / 2, -555 / 2), @@ -72,7 +92,7 @@ def setup_scene(width=400, height=300, environment=None): ) ) - Sc.add( + scene.add( Plane( material=white_diffuse, center=vec3(555 / 2, 555, -555 / 2), @@ -83,7 +103,7 @@ def setup_scene(width=400, height=300, environment=None): ) ) - Sc.add( + scene.add( Plane( material=white_diffuse, center=vec3(555 / 2, 0.0, -555 / 2), @@ -103,9 +123,9 @@ def setup_scene(width=400, height=300, environment=None): shadow=False, ) cb.rotate(θ=15, u=vec3(0, 1, 0)) - Sc.add(cb) + scene.add(cb) - Sc.add( + scene.add( Sphere( material=blue_glass, center=vec3(370.5, 165 / 2, -65 - 185 / 2), @@ -116,4 +136,4 @@ def setup_scene(width=400, height=300, environment=None): importance_sampled=True, ) - return Sc + return scene diff --git a/Programming/TRAK/scenes/soap_bubble.py b/Programming/TRAK/scenes/soap_bubble.py index 9a2386e7..731b3a00 100644 --- a/Programming/TRAK/scenes/soap_bubble.py +++ b/Programming/TRAK/scenes/soap_bubble.py @@ -1,14 +1,22 @@ -from sightpy import * +"""One soap bubble in front of a blurred backdrop. + +The colours come from thin-film interference, not from a texture: the shell +is a few hundred nanometres thick, so different wavelengths cancel at +different angles. +""" + +import numpy as np +from sightpy import Scene, Sphere, ThinFilmInterference, rgb, vec3 -# define materials to use -def setup_scene(width=400, height=300, environment="lake.png"): - # Set Scene - - Sc = Scene(ambient_color=rgb(0.01, 0.01, 0.01)) +def setup_scene( + width: int = 400, height: int = 300, environment: str = "lake.png" +) -> Scene: + """Build the soap-bubble scene at `width` x `height` over `environment`.""" + scene = Scene(ambient_color=rgb(0.01, 0.01, 0.01)) angle = -np.pi * 0.5 - Sc.add_Camera( + scene.add_Camera( screen_height=height, screen_width=width, look_from=vec3(4.0 * np.sin(angle), 0.00, 4.0 * np.cos(angle)), @@ -16,7 +24,7 @@ def setup_scene(width=400, height=300, environment="lake.png"): ) soap_bubble = ThinFilmInterference(thickness=330, noise=60.0) - Sc.add( + scene.add( Sphere( material=soap_bubble, center=vec3(1.0, 0.0, 1.5), @@ -26,6 +34,6 @@ def setup_scene(width=400, height=300, environment="lake.png"): ) ) - Sc.add_Background(environment, blur=10.0) + scene.add_Background(environment, blur=10.0) - return Sc + return scene diff --git a/Programming/TRAK/scenes/three_spheres.py b/Programming/TRAK/scenes/three_spheres.py index 0df887ce..8ef69708 100644 --- a/Programming/TRAK/scenes/three_spheres.py +++ b/Programming/TRAK/scenes/three_spheres.py @@ -1,9 +1,17 @@ -from sightpy import * +"""Three coloured glass spheres on a checkered floor, over a photo backdrop. + +Each sphere filters a different channel, so the floor seen through them shows +what the refractive material does to each colour separately. +""" + +import numpy as np +from sightpy import Glossy, Plane, Refractive, Scene, Sphere, image, rgb, vec3 -# define materials to use -def setup_scene(width=400, height=300, environment="miramar.jpeg"): - +def setup_scene( + width: int = 400, height: int = 300, environment: str = "miramar.jpeg" +) -> Scene: + """Build the three-sphere scene at `width` x `height` over `environment`.""" blue_glass = Refractive( n=vec3(1.5 + 4e-8j, 1.5 + 4e-8j, 1.5 + 0.0j) ) # n = index of refraction @@ -20,19 +28,19 @@ def setup_scene(width=400, height=300, environment="miramar.jpeg"): # Set Scene - Sc = Scene(ambient_color=rgb(0.05, 0.05, 0.05)) + scene = Scene(ambient_color=rgb(0.05, 0.05, 0.05)) angle = np.pi / 2 * 0.3 - Sc.add_Camera( + scene.add_Camera( look_from=vec3(2.5 * np.sin(angle), 0.25, 2.5 * np.cos(angle) - 1.5), look_at=vec3(0.0, 0.25, -1.5), screen_width=width, screen_height=height, ) - Sc.add_DirectionalLight(Ldir=vec3(0.52, 0.45, -0.5), color=rgb(0.15, 0.15, 0.15)) + scene.add_DirectionalLight(Ldir=vec3(0.52, 0.45, -0.5), color=rgb(0.15, 0.15, 0.15)) - Sc.add( + scene.add( Sphere( material=blue_glass, center=vec3(-1.2, 0.0, -1.5), @@ -41,7 +49,7 @@ def setup_scene(width=400, height=300, environment="miramar.jpeg"): max_ray_depth=3, ) ) - Sc.add( + scene.add( Sphere( material=green_glass, center=vec3(0.0, 0.0, -1.5), @@ -50,7 +58,7 @@ def setup_scene(width=400, height=300, environment="miramar.jpeg"): max_ray_depth=3, ) ) - Sc.add( + scene.add( Sphere( material=red_glass, center=vec3(1.2, 0.0, -1.5), @@ -60,7 +68,7 @@ def setup_scene(width=400, height=300, environment="miramar.jpeg"): ) ) - Sc.add( + scene.add( Plane( material=floor, center=vec3(0, -0.5, -3.0), @@ -73,6 +81,6 @@ def setup_scene(width=400, height=300, environment="miramar.jpeg"): ) # see sightpy/backgrounds - Sc.add_Background(environment) + scene.add_Background(environment) - return Sc + return scene diff --git a/Programming/TRAK/scenes/transparent_cuboid.py b/Programming/TRAK/scenes/transparent_cuboid.py index c9e32e18..106829ea 100644 --- a/Programming/TRAK/scenes/transparent_cuboid.py +++ b/Programming/TRAK/scenes/transparent_cuboid.py @@ -1,8 +1,16 @@ -from sightpy import * +"""A rotated block of green glass standing on a checkered floor. + +The floor is deliberately close and low-repeat, so the refraction through the +block is easy to compare against the undistorted tiles beside it. +""" + +from sightpy import Cuboid, Glossy, Plane, Refractive, Scene, image, rgb, vec3 -# define materials to use -def setup_scene(width=400, height=300, environment="stormydays.png"): +def setup_scene( + width: int = 400, height: int = 300, environment: str = "stormydays.png" +) -> Scene: + """Build the glass-cuboid scene at `width` x `height` over `environment`.""" floor = Glossy( diff_color=image("checkered_floor.png", repeat=2.0), roughness=0.2, @@ -12,17 +20,17 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): ) # n = index of refraction green_glass = Refractive(n=vec3(1.5 + 4e-8j, 1.5 + 0.0j, 1.5 + 4e-8j)) - Sc = Scene() - Sc.add_Camera( + scene = Scene() + scene.add_Camera( look_from=vec3(0.0, 0.25, 1.0), look_at=vec3(0.0, 0.25, -3.0), screen_width=width, screen_height=height, ) - Sc.add_DirectionalLight(Ldir=vec3(0.0, 0.5, 0.5), color=rgb(0.5, 0.5, 0.5)) + scene.add_DirectionalLight(Ldir=vec3(0.0, 0.5, 0.5), color=rgb(0.5, 0.5, 0.5)) - Sc.add( + scene.add( Plane( material=floor, center=vec3(0, -0.5, -3.0), @@ -44,9 +52,9 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): max_ray_depth=5, ) cb.rotate(θ=30, u=vec3(0, 1, 0)) - Sc.add(cb) + scene.add(cb) # see sightpy/backgrounds - Sc.add_Background(environment) + scene.add_Background(environment) - return Sc + return scene diff --git a/Programming/TRAK/scenes/two_spheres.py b/Programming/TRAK/scenes/two_spheres.py index 5f08b654..8d1c8c84 100644 --- a/Programming/TRAK/scenes/two_spheres.py +++ b/Programming/TRAK/scenes/two_spheres.py @@ -1,8 +1,17 @@ -from sightpy import * +"""Two metal spheres -- gold and bluish -- on a checkered floor. + +Both use the Glossy material, so this is the scene where the complex index of +refraction and the roughness parameter are easiest to read off the render. +""" + +import numpy as np +from sightpy import Glossy, Plane, Scene, Sphere, image, rgb, vec3 -# define materials to use -def setup_scene(width=400, height=300, environment="stormydays.png"): +def setup_scene( + width: int = 400, height: int = 300, environment: str = "stormydays.png" +) -> Scene: + """Build the two-sphere scene at `width` x `height` over `environment`.""" gold_metal = Glossy( diff_color=rgb(1.0, 0.572, 0.184), n=vec3(0.15 + 3.58j, 0.4 + 2.37j, 1.54 + 1.91j), @@ -27,19 +36,19 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): ) # Set Scene - Sc = Scene(ambient_color=rgb(0.05, 0.05, 0.05)) + scene = Scene(ambient_color=rgb(0.05, 0.05, 0.05)) angle = -np.pi / 2 * 0.3 - Sc.add_Camera( + scene.add_Camera( look_from=vec3(2.5 * np.sin(angle), 0.25, 2.5 * np.cos(angle) - 1.5), look_at=vec3(0.0, 0.25, -3.0), screen_width=width, screen_height=height, ) - Sc.add_DirectionalLight(Ldir=vec3(0.52, 0.45, -0.5), color=rgb(0.15, 0.15, 0.15)) + scene.add_DirectionalLight(Ldir=vec3(0.52, 0.45, -0.5), color=rgb(0.15, 0.15, 0.15)) - Sc.add( + scene.add( Sphere( material=gold_metal, center=vec3(-0.75, 0.1, -3.0), @@ -47,7 +56,7 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): max_ray_depth=3, ) ) - Sc.add( + scene.add( Sphere( material=bluish_metal, center=vec3(1.25, 0.1, -3.0), @@ -56,7 +65,7 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): ) ) - Sc.add( + scene.add( Plane( material=floor, center=vec3(0, -0.5, -3.0), @@ -69,6 +78,6 @@ def setup_scene(width=400, height=300, environment="stormydays.png"): ) # see sightpy/backgrounds - Sc.add_Background(environment) + scene.add_Background(environment) - return Sc + return scene diff --git a/Programming/TRAK/utils.py b/Programming/TRAK/utils.py old mode 100644 new mode 100755 index e59a979e..8464e263 --- a/Programming/TRAK/utils.py +++ b/Programming/TRAK/utils.py @@ -1,17 +1,25 @@ +#!/usr/bin/env python3 +"""Configuration and CLI-argument helpers shared by the renderer entry point.""" + from configparser import ConfigParser -def load_config(config_path): +def load_config(config_path: str) -> ConfigParser: + """Read `config_path` and return the parsed configuration. + + A missing file is not an error: ``ConfigParser.read`` ignores it and the + caller falls back to argparse defaults. + """ config = ConfigParser() config.read(config_path) return config -def parse_resolution(resolution): +def parse_resolution(resolution: str) -> tuple[int, int]: + """Split a ``WIDTHxHEIGHT`` string into its two integer components.""" try: width, height = map(int, resolution.lower().split("x")) - return width, height - except ValueError: - raise ValueError( - "Resolution must be in the format WIDTHxHEIGHT, e.g., 1920x1080." - ) + except ValueError as err: + msg = "Resolution must be in the format WIDTHxHEIGHT, e.g., 1920x1080." + raise ValueError(msg) from err + return width, height diff --git a/TODO-refactor-remaining.md b/TODO-refactor-remaining.md index a6d1ed04..7d076508 100644 --- a/TODO-refactor-remaining.md +++ b/TODO-refactor-remaining.md @@ -32,14 +32,13 @@ own docs and the rest of this fleet carry. A test asserts by raising. |---|---| | Baseline, before any of this | 9,281 | | After the format + safe-fix pass and excluding `TRAK/sightpy` | 3,686 | -| **Now** | **1,183** | +| **Now** | **651** | Courses at zero: `WDWR`, `SPD`, `MOM`, `ECRYPT`, `ECRYPT_PROJECT`, `ECOTE`, -`NLP`, `EARIN`, `PORR`. What is left, by course: +`NLP`, `EARIN`, `PORR`, `TRAK`. What is left, by course: | Course | ruff | Can it be verified by running? | |---|---|---| -| `TRAK` | 531 | **Yes** -- `run` ray-traces cornell_box in 37s. `sightpy/` is vendored and excluded; this is kuhy's own code | | `twm_4` | 228 | No: `jupyter` on PATH is a dead pipx shim. `pip install jupyter` in a venv would fix that | | `ERSMS-project` | 156 | No: docker-compose course | | `PSD` | 152 | No: `run.sh` predates the harness and installs packages with sudo | @@ -53,10 +52,10 @@ deliverable. This exact bug has already been found twice (TRAK's ## mypy has NOT been run against any of this -`pyproject.toml` sets `strict = true` and `CLAUDE.md` says the repo runs mypy, -but mypy is not installed on this machine (pip is externally managed; use a -venv) and no commit so far has type-checked. **"Zero findings" in the commit -messages means zero RUFF findings.** Measured after installing it: +`pyproject.toml` sets `strict = true` and `CLAUDE.md` says the repo runs mypy. +mypy IS installed now (`~/.local/bin/mypy`) -- the earlier note here saying it +was not is out of date. **"Zero findings" in a commit message still means zero +RUFF findings unless it says otherwise**; only `TRAK` has been type-checked. | Course | mypy --strict errors | |---|---| @@ -64,7 +63,7 @@ messages means zero RUFF findings.** Measured after installing it: | `ECOTE` | 14 | | `NLP` | 29 | | `EARIN/lab1` alone | 46 | -| `TRAK` | 1 | +| `TRAK` | **0** -- cleared 2026-08-29 | | `PSD/zin1` | 5 | Two structural notes for whoever runs it next: @@ -74,15 +73,24 @@ Two structural notes for whoever runs it next: - Most errors are `Returning Any` from untyped third-party calls (openai, pandas) and genuinely-wrong unions this pass introduced by hand, e.g. `str | bool` returns in `gpt_chunks.reformat`. They are real. +- Point it at the course's own venv or every import is unresolvable: + `mypy --strict --python-executable Programming//.venv/bin/python `. + Without that, TRAK reported 4 phantom `matplotlib` import errors on top of + the real ones. +- `pyproject.toml` now carries a `[[tool.mypy.overrides]]` block with + `ignore_missing_imports` for `Imath`, `OpenEXR`, `OpenGL` and `matplotlib`. + Those four ship neither stubs nor `py.typed`; without it mypy cannot import + them at all and reports nothing else about the file. It is not a suppression + of our own code -- do not extend it to anything we wrote. ## Still open -### 1. The 250-line cap: 51 files +### 1. The 250-line cap: 50 files | Type | Count | Note | |---|---|---| | `.txt` | 27 | Blocked -- see `TODO-file-length-250.md` | -| `.py` | 11 | Several GREW during the lint pass, because docstrings and annotations are lines: `NLP/gpt_chunks.py` 911 -> 973, `EARIN/lab2/main.py` 581 -> 647, `PORR/.../linear_algebra_utils.py` 422 -> 553. Split them AFTER the lint work, not during | +| `.py` | 10 | Several GREW during the lint pass, because docstrings and annotations are lines: `NLP/gpt_chunks.py` 911 -> 973, `EARIN/lab2/main.py` 581 -> 647, `PORR/.../linear_algebra_utils.py` 422 -> 553. Split them AFTER the lint work, not during. TRAK is already clear: its three over-cap files were split in the same pass that linted them, each verified by a byte-identical render | | `.cpp` | 9 | | | `.c` | 4 | | | `.java` | 2 | | @@ -120,6 +128,15 @@ And run mypy per directory afterwards: it found 12 latent TypeErrors in EARIN that ruff could not see -- call sites still passing positionally to arguments the pass had made keyword-only. +Prefer a check the program itself settles. TRAK's `ray_tracing0` path has no +RNG at all, so `main.py --algorithm ray_tracing0 --output ../build/x.png` +renders byte-identically every time: refactoring 109 findings out of +`rendering.py` and splitting it in two was verified by one md5 comparison. The +cornell_box path is stochastic and cannot be used that way -- its ray count +drifts by ~20 between runs on its own. Where an RNG is unavoidable, do the +restructure and the `NPY002` swap as SEPARATE steps: the first must be +bit-identical under an external seed, only the second may move the pixels. + ### 6. Disk `Programming/PORR/code/` holds 3.8 GB of generated .npz caches (poli3 alone is diff --git a/pyproject.toml b/pyproject.toml index 72d5af28..32c7d28e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -73,6 +73,13 @@ exclude = [ "^Programming/TRAK/sightpy/", ] +# Third-party packages that ship neither stubs nor a py.typed marker. This is +# not a suppression of our own code: without it mypy cannot import them at all +# and reports nothing else about the file that uses them. +[[tool.mypy.overrides]] +module = ["Imath.*", "OpenEXR.*", "OpenGL.*", "matplotlib.*"] +ignore_missing_imports = true + [tool.coverage.run] branch = true omit = [