diff --git a/.gitignore b/.gitignore index f9949030..126fa691 100644 --- a/.gitignore +++ b/.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/ \ No newline at end of file +#.idea/ + +*.exr \ No newline at end of file diff --git a/code/mapa_srodowiska/README.md b/code/mapa_srodowiska/README.md new file mode 100644 index 00000000..a0efbc03 --- /dev/null +++ b/code/mapa_srodowiska/README.md @@ -0,0 +1,2 @@ +Download file in (!) EXR (!) format +https://polyhaven.com/a/lilienstein \ No newline at end of file diff --git a/code/mapa_srodowiska/lilienstein_1k.exr b/code/mapa_srodowiska/lilienstein_1k.exr new file mode 100644 index 00000000..ebcdb200 Binary files /dev/null and b/code/mapa_srodowiska/lilienstein_1k.exr differ diff --git a/code/mapa_srodowiska/main.py b/code/mapa_srodowiska/main.py new file mode 100644 index 00000000..311581f8 --- /dev/null +++ b/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) diff --git a/code/photonmapping/main.py b/code/photonmapping/main.py new file mode 100644 index 00000000..9aec82d6 --- /dev/null +++ b/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()