From 732ffb38f1da2355c670acdcc6ba450d1069e3a3 Mon Sep 17 00:00:00 2001 From: Krzysztof Rudnicki Date: Fri, 27 Dec 2024 11:36:24 +0100 Subject: [PATCH] feat: add primitive raytracing --- .python-version | 1 + code/main.py | 215 ++++++++++++++++++++++++++++-------------- code/requirements.txt | 3 +- 3 files changed, 146 insertions(+), 73 deletions(-) create mode 100644 .python-version diff --git a/.python-version b/.python-version new file mode 100644 index 00000000..2d4715b6 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.11.11 diff --git a/code/main.py b/code/main.py index 3bff121f..3a45de39 100644 --- a/code/main.py +++ b/code/main.py @@ -1,87 +1,158 @@ -from OpenGL.GL import * -from OpenGL.GLUT import * -from OpenGL.GL.shaders import compileProgram, compileShader import numpy as np +import matplotlib.pyplot as plt -# Vertex shader -VERTEX_SHADER = """ -#version 330 -in vec3 position; -in vec3 color; -out vec3 vertexColor; -void main() { - gl_Position = vec4(position, 1.0); - vertexColor = color; -} -""" +w = 400 +h = 300 -# Fragment shader -FRAGMENT_SHADER = """ -#version 330 -in vec3 vertexColor; -out vec4 fragColor; -void main() { - fragColor = vec4(vertexColor, 1.0); -} -""" +def normalize(x): + x /= np.linalg.norm(x) + return x -# Define vertices and colors -vertices = np.array([ - # Positions # Colors - 0.0, 0.5, 0.0, 1.0, 0.0, 0.0, # Top (red) - -0.5, -0.5, 0.0, 0.0, 1.0, 0.0, # Bottom-left (green) - 0.5, -0.5, 0.0, 0.0, 0.0, 1.0 # Bottom-right (blue) -], dtype=np.float32) +def intersect_plane(O, D, P, N): + # Return the distance from O to the intersection of the ray (O, D) with the + # plane (P, N), or +inf if there is no intersection. + # O and P are 3D points, D and N (normal) are normalized vectors. + denom = np.dot(D, N) + if np.abs(denom) < 1e-6: + return np.inf + d = np.dot(P - O, N) / denom + if d < 0: + return np.inf + return d -# Initialize OpenGL -def init(): - global shader, VAO +def intersect_sphere(O, D, S, R): + # Return the distance from O to the intersection of the ray (O, D) with the + # sphere (S, R), or +inf if there is no intersection. + # O and S are 3D points, D (direction) is a normalized vector, R is a scalar. + a = np.dot(D, D) + OS = O - S + b = 2 * np.dot(D, OS) + c = np.dot(OS, OS) - R * R + disc = b * b - 4 * a * c + if disc > 0: + distSqrt = np.sqrt(disc) + q = (-b - distSqrt) / 2.0 if b < 0 else (-b + distSqrt) / 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 - # Compile shaders - shader = compileProgram( - compileShader(VERTEX_SHADER, GL_VERTEX_SHADER), - compileShader(FRAGMENT_SHADER, GL_FRAGMENT_SHADER) - ) +def intersect(O, D, obj): + if obj['type'] == 'plane': + return intersect_plane(O, D, obj['position'], obj['normal']) + elif obj['type'] == 'sphere': + return intersect_sphere(O, D, obj['position'], obj['radius']) - # Generate VAO and VBO - VAO = glGenVertexArrays(1) - VBO = glGenBuffers(1) +def get_normal(obj, M): + # Find normal. + if obj['type'] == 'sphere': + N = normalize(M - obj['position']) + elif obj['type'] == 'plane': + N = obj['normal'] + return N + +def get_color(obj, M): + color = obj['color'] + if not hasattr(color, '__len__'): + color = color(M) + return color - glBindVertexArray(VAO) +def trace_ray(rayO, rayD): + # Find first point of intersection with the scene. + t = np.inf + for i, obj in enumerate(scene): + t_obj = intersect(rayO, rayD, obj) + if t_obj < t: + t, obj_idx = t_obj, i + # Return None if the ray does not intersect any object. + if t == np.inf: + return + # Find the object. + obj = scene[obj_idx] + # Find the point of intersection on the object. + M = rayO + rayD * t + # Find properties of the object. + N = get_normal(obj, M) + color = get_color(obj, M) + toL = normalize(L - M) + toO = normalize(O - M) + # Shadow: find if the point is shadowed or not. + l = [intersect(M + N * .0001, toL, obj_sh) + for k, obj_sh in enumerate(scene) if k != obj_idx] + if l and min(l) < np.inf: + return + # Start computing the color. + col_ray = ambient + # Lambert shading (diffuse). + col_ray += obj.get('diffuse_c', diffuse_c) * max(np.dot(N, toL), 0) * color + # Blinn-Phong shading (specular). + col_ray += obj.get('specular_c', specular_c) * max(np.dot(N, normalize(toL + toO)), 0) ** specular_k * color_light + return obj, M, N, col_ray - glBindBuffer(GL_ARRAY_BUFFER, VBO) - glBufferData(GL_ARRAY_BUFFER, vertices.nbytes, vertices, GL_STATIC_DRAW) +def add_sphere(position, radius, color): + return dict(type='sphere', position=np.array(position), + radius=np.array(radius), color=np.array(color), reflection=.5) + +def add_plane(position, normal): + return dict(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=.75, specular_c=.5, reflection=.25) + +# List of objects. +color_plane0 = 1. * np.ones(3) +color_plane1 = 0. * np.ones(3) +scene = [add_sphere([.75, .1, 1.], .6, [1., 0., 0.]), + add_sphere([-.75, .1, 2.25], .6, [0., 1., 0.]), + add_sphere([-2.75, .1, 3.5], .6, [0., 0., 1.]), + add_plane([0., -.5, 0.], [0., 1., 0.]), + ] - # Position attribute - glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 6 * vertices.itemsize, None) - glEnableVertexAttribArray(0) +# Light position and color. +L = np.array([5., 5., -10.]) +color_light = np.ones(3) - # Color attribute - glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 6 * vertices.itemsize, ctypes.c_void_p(3 * vertices.itemsize)) - glEnableVertexAttribArray(1) +# Default light and material parameters. +ambient = .05 +diffuse_c = 1. +specular_c = 1. +specular_k = 50 - glBindBuffer(GL_ARRAY_BUFFER, 0) - glBindVertexArray(0) +depth_max = 5 # Maximum number of light reflections. +col = np.zeros(3) # Current color. +O = np.array([0., 0.35, -1.]) # Camera. +Q = np.array([0., 0., 0.]) # Camera pointing to. +img = np.zeros((h, w, 3)) -# Render function -def display(): - glClear(GL_COLOR_BUFFER_BIT) - glUseProgram(shader) - glBindVertexArray(VAO) - glDrawArrays(GL_TRIANGLES, 0, 3) - glBindVertexArray(0) - glUseProgram(0) - glutSwapBuffers() +r = float(w) / h +# Screen coordinates: x0, y0, x1, y1. +S = (-1., -1. / r + .25, 1., 1. / r + .25) -# Main function -def main(): - glutInit() - glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB) - glutCreateWindow(b"PyOpenGL Triangle") - glutInitWindowSize(800, 600) - glutDisplayFunc(display) - init() - glutMainLoop() +# Loop through all pixels. +for i, x in enumerate(np.linspace(S[0], S[2], w)): + if i % 10 == 0: + print(i / float(w) * 100, "%") + for j, y in enumerate(np.linspace(S[1], S[3], h)): + col[:] = 0 + Q[:2] = (x, y) + D = normalize(Q - O) + depth = 0 + rayO, rayD = O, D + reflection = 1. + # Loop through initial and secondary rays. + while depth < depth_max: + traced = trace_ray(rayO, rayD) + if not traced: + break + obj, M, N, col_ray = traced + # Reflection: create a new ray. + rayO, rayD = M + N * .0001, normalize(rayD - 2 * np.dot(rayD, N) * N) + depth += 1 + col += reflection * col_ray + reflection *= obj.get('reflection', 1.) + img[h - j - 1, i, :] = np.clip(col, 0, 1) -if __name__ == "__main__": - main() +plt.imsave('fig.png', img) \ No newline at end of file diff --git a/code/requirements.txt b/code/requirements.txt index 82f9807d..ac121f2d 100644 --- a/code/requirements.txt +++ b/code/requirements.txt @@ -1,4 +1,5 @@ glfw numpy pyrr -PyOpenGL \ No newline at end of file +PyOpenGL +matplotlib \ No newline at end of file