Kuhyx (#1)

* feat: load exr map

* feat: scale and flip the hdr map

* feat: added camera

* feat: add simple slow but working photonmapping

.gitignore

@@ -161,4 +161,6 @@ cython_debug/
 #  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
 #  and can be added to the global gitignore or merged into this file.  For a more nuclear
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
-#.idea/
+#.idea/
+
+*.exr

code/mapa_srodowiska/README.md

@@ -0,0 +1,2 @@
+Download file in (!) EXR (!) format
+https://polyhaven.com/a/lilienstein

code/mapa_srodowiska/main.py

@@ -0,0 +1,162 @@
+import OpenEXR
+import Imath
+import os
+import OpenGL.GL as gl
+import OpenGL.GLUT as glut
+import OpenGL.GLU as glu
+import math
+
+# 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
+
+def load_hdr_environment_map(filepath):
+    """
+    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.
+    """
+    if not os.path.exists(filepath):
+        raise FileNotFoundError(f"File not found: {filepath}")
+    
+    # Check file permissions
+    if not os.access(filepath, os.R_OK):
+        raise PermissionError(f"File is not readable: {filepath}")
+    
+    # Open the EXR file
+    try:
+        exr_file = OpenEXR.InputFile(filepath)
+    except Exception as e:
+        raise OSError(f"Unable to open '{filepath}' for read: {str(e)}")
+    
+    # Get the image dimensions
+    header = exr_file.header()
+    dw = header['dataWindow']
+    width = dw.max.x - dw.min.x + 1
+    height = dw.max.y - dw.min.y + 1
+
+    # Define the channel names (R, G, B)
+    channels = ['R', 'G', 'B']
+    
+    # Read the channel data
+    channel_data = {
+        channel: exr_file.channel(channel, Imath.PixelType(Imath.PixelType.FLOAT))
+        for channel in channels
+    }
+    
+    # Combine channel data into a single bytes object
+    hdr_image = bytearray(width * height * 3 * 4)  # 3 channels, 4 bytes per float
+    for i, channel in enumerate(channels):
+        channel_buffer = channel_data[channel]
+        for j in range(height):
+            for k in range(width):
+                index = (j * width + k) * 3 * 4 + i * 4
+                hdr_image[index:index + 4] = channel_buffer[(j * width + k) * 4:(j * width + k + 1) * 4]
+
+    # Flip the image vertically
+    flipped_hdr_image = bytearray(width * height * 3 * 4)
+    row_size = width * 3 * 4
+    for j in range(height):
+        src_index = j * row_size
+        dst_index = (height - 1 - j) * row_size
+        flipped_hdr_image[dst_index:dst_index + row_size] = hdr_image[src_index:src_index + row_size]
+
+    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)
+    
+    # Enable texture mapping
+    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)
+    
+    # Load the texture
+    gl.glTexImage2D(gl.GL_TEXTURE_2D, 0, gl.GL_RGB32F, width, height, 0, 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()
+    
+    # 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)
+    
+    # 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)
+    
+    # Disable texture mapping
+    gl.glDisable(gl.GL_TEXTURE_2D)
+    
+    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_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 main(filepath):    
+    try:
+        width, height, hdr_map = load_hdr_environment_map(filepath)
+        print("HDR Map Loaded. Dimensions:", width, "x", height)
+        
+        # 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 display callback
+        glut.glutDisplayFunc(lambda: display_hdr_image(width, height, hdr_map))
+        
+        # Set mouse callbacks
+        glut.glutMotionFunc(mouse_motion)
+        glut.glutMouseFunc(mouse_button)
+        
+        # Start the GLUT main loop
+        glut.glutMainLoop()
+    except Exception as e:
+        print(e)
+
+# Example usage
+if __name__ == "__main__":
+    filepath = "lilienstein_4k.exr"
+    main(filepath)

code/photonmapping/main.py

@@ -0,0 +1,221 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# 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
+        else:
+            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
+
+# Scene setup
+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]
+
+# Light source
+light_position = Vector3(-5, 5, -5)
+light_power = np.array([1, 1, 1]) * 1000  # Intense white light
+
+# Photon map
+photon_map = []
+
+# Parameters
+num_photons = 10000  # Number of photons to emit
+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):
+    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)
+    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 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):
+    closest_t = np.inf
+    hit_object = None
+    hit_info = None
+    # Find the nearest intersection
+    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)
+    else:
+        return np.array([0, 0, 0])  # Background color
+
+def compute_direct_light(point, normal):
+    # 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
+    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])
+    else:
+        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
+    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 * num_photons)
+
+# Main execution
+if __name__ == '__main__':
+    print("Emitting photons...")
+    emit_photons()
+    
+    print("Rendering image...")
+    width = 200
+    height = 100
+    image = render_image(width, height)
+    
+    # Display the image
+    plt.imshow(image)
+    plt.axis('off')
+    plt.show()