diff --git a/code/__init__.py b/code/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/code/output_ray_traced.png b/code/output_ray_traced.png new file mode 100644 index 00000000..7da47ec3 Binary files /dev/null and b/code/output_ray_traced.png differ diff --git a/code/requirements.txt b/code/requirements.txt index 4f6b1783..3085b712 100644 --- a/code/requirements.txt +++ b/code/requirements.txt @@ -6,4 +6,6 @@ matplotlib flake8 black autopep8 -flake8-max-function-length \ No newline at end of file +flake8-max-function-length +pillow +progressbar \ No newline at end of file diff --git a/code/scenes/EXAMPLE1.png b/code/scenes/EXAMPLE1.png new file mode 100644 index 00000000..c4215e8d Binary files /dev/null and b/code/scenes/EXAMPLE1.png differ diff --git a/code/sightpy/__init__.py b/code/sightpy/__init__.py new file mode 100644 index 00000000..19c250df --- /dev/null +++ b/code/sightpy/__init__.py @@ -0,0 +1,17 @@ +from .utils.constants import * +from .utils.vector3 import * +from .utils.colour_functions import * +from .utils.image_functions import * + +from .ray import * +from .scene import * +from .geometry import * +from .lights import * +from .materials import * +from .textures.texture import * +from .animation import * + + + + + diff --git a/code/sightpy/animation.py b/code/sightpy/animation.py new file mode 100644 index 00000000..30411a9b --- /dev/null +++ b/code/sightpy/animation.py @@ -0,0 +1,54 @@ +from .scene import Scene +import numpy as np +from pathlib import Path + + +def create_animation(scene,samples_per_pixel, fps, start_time, final_time, update_scene, name): + + """ + this function render a list of frames and saves them in ./frames folder. You can make an animation the using ffmpeg running + from the command prompt: + """ + #ffmpeg -r 60 -f image2 -s 854x480 -i your_image_%d.png -vcodec libx264 -crf 1 -pix_fmt yuv420p your_video.mp4 + #fps #resoluion #crf = quality (less is better) + + + number_of_frames = int(fps*(final_time - start_time)) + dt = (final_time - start_time)/number_of_frames + t = start_time + + try: + Path("./frames").mkdir() + + except FileExistsError: + pass + + + for i in range(0,number_of_frames): + update_scene(scene, t) + img = scene.render(samples_per_pixel) + t += dt + img.save("frames/" + name + "_" + str(i) + ".png") + + + + +def create_animation_using_opencv(scene, samples_per_pixel , fps, start_time, final_time, update_scene, name): + + import cv2 + number_of_frames = int(fps*(final_time - start_time)) + dt = (final_time - start_time)/number_of_frames + t = start_time + + + videodims = (scene.camera.screen_width, scene.camera.screen_height) + fourcc = cv2.VideoWriter_fourcc('M', 'J', 'P', 'G') + video = cv2.VideoWriter(name,fourcc, fps,videodims) + + for i in range(0,number_of_frames): + update_scene(scene, t) + frame = scene.render(samples_per_pixel) + video.write(cv2.cvtColor(np.array(frame), cv2.COLOR_RGB2BGR)) + t += dt + + video.release() diff --git a/code/sightpy/backgrounds/__init__.py b/code/sightpy/backgrounds/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/code/sightpy/backgrounds/lake.png b/code/sightpy/backgrounds/lake.png new file mode 100644 index 00000000..f7541a2c Binary files /dev/null and b/code/sightpy/backgrounds/lake.png differ diff --git a/code/sightpy/backgrounds/lightmaps/lake.png b/code/sightpy/backgrounds/lightmaps/lake.png new file mode 100644 index 00000000..349a96a6 Binary files /dev/null and b/code/sightpy/backgrounds/lightmaps/lake.png differ diff --git a/code/sightpy/backgrounds/miramar.jpeg b/code/sightpy/backgrounds/miramar.jpeg new file mode 100644 index 00000000..1a720c20 Binary files /dev/null and b/code/sightpy/backgrounds/miramar.jpeg differ diff --git a/code/sightpy/backgrounds/panorama.py b/code/sightpy/backgrounds/panorama.py new file mode 100644 index 00000000..3eb9867f --- /dev/null +++ b/code/sightpy/backgrounds/panorama.py @@ -0,0 +1,18 @@ +from ..geometry import Sphere_Collider, Primitive +from ..materials import Material +from ..utils.vector3 import vec3 +from ..utils.constants import SKYBOX_DISTANCE +from ..utils.image_functions import load_image, load_image_as_linear_sRGB +from .util.blur_background import blur_skybox +from .skybox import SkyBox_Material + +class Panorama(Primitive): + def __init__(self, panorama, center = vec3(0.,0.,0.), light_intensity = 0.0, blur = 0.0): + super().__init__(center, SkyBox_Material(panorama, light_intensity, blur), shadow = False) + l = SKYBOX_DISTANCE + self.light_intensity = light_intensity + self.collider_list += [Sphere_Collider(assigned_primitive = self, center = center , radius = SKYBOX_DISTANCE)] + + + def get_uv(self, hit): + return hit.collider.get_uv(hit) \ No newline at end of file diff --git a/code/sightpy/backgrounds/skybox.py b/code/sightpy/backgrounds/skybox.py new file mode 100644 index 00000000..93802d37 --- /dev/null +++ b/code/sightpy/backgrounds/skybox.py @@ -0,0 +1,57 @@ +from ..geometry import Cuboid_Collider, Primitive +from ..materials import Material +from ..utils.vector3 import vec3 +from ..utils.constants import SKYBOX_DISTANCE +from ..utils.image_functions import load_image, load_image_as_linear_sRGB +from .util.blur_background import blur_skybox + +class SkyBox(Primitive): + def __init__(self, cubemap, center = vec3(0.,0.,0.), light_intensity = 0.0, blur = 0.0): + super().__init__(center, SkyBox_Material(cubemap, light_intensity, blur), shadow = False) + l = SKYBOX_DISTANCE + self.light_intensity = light_intensity + #BOTTOM + self.collider_list += [Cuboid_Collider(assigned_primitive = self, center = center, width = 2*l, height =2*l ,length =2*l )] + + + def get_uv(self, hit): + u,v = hit.collider.get_uv(hit) + u,v = u/4,v/3 + return u,v + + +class SkyBox_Material(Material): + def __init__(self, cubemap, light_intensity, blur): + self.texture = load_image_as_linear_sRGB("sightpy/backgrounds/" + cubemap) + + if light_intensity != 0.0: + self.lightmap = load_image("sightpy/backgrounds/lightmaps/" + cubemap) + + if blur != 0.0: + self.blur_image = blur_skybox(load_image("sightpy/backgrounds/" + cubemap), blur, cubemap) + + self.blur = blur + self.light_intensity = light_intensity + self.repeat = 1.0 + + def get_texture_color(self, hit, ray): + u,v = hit.get_uv() + + if (self.blur != 0.0) : + im = self.blur_image[-((v * self.blur_image.shape[0]*self.repeat ).astype(int)% self.blur_image.shape[0]) , (u * self.blur_image.shape[1]*self.repeat).astype(int) % self.blur_image.shape[1] ].T + else: + im = self.texture[-((v * self.texture.shape[0]*self.repeat ).astype(int)% self.texture.shape[0]) , (u * self.texture.shape[1]*self.repeat).astype(int) % self.texture.shape[1] ].T + + if (ray.depth != 0) and (self.light_intensity != 0.0): + ls = self.lightmap[-((v * self.texture.shape[0]*self.repeat ).astype(int)% self.texture.shape[0]) , (u * self.texture.shape[1]*self.repeat).astype(int) % self.texture.shape[1] ].T + color = vec3(im[0] + self.light_intensity * ls[0], im[1] + self.light_intensity * ls[1], im[2] + self.light_intensity * ls[2]) + + else: + color = vec3(im[0] , im[1] , im[2] ) + return color + + + def get_color(self, scene, ray, hit): + hit.point = (ray.origin + ray.dir * hit.distance) + return hit.material.get_texture_color(hit,ray) + diff --git a/code/sightpy/backgrounds/stormydays.png b/code/sightpy/backgrounds/stormydays.png new file mode 100644 index 00000000..0e76d1a1 Binary files /dev/null and b/code/sightpy/backgrounds/stormydays.png differ diff --git a/code/sightpy/backgrounds/util/__init__.py b/code/sightpy/backgrounds/util/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/code/sightpy/backgrounds/util/blur_background.py b/code/sightpy/backgrounds/util/blur_background.py new file mode 100644 index 00000000..02f6ca5d --- /dev/null +++ b/code/sightpy/backgrounds/util/blur_background.py @@ -0,0 +1,158 @@ +from PIL import Image, ImageFilter +import numpy as np +from ...utils.colour_functions import sRGB_to_sRGB_linear + + + +def to_image(arr): + img = [Image.fromarray((255 * arr[:,:,i]).astype(np.uint8), "L") for i in range(0,3)] + return Image.merge("RGB", img) + +def to_array(img): + return np.asarray(img)/256. + + + + + +def blur_skybox(img_array, blur, cubemap): + + print("blurring " + cubemap) + + N = int(img_array.shape[0]/3) + + left = img_array[ 1*N:2*N, 0*N:1*N] + front = img_array[N:2*N, N:2*N] + right = img_array[ N:2*N, 2*N:3*N] + back = img_array[N:2*N, 3*N:4*N] + top = img_array[0:1*N, 1*N:2*N] + bottom = img_array[2*N:3*N, 1*N:2*N] + + + cubelist = [[None,top,None, None], + [left,front,right, back], + [None,bottom,None, None]] + + + + back_blur = np.zeros((N*3, N*3,3)) + + + back_blur[ 1*N:2*N, 0*N:1*N] = cubelist[1][0-2] #left + back_blur[N:2*N, N:2*N] = cubelist[1][1-2] #front + back_blur[ N:2*N, 2*N:3*N] = cubelist[1][2-2] # right + back_blur[2*N:3*N, 1*N:2*N] = np.rot90(cubelist[2][1] , k=2) #bottom + back_blur[0:1*N, 1*N:2*N] = np.rot90(cubelist[0][1], k=2) #top + + + back_blur = to_image(back_blur) + back_blur = (back_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + back_blur = to_array(back_blur) + + back_blur = back_blur[N:2*N, N:2*N] + + + top_blur = np.zeros((N*3, N*3,3)) + top_blur[ 1*N:2*N, 0*N:1*N] = np.rot90(cubelist[1][0], k=-1) #left + top_blur[N:2*N, N:2*N] = cubelist[1-1][1] #front + top_blur[ N:2*N, 2*N:3*N] = np.rot90(cubelist[1][2], k=1) # right + top_blur[2*N:3*N, 1*N:2*N] = cubelist[1][1] #bottom + top_blur[0:1*N, 1*N:2*N] = np.rot90(cubelist[1][3], k=2) #top + top_blur = to_image(top_blur) + top_blur = (top_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + top_blur = to_array(top_blur) + + top_blur = top_blur[N:2*N, N:2*N] + + + + + bottom_blur = np.zeros((N*3, N*3,3)) + + + bottom_blur[ 1*N:2*N, 0*N:1*N] = np.rot90(cubelist[1][0], k=1) #left + bottom_blur[N:2*N, N:2*N] = cubelist[1+1][1] #front + bottom_blur[ N:2*N, 2*N:3*N] = np.rot90(cubelist[1][2], k=-1) # right + bottom_blur[2*N:3*N, 1*N:2*N] = np.rot90(cubelist[1][3], k=2) #bottom + bottom_blur[0:1*N, 1*N:2*N] = cubelist[1][1] #top + + bottom_blur = to_image(bottom_blur) + bottom_blur = (bottom_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + bottom_blur = to_array(bottom_blur) + + bottom_blur = bottom_blur[N:2*N, N:2*N] + + + + right_blur = np.zeros((N*3, N*3,3)) + + right_blur[ 1*N:2*N, 0*N:1*N] = cubelist[1][0+1] #left + right_blur[N:2*N, N:2*N] = cubelist[1][1+1] #front + right_blur[ N:2*N, 2*N:3*N] = cubelist[1][2+1] # right + right_blur[2*N:3*N, 1*N:2*N] = np.rot90(cubelist[2][1] ) #bottom + right_blur[0:1*N, 1*N:2*N] = np.rot90(cubelist[0][1], k=-1) #top + + right_blur = to_image(right_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + + right_blur = to_array(right_blur) + + right_blur = right_blur[N:2*N, N:2*N] + + + + front_blur = np.zeros((N*3, N*3,3)) + + + front_blur [ 1*N:2*N, 0*N:1*N] = cubelist[1][0] #left + front_blur [N:2*N, N:2*N] = cubelist[1][1] #front + front_blur [ N:2*N, 2*N:3*N] = cubelist[1][2] # right + front_blur [2*N:3*N, 1*N:2*N] = cubelist[2][1] #bottom + front_blur [0:1*N, 1*N:2*N] = cubelist[0][1] #top + + front_blur = to_image(front_blur) + + front_blur = (front_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + front_blur = to_array(front_blur) + + front_blur = front_blur[N:2*N, N:2*N] + + + + + left_blur = np.zeros((N*3, N*3,3)) + + left_blur[ 1*N:2*N, 0*N:1*N] = cubelist[1][0-1] #left + left_blur[N:2*N, N:2*N] = cubelist[1][1-1] #front + left_blur[ N:2*N, 2*N:3*N] = cubelist[1][2-1] + left_blur[2*N:3*N, 1*N:2*N] = np.rot90(cubelist[2][1], k=-1 ) #bottom + left_blur[0:1*N, 1*N:2*N] = np.rot90(cubelist[0][1], k=1) #top + + + left_blur = to_image(left_blur).filter(ImageFilter.GaussianBlur(radius=blur)) + + + left_blur = to_array(left_blur) + + left_blur = left_blur[N:2*N, N:2*N] + + + + + + + skybox_blurred = np.zeros((N*3, N*4,3)) + + skybox_blurred[ 1*N:2*N, 0*N:1*N] = left_blur + skybox_blurred[N:2*N, N:2*N] = front_blur + skybox_blurred[ N:2*N, 2*N:3*N] = right_blur + skybox_blurred[N:2*N, 3*N:4*N] = back_blur + skybox_blurred[0:1*N, 1*N:2*N] = top_blur + skybox_blurred[2*N:3*N, 1*N:2*N] = bottom_blur + + return sRGB_to_sRGB_linear(skybox_blurred) + diff --git a/code/sightpy/camera.py b/code/sightpy/camera.py new file mode 100644 index 00000000..0f29133a --- /dev/null +++ b/code/sightpy/camera.py @@ -0,0 +1,52 @@ +from .utils.vector3 import vec3, rgb +from .utils.random import random_in_unit_disk +import numpy as np +from .ray import Ray + +class Camera(): + def __init__(self, look_from, look_at, screen_width = 400 ,screen_height = 300, field_of_view = 90., aperture = 0., focal_distance = 1.): + self.screen_width = screen_width + self.screen_height = screen_height + self.aspect_ratio = float(screen_width) / screen_height + + self.look_from = look_from + self.look_at = look_at + self.camera_width = np.tan(field_of_view * np.pi/180 /2.)*2. + self.camera_height = self.camera_width/self.aspect_ratio + + #camera reference basis in world coordinates + self.cameraFwd = (look_at - look_from).normalize() + self.cameraRight = (self.cameraFwd.cross(vec3(0.,1.,0.))).normalize() + self.cameraUp = self.cameraRight.cross(self.cameraFwd) + + + + #if you use a lens_radius >= 0.0 make sure that samples_per_pixel is a large number. Otherwise you'll get a lot of noise + self.lens_radius = aperture / 2. + self.focal_distance = focal_distance + + + + # Pixels coordinates in camera basis: + self.x = np.linspace(-self.camera_width/2., self.camera_width/2., self.screen_width) + self.y = np.linspace(self.camera_height/2., -self.camera_height/2., self.screen_height) + + # we are going to cast a total of screen_width * screen_height * samples_per_pixel rays + # xx,yy store the origin of each ray in a 3d array where the first and second dimension are the x,y coordinates of each pixel + # and the third dimension is the sample index of each pixel + xx,yy = np.meshgrid(self.x,self.y) + self.x = xx.flatten() + self.y = yy.flatten() + + def get_ray(self,n): # n = index of refraction of scene main medium (for air n = 1.) + + # in each pixel, take a random position to avoid aliasing. + x = self.x + (np.random.rand(len(self.x )) - 0.5)*self.camera_width /(self.screen_width) + y = self.y + (np.random.rand(len(self.y )) - 0.5)*self.camera_height /(self.screen_height) + + # set ray direction in world space: + rx, ry = random_in_unit_disk(x.shape[0]) + ray_origin = self.look_from + self.cameraRight *rx* self.lens_radius + self.cameraUp *ry* self.lens_radius + ray_dir = (self.look_from + self.cameraUp*y*self.focal_distance + self.cameraRight*x*self.focal_distance + self.cameraFwd*self.focal_distance - ray_origin ).normalize() + return Ray(origin=ray_origin, dir=ray_dir, depth=0, n=n, reflections = 0, transmissions = 0, diffuse_reflections = 0) + diff --git a/code/sightpy/geometry/__init__.py b/code/sightpy/geometry/__init__.py new file mode 100644 index 00000000..cdeaa1bd --- /dev/null +++ b/code/sightpy/geometry/__init__.py @@ -0,0 +1,13 @@ +from .primitive import * +from .collider import * +from .sphere import * +from .plane import * +from .triangle import * +from .triangle_mesh import * +from .cuboid import * + + + + + + diff --git a/code/sightpy/geometry/collider.py b/code/sightpy/geometry/collider.py new file mode 100644 index 00000000..97a500fa --- /dev/null +++ b/code/sightpy/geometry/collider.py @@ -0,0 +1,17 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from abc import abstractmethod + +class Collider: + def __init__(self,assigned_primitive, center): + self.assigned_primitive = assigned_primitive + self.center = center + + @abstractmethod + def intersect(self, O, D): + pass + + @abstractmethod + def get_Normal(self, hit): + pass \ No newline at end of file diff --git a/code/sightpy/geometry/cuboid.py b/code/sightpy/geometry/cuboid.py new file mode 100644 index 00000000..086b5231 --- /dev/null +++ b/code/sightpy/geometry/cuboid.py @@ -0,0 +1,145 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from ..geometry import Primitive, Collider + +class Cuboid(Primitive): + def __init__(self,center, material, width,height, length,max_ray_depth = 5, shadow = True): + super().__init__(center, material, max_ray_depth, shadow = shadow) + self.width = width + self.height = height + self.length = length + self.bounded_sphere_radius = np.sqrt((self.width/2)**2 + (self.height/2)**2 + (self.length/2)**2) + + self.collider_list += [Cuboid_Collider(assigned_primitive = self, center = center, width = width, height =height ,length =length )] + + + def get_uv(self, hit): + u,v = hit.collider.get_uv(hit) + u,v = u/4,v/3 + return u,v + + +""" + This was the old approach, but remplaced by a Box collider that is more efficient + #we model a cuboid as six planes + + + #BOTTOM #BOTTOM + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(0.0,-h, 0.0), u_axis = vec3(1.0, 0.0, 0.0), v_axis = vec3(0.0, 0.0, 1.0), w = w, h = l, uv_shift = (1,0))] + #TOP #TOP + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(0.0,h, 0.0), u_axis = vec3(1.0, 0.0, 0.0), v_axis = vec3(0.0, 0.0, -1.0), w = w, h = l, uv_shift= (1,2))] + #RIGHT #RIGHT + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(w,0.0, 0.0), u_axis = vec3(0.0, 0.0, -1.0), v_axis = vec3(0.0, 1.0, 0.0), w = l, h = h, uv_shift= (2,1))] + #LEFT #LEFT + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(-w,0.0, 0.0), u_axis = vec3(0.0, 0.0, 1.0), v_axis = vec3(0.0, 1.0, 0.0), w = l, h = h, uv_shift= (0,1))] + #FRONT #FRONT + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(0,0, l), u_axis = vec3(1.0, 0.0, 0.0), v_axis = vec3(0.0, 1.0, 0.0), w = w, h = h, uv_shift= (1,1))] + #BACK #BACK + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center + vec3(0,0, -l), u_axis = vec3(-1.0, 0.0, 0.0), v_axis = vec3(0.0, 1.0, 0.0), w = w, h = h, uv_shift= (3,1))] +""" + + +class Cuboid_Collider(Collider): + def __init__(self, width, height,length,**kwargs): + super().__init__(**kwargs) + + self.lb = self.center - vec3(width/2, height/2, length/2) + self.rt = self.center + vec3(width/2, height/2, length/2) + + self.lb_local_basis = self.lb + self.rt_local_basis = self.rt + + self.width = width + self.height = height + self.length = length + + # basis vectors + self.ax_w = vec3(1.,0.,0.) + self.ax_h = vec3(0.,1.,0.) + self.ax_l = vec3(0.,0.,1.) + + self.inverse_basis_matrix = np.array([[self.ax_w.x, self.ax_h.x, self.ax_l.x], + [self.ax_w.y, self.ax_h.y, self.ax_l.y], + [self.ax_w.z, self.ax_h.z, self.ax_l.z]]) + + self.basis_matrix = self.inverse_basis_matrix.T + + + def rotate(self,M, center): + self.ax_w = self.ax_w.matmul(M) + self.ax_h = self.ax_h.matmul(M) + self.ax_l = self.ax_l.matmul(M) + + self.inverse_basis_matrix = np.array([[self.ax_w.x, self.ax_h.x, self.ax_l.x], + [self.ax_w.y, self.ax_h.y, self.ax_l.y], + [self.ax_w.z, self.ax_h.z, self.ax_l.z]]) + + self.basis_matrix = self.inverse_basis_matrix.T + + self.lb = center + (self.lb-center).matmul(M) + self.rt = center + (self.rt-center).matmul(M) + + self.lb_local_basis = self.lb.matmul(self.basis_matrix) + self.rt_local_basis = self.rt.matmul(self.basis_matrix) + + def intersect(self, O, D): + + + O_local_basis = O.matmul(self.basis_matrix) + D_local_basis = D.matmul(self.basis_matrix) + + dirfrac = 1.0 / D_local_basis + + # lb is the corner of AABB with minimal coordinates - left bottom, rt is maximal corner + t1 = (self.lb_local_basis.x - O_local_basis.x)*dirfrac.x; + t2 = (self.rt_local_basis.x - O_local_basis.x)*dirfrac.x; + t3 = (self.lb_local_basis.y - O_local_basis.y)*dirfrac.y; + t4 = (self.rt_local_basis.y - O_local_basis.y)*dirfrac.y; + t5 = (self.lb_local_basis.z - O_local_basis.z)*dirfrac.z; + t6 = (self.rt_local_basis.z - O_local_basis.z)*dirfrac.z; + + tmin = np.maximum(np.maximum(np.minimum(t1, t2), np.minimum(t3, t4)), np.minimum(t5, t6)) + tmax = np.minimum(np.minimum(np.maximum(t1, t2), np.maximum(t3, t4)), np.maximum(t5, t6)) + + # if tmax < 0, ray (line) is intersecting AABB, but the whole AABB is behind us + # if tmin > tmax, ray doesn't intersect AAB + mask1 = (tmax < 0) | (tmin > tmax) + + # if tmin < 0 then the ray origin is inside of the AABB and tmin is behind the start of the ray so tmax is the first intersection + mask2 = tmin < 0 + return np.select([mask1,mask2,True] , [FARAWAY , [tmax, np.tile(UPDOWN, tmin.shape)] , [tmin, np.tile(UPWARDS, tmin.shape)]]) + + + def get_Normal(self, hit): + + P = (hit.point-self.center).matmul(self.basis_matrix) + absP = vec3(1./self.width, 1./self.height, 1./self.length)*np.abs(P) + Pmax = np.maximum(np.maximum(absP.x, absP.y), absP.z) + P.x = np.where(Pmax == absP.x, np.sign(P.x), 0.) + P.y = np.where(Pmax == absP.y, np.sign(P.y), 0.) + P.z = np.where(Pmax == absP.z, np.sign(P.z), 0.) + + return P.matmul(self.inverse_basis_matrix) + + def get_uv(self, hit): + hit.N = self.get_Normal(hit) + M_C = hit.point - self.center + + BOTTOM = (hit.N == vec3(0.,-1.,0.)) + TOP = (hit.N == vec3(0., 1.,0.)) + RIGHT = (hit.N == vec3(1.,0.,0.)) + LEFT = (hit.N == vec3(-1.,0.,0.) ) + FRONT = (hit.N == vec3(0.,0.,1.)) + BACK = (hit.N == vec3(0.,0.,-1.)) + + #0.985 to avoid corners + u = np.select([BOTTOM , TOP, RIGHT, LEFT , FRONT , BACK], [((self.ax_w.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1), ((self.ax_w.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1), ((self.ax_l.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 2), (((self.ax_l*-1).dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 0), (((self.ax_w*-1).dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 3), (( self.ax_w.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1)]) + v = np.select([BOTTOM , TOP, RIGHT, LEFT , FRONT , BACK], [(((self.ax_l*-1).dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 0), ((self.ax_l.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 2), ((self.ax_h.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1), (((self.ax_h).dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1), (((self.ax_h).dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1), (( self.ax_h.dot(M_C)/self.width*2 *0.985 + 1 ) /2 + 1)]) + return u,v + + + + + + diff --git a/code/sightpy/geometry/plane.py b/code/sightpy/geometry/plane.py new file mode 100644 index 00000000..9af46b1c --- /dev/null +++ b/code/sightpy/geometry/plane.py @@ -0,0 +1,81 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from ..geometry import Primitive, Collider + + +class Plane(Primitive): + def __init__(self,center, material, width,height, u_axis, v_axis, max_ray_depth = 5, shadow = True): + super().__init__(center, material, max_ray_depth,shadow = shadow) + self.collider_list += [Plane_Collider(assigned_primitive = self, center = center, u_axis = u_axis, v_axis = v_axis, w= width/2, h=height/2)] + self.width = width + self.height = height + self.bounded_sphere_radius = np.sqrt((width/2)**2 + (height/2)**2) + + def get_uv(self, hit): + return hit.collider.get_uv(hit) + + +class Plane_Collider(Collider): + def __init__(self, u_axis, v_axis, w, h, uv_shift = (0.,0.),**kwargs): + super().__init__(**kwargs) + self.normal = u_axis.cross(v_axis).normalize() + + + self.w = w + self.h = h + self.u_axis = u_axis + self.v_axis = v_axis + self.uv_shift = uv_shift + self.inverse_basis_matrix = np.array([[self.u_axis.x, self.v_axis.x, self.normal.x], + [self.u_axis.y, self.v_axis.y, self.normal.y], + [self.u_axis.z, self.v_axis.z, self.normal.z]]) + self.basis_matrix = self.inverse_basis_matrix.T + + + + + def intersect(self, O, D): + N = self.normal + + NdotD = N.dot(D) + NdotD = np.where(NdotD == 0., NdotD + 0.0001, NdotD) #avoid zero division + + NdotC_O = N.dot(self.center - O) + d = D * NdotC_O / NdotD + M = O + d # intersection point + dis = d.length() + + M_C = M - self.center + + #plane basis coordinates + u = self.u_axis.dot(M_C) + v = self.v_axis.dot(M_C) + + + hit_inside = (np.abs(u) <= self.w) & (np.abs(v) <= self.h) & (NdotC_O * NdotD > 0) + hit_UPWARDS = (NdotD < 0) + hit_UPDOWN = np.logical_not(hit_UPWARDS) + + + pred1 = hit_inside & hit_UPWARDS + pred2 = hit_inside & hit_UPDOWN + pred3 = True + return np.select([pred1,pred2,pred3] , [[dis, np.tile(UPWARDS, dis.shape) ], [dis,np.tile(UPDOWN, dis.shape)], FARAWAY]) + + def rotate(self,M, center): + self.u_axis = self.u_axis.matmul(M) + self.v_axis = self.v_axis.matmul(M) + self.normal = self.normal.matmul(M) + self.center = center + (self.center-center).matmul(M) + + def get_uv(self, hit): + M_C = hit.point - self.center + u = ((self.u_axis.dot(M_C)/self.w + 1 ) /2 + self.uv_shift[0]) + v = ((self.v_axis.dot(M_C)/self.h + 1 ) /2 + self.uv_shift[1]) + return u,v + + + def get_Normal(self, hit): + return self.normal + diff --git a/code/sightpy/geometry/primitive.py b/code/sightpy/geometry/primitive.py new file mode 100644 index 00000000..9895a942 --- /dev/null +++ b/code/sightpy/geometry/primitive.py @@ -0,0 +1,30 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3 +import numpy as np + + +class Primitive: + def __init__(self, center, material, max_ray_depth = 1, shadow = True): + self.center = center + self.material = material + self.material.assigned_primitive = self + self.shadow = shadow + self.collider_list = [] + self.max_ray_depth = max_ray_depth + + def rotate(self, θ, u): + + u = u.normalize() + θ = θ/180 *np.pi + cosθ = np.cos(θ) + sinθ = np.sqrt(1-cosθ**2) * np.sign(θ) + + #rotation matrix along u axis + M = np.array([ + [cosθ + u.x*u.x * (1-cosθ), u.x*u.y*(1-cosθ) - u.z*sinθ, u.x*u.z*(1-cosθ) +u.y*sinθ], + [u.y*u.x*(1-cosθ) + u.z*sinθ, cosθ + u.y**2 * (1-cosθ), u.y*u.z*(1-cosθ) -u.x*sinθ], + [u.z*u.x*(1-cosθ) -u.y*sinθ, u.z*u.y*(1-cosθ) + u.x*sinθ, cosθ + u.z*u.z*(1-cosθ)] + ]) + for c in self.collider_list: + c.rotate(M, self.center) + diff --git a/code/sightpy/geometry/sphere.py b/code/sightpy/geometry/sphere.py new file mode 100644 index 00000000..1d748a02 --- /dev/null +++ b/code/sightpy/geometry/sphere.py @@ -0,0 +1,52 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from ..geometry import Primitive, Collider + +class Sphere(Primitive): + def __init__(self,center, material, radius, max_ray_depth = 5, shadow = True): + super().__init__(center, material, max_ray_depth, shadow = shadow) + self.collider_list += [Sphere_Collider(assigned_primitive = self, center = center, radius = radius)] + self.bounded_sphere_radius = radius + + def get_uv(self, hit): + return hit.collider.get_uv(hit) + + +class Sphere_Collider(Collider): + def __init__(self, radius, **kwargs): + super().__init__(**kwargs) + self.radius = radius + + + def intersect(self, O, D): + + b = 2 * D.dot(O - self.center) + c = self.center.square_length() + O.square_length() - 2 * self.center.dot(O) - (self.radius * self.radius) + disc = (b ** 2) - (4 * c) + sq = np.sqrt(np.maximum(0, disc)) + h0 = (-b - sq) / 2 + h1 = (-b + sq) / 2 + h = np.where((h0 > 0) & (h0 < h1), h0, h1) + pred = (disc > 0) & (h > 0) + M = (O + D * h) + NdotD = ((M - self.center) * (1. / self.radius) ).dot(D) + + pred1 = (disc > 0) & (h > 0) & (NdotD > 0) + pred2 = (disc > 0) & (h > 0) & (NdotD < 0) + pred3 = True + + #return an array with hit distance and the hit orientation + return np.select([pred1,pred2,pred3] , [[h, np.tile(UPDOWN, h.shape)], [h,np.tile(UPWARDS, h.shape)], FARAWAY]) + + def get_Normal(self, hit): + # M = intersection point + return (hit.point - self.center) * (1. / self.radius) + + def get_uv(self, hit): + M_C = (hit.point - self.center) / self.radius + phi = np.arctan2(M_C.z, M_C.x) + theta = np.arcsin(M_C.y) + u = (phi + np.pi) / (2*np.pi) + v = (theta + np.pi/2) / np.pi + return u,v \ No newline at end of file diff --git a/code/sightpy/geometry/surface.py b/code/sightpy/geometry/surface.py new file mode 100644 index 00000000..600bbb92 --- /dev/null +++ b/code/sightpy/geometry/surface.py @@ -0,0 +1,30 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3 +import numpy as np + + +class Surface: + def __init__(self, center, material, shadow = True): + self.center = center + self.material = material + self.material.assigned_surface = self + self.shadow = shadow + self.collider_list = [] + + + def rotate(self, θ, u): + + u = u.normalize() + θ = θ/180 *np.pi + cosθ = np.cos(θ) + sinθ = np.sqrt(1-cosθ**2) * np.sign(θ) + + #rotation matrix along u axis + M = np.array([ + [cosθ + u.x*u.x * (1-cosθ), u.x*u.y*(1-cosθ) - u.z*sinθ, u.x*u.z*(1-cosθ) +u.y*sinθ], + [u.y*u.x*(1-cosθ) + u.z*sinθ, cosθ + u.y**2 * (1-cosθ), u.y*u.z*(1-cosθ) -u.x*sinθ], + [u.z*u.x*(1-cosθ) -u.y*sinθ, u.z*u.y*(1-cosθ) + u.x*sinθ, cosθ + u.z*u.z*(1-cosθ)] + ]) + for c in self.collider_list: + c.rotate(M, self.center) + diff --git a/code/sightpy/geometry/triangle.py b/code/sightpy/geometry/triangle.py new file mode 100644 index 00000000..bee1cb87 --- /dev/null +++ b/code/sightpy/geometry/triangle.py @@ -0,0 +1,77 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from ..geometry.primitive import Primitive +from ..geometry.collider import Collider + + +class Triangle(Primitive): + def __init__(self,center, material, p1 , p2, p3, max_ray_depth,shadow = True): + super().__init__(center, material, max_ray_depth, shadow = shadow) + self.collider_list += [Triangle_Collider(assigned_primitive = self, p1 =p1, p2 = p2, p3 = p3)] + + def get_uv(self, M, collider): + return collider.get_uv(M) + + +class Triangle_Collider(Collider): + def __init__(self,assigned_surface, p1, p2, p3): + + + self.assigned_primitive = assigned_surface + self.p1 = p1 + self.p2 = p2 + self.p3 = p3 + self.normal = ((self.p2 - self.p1).cross( self.p3 - self.p1)).normalize() + + + self.centroid = (self.p1 + self.p2 + self.p3)/3 #one possible definition of center. Used for intersect(). + + + self.n31 = (self.p3 - self.p1).cross(self.normal) + self.n12 = (self.p1- self.p2).cross(self.normal) + self.n23 = (self.p2 - self.p3).cross(self.normal) + + + + def intersect(self, O, D): + N = self.normal + + NdotD = N.dot(D) + NdotD = np.where(NdotD == 0., NdotD + 0.0001, NdotD) #avoid zero division + + NdotC_O = N.dot(self.centroid - O) + d = D * NdotC_O / NdotD + M = O + d # intersection point + dis = d.length() + M_C = M - self.centroid + hit_inside = (self.n31.dot(M-self.p1) >= 0) & (self.n12.dot(M-self.p2) >= 0)& (self.n23.dot(M-self.p3) >= 0) & (NdotC_O * NdotD > 0) + hit_UPWARDS = (NdotD < 0) + hit_UPDOWN = np.logical_not(hit_UPWARDS) + + + pred1 = hit_inside & hit_UPWARDS + pred2 = hit_inside & hit_UPDOWN + pred3 = True + return np.select([pred1,pred2,pred3] , [[dis, np.tile(UPWARDS, dis.shape) ], [dis,np.tile(UPDOWN, dis.shape)], FARAWAY]) + + def rotate(self,M, center): + self.p1 = center + (self.p1 -center).matmul(M) + self.p2 = center + (self.p2 -center).matmul(M) + self.p3 = center + (self.p3 -center).matmul(M) + + self.n31 = self.n31.matmul(M) + self.n12 = self.n12.matmul(M) + self.n23 = self.n23.matmul(M) + self.normal = self.normal.matmul(M) + self.centroid = center + (self.centroid-center).matmul(M) + + def get_uv(self, hit): + M_C = hit.point - self.center + u = ((self.pu.dot(M_C)/self.w + 1 ) /2 + self.uv_shift[0]) + v = ((self.pv.dot(M_C)/self.h + 1 ) /2 + self.uv_shift[1]) + return u,v + + + def get_Normal(self, hit): + return self.normal diff --git a/code/sightpy/geometry/triangle_mesh.py b/code/sightpy/geometry/triangle_mesh.py new file mode 100644 index 00000000..7e95d55a --- /dev/null +++ b/code/sightpy/geometry/triangle_mesh.py @@ -0,0 +1,39 @@ +import numpy as np +from ..utils.constants import * +from ..utils.vector3 import vec3 +from ..geometry import Primitive, Triangle_Collider + + +# WORK IN PROGRESS. +# We need to implement a bounding volume hierarchy to make TriangleMesh collision efficient. +# Without a bounding volume hierarchy a model with 200 triangles takes around 3 minutes to be rendered + +class TriangleMesh(Primitive): + def __init__(self,file_name, center, material, max_ray_depth,shadow = True): + super().__init__(center, material,max_ray_depth, shadow = shadow) + self.collider_list += [] + vs = [] + fs = [] + with open(file_name, 'r') as f: + r = f.read() + r = r.split('\n') + for i in r: + i = i.split() + if not i: + continue + elif i[0] == 'v': + x = float(i[1]) + y = float(i[2]) + z = float(i[3]) + vs.append(vec3(x, y, z)) + elif i[0] == 'f': + f1 = int(i[1].split('/')[0]) - 1 + f2 = int(i[2].split('/')[0]) - 1 + f3 = int(i[3].split('/')[0]) - 1 + fs.append([f1, f2, f3]) + for i in fs: + p1 = vs[i[0]] + center + p2 = vs[i[1]] + center + p3 = vs[i[2]] + center + self.collider_list += [colliders.Triangle_Collider(assigned_primitive = self, p1 =p1, p2 = p2, p3 = p3)] + diff --git a/code/sightpy/lights.py b/code/sightpy/lights.py new file mode 100644 index 00000000..3335da13 --- /dev/null +++ b/code/sightpy/lights.py @@ -0,0 +1,49 @@ +from .utils.constants import SKYBOX_DISTANCE +import numpy as np +from abc import abstractmethod + + +# lights only have effect on Glossy materials +class Light: + def __init__(self, pos, color): + self.pos = pos + self.color = color + + @abstractmethod + def get_L(self): + pass + + @abstractmethod + def get_irradiance(self, dist_light, NdotL): + pass + + @abstractmethod + def get_distance(self, M): + pass + + +class PointLight(Light): + def __init__(self, pos, color): + self.pos = pos + self.color = color + def get_L(self): + return (self.pos - M)*(1./(dist_light)) + + def get_distance(self, M): + return np.sqrt((self.pos - M).dot(self.pos - M)) + + def get_irradiance(self,dist_light, NdotL): + return self.color * NdotL/(dist_light**2.) * 100 + +class DirectionalLight(Light): + def __init__(self, Ldir, color): + self.Ldir = Ldir + self.color = color + def get_L(self): + return self.Ldir + + def get_distance(self, M): + return SKYBOX_DISTANCE + + def get_irradiance(self, dist_light, NdotL): + return self.color * NdotL \ No newline at end of file diff --git a/code/sightpy/materials/__init__.py b/code/sightpy/materials/__init__.py new file mode 100644 index 00000000..88b13c0b --- /dev/null +++ b/code/sightpy/materials/__init__.py @@ -0,0 +1,7 @@ +from .material import Material + +from .glossy import Glossy +from .refractive import Refractive +from .thin_film_interference import ThinFilmInterference +from .diffuse import Diffuse +from .emissive import Emissive \ No newline at end of file diff --git a/code/sightpy/materials/diffuse.py b/code/sightpy/materials/diffuse.py new file mode 100644 index 00000000..1a11b3bc --- /dev/null +++ b/code/sightpy/materials/diffuse.py @@ -0,0 +1,99 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3, rgb, extract +from ..utils.random import spherical_caps_pdf, cosine_pdf, mixed_pdf +from functools import reduce as reduce +from ..ray import Ray, get_raycolor +from .. import lights +import numpy as np +from . import Material +from ..textures import * + +class Diffuse(Material): + def __init__(self, diff_color, diffuse_rays = 20, ambient_weight = 0.5, **kwargs): + super().__init__(**kwargs) + + if isinstance(diff_color, vec3): + self.diff_texture = solid_color(diff_color) + elif isinstance(diff_color, texture): + self.diff_texture = diff_color + + self.diffuse_rays = diffuse_rays + self.max_diffuse_reflections = 2 + self.ambient_weight = ambient_weight + + def get_color(self, scene, ray, hit): + + hit.point = (ray.origin + ray.dir * hit.distance) # intersection point + N = hit.material.get_Normal(hit) # normal + + diff_color = self.diff_texture.get_color(hit) + + + color = rgb(0.,0.,0.) + + if ray.diffuse_reflections < 1: + + + nudged = hit.point + N * .000001 + N_repeated = N.repeat(self.diffuse_rays) + + + if ray.n.shape() == 1: + n_repeated = ray.n + else: + n_repeated = ray.n.repeat(self.diffuse_rays) + + nudged_repeated = nudged.repeat(self.diffuse_rays) + hit_repeated = hit.point.repeat(self.diffuse_rays) + + + size = N.shape()[0] * self.diffuse_rays + + pdf1 = cosine_pdf(size, N_repeated) + pdf2 = spherical_caps_pdf(size, nudged_repeated, scene.importance_sampled_list) + + s_pdf = None + if scene.importance_sampled_list == []: + s_pdf = cosine_pdf(size, N_repeated) + else: + s_pdf = mixed_pdf(size, pdf1, pdf2, self.ambient_weight) + + ray_dir = s_pdf.generate() + PDF_val = s_pdf.value(ray_dir) + + NdotL = np.clip(ray_dir.dot(N_repeated),0.,1.) + color_temp = get_raycolor(Ray(nudged_repeated, ray_dir, ray.depth + 1, n_repeated, ray.reflections + 1, ray.transmissions, ray.diffuse_reflections + 1), scene) + color_temp = color_temp * NdotL / PDF_val / (np.pi) # diff_color/np.pi = Lambertian BRDF + color += diff_color * color_temp.reshape(N.shape()[0], self.diffuse_rays).mean(axis = 1) + + return color + + elif ray.diffuse_reflections < self.max_diffuse_reflections: + + """ + when ray.diffuse_reflections > 1 we just call one diffuse ray to solve rendering equation (otherwise is too slow) + """ + + nudged = hit.point + N * .000001 + size = N.shape()[0] + s_pdf = None + + pdf1 = cosine_pdf(size, N) + pdf2 = spherical_caps_pdf(size, nudged, scene.importance_sampled_list) + + if scene.importance_sampled_list == []: + s_pdf = cosine_pdf(size, N) + else: + s_pdf = mixed_pdf(size, pdf1, pdf2, self.ambient_weight) + + ray_dir = s_pdf.generate() + PDF_val = s_pdf.value(ray_dir) + + NdotL = np.clip(N.dot(ray_dir),0.,1.) + color_temp = diff_color * get_raycolor(Ray(nudged, ray_dir, ray.depth + 1, ray.n, ray.reflections + 1, ray.transmissions, ray.diffuse_reflections + 1), scene) + color = color_temp * NdotL / PDF_val / (np.pi) + + return color + + else: + return color diff --git a/code/sightpy/materials/emissive.py b/code/sightpy/materials/emissive.py new file mode 100644 index 00000000..3f375c64 --- /dev/null +++ b/code/sightpy/materials/emissive.py @@ -0,0 +1,23 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3, rgb, extract +from functools import reduce as reduce +from ..ray import Ray, get_raycolor +from .. import lights +import numpy as np +from . import Material +from ..textures import * + +class Emissive(Material): + def __init__(self, color, **kwargs): + + if isinstance(color, vec3): + self.texture_color = solid_color(color) + elif isinstance(color, texture): + self.texture_color = color + + super().__init__(**kwargs) + + + def get_color(self, scene, ray, hit): + diff_color = self.texture_color.get_color(hit) + return diff_color \ No newline at end of file diff --git a/code/sightpy/materials/glossy.py b/code/sightpy/materials/glossy.py new file mode 100644 index 00000000..98162167 --- /dev/null +++ b/code/sightpy/materials/glossy.py @@ -0,0 +1,91 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3, rgb, extract +from functools import reduce as reduce +from ..ray import Ray, get_raycolor +from .. import lights +import numpy as np +from . import Material +from ..textures import * + +class Glossy(Material): + def __init__(self, diff_color, roughness, spec_coeff, diff_coeff, n, **kwargs): + super().__init__(**kwargs) + + if isinstance(diff_color, vec3): + self.diff_texture = solid_color(diff_color) + elif isinstance(diff_color, texture): + self.diff_texture = diff_color + + self.roughness = roughness + self.diff_coeff = diff_coeff + self.spec_coeff = spec_coeff + self.n = n # index of refraction + + + def get_color(self, scene, ray, hit): + + hit.point = (ray.origin + ray.dir * hit.distance) # intersection point + N = hit.material.get_Normal(hit) # normal + + diff_color = self.diff_texture.get_color(hit)* self.diff_coeff + + # Ambient + color = scene.ambient_color * diff_color + V = ray.dir*-1. + nudged = hit.point + N * .000001 # M nudged to avoid itself + + for light in scene.Light_list: + + L = light.get_L() # direction to light + dist_light = light.get_distance(hit.point) # distance to light + NdotL = np.maximum(N.dot(L), 0.) + lv = light.get_irradiance(dist_light, NdotL) # amount of intensity that falls on the surface + + + # direction to ray origin + + H = (L + V).normalize() # Half-way vector + + + + + # Shadow: find if the point is shadowed or not. + # This amounts to finding out if M can see the light + # Shoot a ray from M to L and check what object is the nearest + if not scene.shadowed_collider_list == []: + inters = [s.intersect(nudged, L) for s in scene.shadowed_collider_list] + light_distances, light_hit_orientation = zip(*inters) + light_nearest = reduce(np.minimum, light_distances) + seelight = (light_nearest >= dist_light) + else: + seelight = 1. + + # Lambert shading (diffuse) + color += diff_color * lv * seelight + + if self.roughness != 0.0: + #Fresnel Factor for specular highlight (Schlick’s approximation) + F0 = np.abs((ray.n - self.n)/(ray.n + self.n))**2 + cosθ = np.clip(V.dot(H), 0.0, 1.) + F = F0 + (1. - F0)*(1.- cosθ)**5 + + + # Phong shading (specular highlight) + a = 2./(self.roughness**2.) - 2. + Dphong = np.power(np.clip(N.dot(H), 0., 1.), a) * (a + 2.)/(2.*np.pi) + + # Cook-Torrance model + color += F * Dphong /(4. * np.clip(N.dot(V) * NdotL, 0.001, 1.) ) * seelight * lv * self.spec_coeff + + # Reflection + if ray.depth < hit.surface.max_ray_depth: + + # Fresnel Factor for reflections (Schlick’s approximation) + + F0 = np.abs((scene.n - self.n)/(scene.n + self.n))**2 + cosθ = np.clip(V.dot(N),0.0,1.) + F = F0 + (1. - F0)*(1.- cosθ)**5 + reflected_ray_dir = (ray.dir - N * 2. * ray.dir.dot(N)).normalize() + color += (get_raycolor(Ray(nudged, reflected_ray_dir, ray.depth + 1, ray.n, ray.reflections + 1, ray.transmissions, ray.diffuse_reflections), scene))*F + + return color \ No newline at end of file diff --git a/code/sightpy/materials/material.py b/code/sightpy/materials/material.py new file mode 100644 index 00000000..ef6b050c --- /dev/null +++ b/code/sightpy/materials/material.py @@ -0,0 +1,34 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3, rgb, extract +from functools import reduce as reduce +from ..ray import Ray, get_raycolor +from .. import lights +from ..utils.image_functions import load_image, load_image_as_linear_sRGB +import numpy as np +from abc import abstractmethod + +class Material(): + def __init__(self,normalmap = None): + + if normalmap != None: + normalmap = load_image("sightpy/normalmaps/" + normalmap) + self.normalmap = normalmap + + + def get_Normal(self, hit): + N_coll = hit.collider.get_Normal(hit) + if self.normalmap is not None: + u,v = hit.get_uv() + im = self.normalmap[-((v * self.normalmap.shape[0]*self.repeat ).astype(int)% self.normalmap.shape[0]) , (u * self.normalmap.shape[1]*self.repeat).astype(int) % self.normalmap.shape[1] ].T + N_map = (vec3(im[0] - 0.5,im[1] - 0.5,im[2] - 0.5)) * 2.0 + return N_map.matmul(hit.collider.inverse_basis_matrix).normalize()*hit.orientation + else: + return N_coll*hit.orientation + + def set_normalmap(self, normalmap,repeat= 1.0): + self.normalmap = load_image("sightpy/normalmaps/" + normalmap) + self.repeat = repeat + + @abstractmethod + def get_color(self, scene, ray, hit): + pass diff --git a/code/sightpy/materials/refractive.py b/code/sightpy/materials/refractive.py new file mode 100644 index 00000000..c4524258 --- /dev/null +++ b/code/sightpy/materials/refractive.py @@ -0,0 +1,85 @@ +from ..utils.constants import * +from ..utils.vector3 import vec3, rgb, extract +from functools import reduce as reduce +from ..ray import Ray, get_raycolor +from .. import lights +import numpy as np +from . import Material + +class Refractive(Material): + def __init__(self, n, **kwargs): + super().__init__(**kwargs) + + self.n = n # index of refraction + + # Instead of defining a index of refraction (n) for each wavelenght (computationally expensive) we aproximate defining the index of refraction + # using a vec3 for red = 630 nm, green 555 nm, blue 475 nm, the most sensitive wavelenghts of human eye. + + # Index a refraction is a complex number. + # The real part is involved in how much light is reflected and model refraction direction via Snell Law. + # The imaginary part of n is involved in how much light is reflected and absorbed. For non-transparent materials like metals is usually between (0.1j,3j) + # and for transparent materials like glass is usually between (0.j , 1e-7j) + + + + def get_color(self, scene, ray, hit): + + hit.point = (ray.origin + ray.dir * hit.distance) # intersection point + N = hit.material.get_Normal(hit) # normal + + color = rgb(0.,0.,0.) + + V = ray.dir*-1. # direction to ray origin + nudged = hit.point + N * .000001 # M nudged to avoid itself + # compute reflection and refraction + # a paper explaining formulas used: + # https://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--reflection_refraction.pdf + # reallistic refraction is expensive. (requires exponential complexity because each ray is divided in two) + + if ray.depth 1.0) by scaling + rgb_max = np.amax(rgb, axis=0) + 0.00001 # avoid division by zero + intensity_cutoff = 1.0 + rgb = np.where(rgb_max > intensity_cutoff, rgb * intensity_cutoff / (rgb_max), rgb) + + return rgb + + +def sRGB_to_sRGB_linear(rgb): + + '''sRGB standard for gamma inverse correction.''' + rgb_linear = np.where( rgb <= 0.03928, rgb / 12.92, np.power((rgb + 0.055) / 1.055, 2.4)) + + return rgb_linear + diff --git a/code/sightpy/utils/constants.py b/code/sightpy/utils/constants.py new file mode 100644 index 00000000..0c39f367 --- /dev/null +++ b/code/sightpy/utils/constants.py @@ -0,0 +1,5 @@ +UPWARDS = 1 +UPDOWN = -1 +FARAWAY = 1.0e39 +SKYBOX_DISTANCE = 1.0e6 + diff --git a/code/sightpy/utils/image_functions.py b/code/sightpy/utils/image_functions.py new file mode 100644 index 00000000..34ea7374 --- /dev/null +++ b/code/sightpy/utils/image_functions.py @@ -0,0 +1,40 @@ +from PIL import Image, ImageFilter +import numpy as np +from pathlib import Path +from .colour_functions import sRGB_to_sRGB_linear + + +def load_image(path): + img = Image.open(Path(path)) + return np.asarray(img)/256. + + +def load_image_with_blur(path, blur = 0.): + img = Image.open(Path(path)) + img = img.filter(ImageFilter.GaussianBlur(radius=blur)) + + return np.asarray(img)/256. + + + + +def load_image_as_linear_sRGB(path, blur = 0.0): + + path = Path(path) + location = str(path.parents[0]) + name = str(path.name)#be sure that image doesn't lose quality + + print("proccesing " + name) + img = Image.open(path) + + if blur != 0.0: + img = img.filter(ImageFilter.GaussianBlur(radius=blur)) + + img_array = np.asarray(img)/256. + img_sRGB_linear_array = sRGB_to_sRGB_linear(img_array) + return img_sRGB_linear_array + + + + + diff --git a/code/sightpy/utils/random.py b/code/sightpy/utils/random.py new file mode 100644 index 00000000..0166fc8b --- /dev/null +++ b/code/sightpy/utils/random.py @@ -0,0 +1,231 @@ +import numpy as np +from ..utils.vector3 import vec3 +from abc import abstractmethod + +def random_in_unit_disk(shape): + r = np.sqrt(np.random.rand(shape)) + phi = np.random.rand(shape)*2*np.pi + return r * np.cos(phi), r * np.sin(phi) + +def random_in_unit_sphere(shape): + + #https://mathworld.wolfram.com/SpherePointPicking.html + phi = np.random.rand(shape)*2*np.pi + u = 2.*np.random.rand(shape) - 1. + r = np.sqrt(1-u**2) + return vec3( r*np.cos(phi), r*np.sin(phi), u) + + + +class PDF: + """Probability density function""" + @abstractmethod + def value(self,ray_dir): + """get probability density function value at direction ray_dir""" + pass + + @abstractmethod + def generate(self): + """generate random ray directions according the probability density function""" + pass + + + + +class hemisphere_pdf(PDF): + """Probability density Function""" + def __init__(self,shape, normal): + self.shape = shape + self.normal = normal + + + def value(self,ray_dir): + return 1./(2.*np.pi) + + def generate(self): + r = random_in_unit_sphere(self.shape) + return vec3.where( self.normal.dot(r) < 0. , r*-1., r ) + + +class cosine_pdf(PDF): + """Probability density Function""" + def __init__(self,shape, normal): + self.shape = shape + self.normal = normal + + + + def value(self,ray_dir): + return np.clip(ray_dir.dot(self.normal),0.,1.)/np.pi + + def generate(self): + ax_w = self.normal + a = vec3.where( np.abs(ax_w.x) > 0.9 , vec3(0,1,0) , vec3(1,0,0)) + ax_v = ax_w.cross(a).normalize() + ax_u = ax_w.cross(ax_v) + + phi = np.random.rand(self.shape)*2*np.pi + r2 = np.random.rand(self.shape) + + z = np.sqrt(1 - r2) + x = np.cos(phi) * np.sqrt(r2) + y = np.sin(phi) * np.sqrt(r2) + + return ax_u*x + ax_v*y + ax_w*z + + + + +class spherical_caps_pdf(PDF): + """Probability density Function""" + def __init__(self,shape, origin, importance_sampled_list): + self.shape = shape + self.origin = origin + self.importance_sampled_list = importance_sampled_list + self.l = len(importance_sampled_list) + + def value(self, ray_dir): + PDF_value = 0. + for i in range(self.l): + PDF_value += np.where( ray_dir.dot(self.ax_w_list[i]) > self.cosθmax_list[i] , 1/((1 - self.cosθmax_list[i])*2*np.pi) , 0. ) + PDF_value = PDF_value/self.l + return PDF_value + + + def generate(self): + shape = self.shape + origin = self.origin + importance_sampled_list = self.importance_sampled_list + l = self.l + + mask = (np.random.rand(shape) * l).astype(int) + mask_list = [None]*l + + cosθmax_list = [None]*l + ax_u_list = [None]*l + ax_v_list = [None]*l + ax_w_list = [None]*l + + for i in range(l): + + ax_w_list[i] = (importance_sampled_list[i].center - origin).normalize() + a = vec3.where( np.abs(ax_w_list[i].x) > 0.9 , vec3(0,1,0) , vec3(1,0,0)) + ax_v_list[i] = ax_w_list[i].cross(a).normalize() + ax_u_list[i] = ax_w_list[i].cross(ax_v_list[i]) + mask_list[i] = mask == i + + + target_distance = np.sqrt((importance_sampled_list[i].center - origin).dot(importance_sampled_list[i].center - origin)) + + cosθmax_list[i] = np.sqrt(1 - np.clip(importance_sampled_list[i].bounded_sphere_radius / target_distance, 0., 1.)**2 ) + + self.cosθmax_list = cosθmax_list + self.ax_w_list = ax_w_list + + phi = np.random.rand(shape)*2*np.pi + r2 = np.random.rand(shape) + + cosθmax = np.select(mask_list, cosθmax_list) + ax_w = vec3.select(mask_list, ax_w_list) + ax_v = vec3.select(mask_list, ax_v_list) + ax_u = vec3.select(mask_list, ax_u_list) + + z = 1. + r2 * (cosθmax - 1.) + x = np.cos(phi) * np.sqrt(1. - z**2) + y = np.sin(phi) * np.sqrt(1. - z**2) + + ray_dir = ax_u*x + ax_v*y + ax_w*z + return ray_dir + + +class mixed_pdf(PDF): + """Probability density Function""" + def __init__(self,shape, pdf1, pdf2, pdf1_weight = 0.5): + + self.pdf1_weight = pdf1_weight + self.pdf2_weight = 1. - pdf1_weight + self.shape = shape + self.pdf1 = pdf1 + self.pdf2 = pdf2 + + + def value(self,ray_dir): + return self.pdf1.value(ray_dir) * self.pdf1_weight + self.pdf2.value(ray_dir) * self.pdf2_weight + + def generate(self): + mask = np.random.rand(self.shape) + return vec3.where( mask < self.pdf1_weight, self.pdf1.generate(), self.pdf2.generate() ) + + + + + + + +def random_in_unit_spherical_caps(shape, origin, importance_sampled_list): + + l = len(importance_sampled_list) + + + mask = (np.random.rand(shape) * l).astype(int) + mask_list = [None]*l + + cosθmax_list = [None]*l + ax_u_list = [None]*l + ax_v_list = [None]*l + ax_w_list = [None]*l + + for i in range(l): + + ax_w_list[i] = (importance_sampled_list[i].center - origin).normalize() + a = vec3.where( np.abs(ax_w_list[i].x) > 0.9 , vec3(0,1,0) , vec3(1,0,0)) + ax_v_list[i] = ax_w_list[i].cross(a).normalize() + ax_u_list[i] = ax_w_list[i].cross(ax_v_list[i]) + mask_list[i] = mask == i + + + target_distance = np.sqrt((importance_sampled_list[i].center - origin).dot(importance_sampled_list[i].center - origin)) + + cosθmax_list[i] = np.sqrt(1 - np.clip(importance_sampled_list[i].bounded_sphere_radius / target_distance, 0., 1.)**2 ) + + + phi = np.random.rand(shape)*2*np.pi + r2 = np.random.rand(shape) + + cosθmax = np.select(mask_list, cosθmax_list) + ax_w = vec3.select(mask_list, ax_w_list) + ax_v = vec3.select(mask_list, ax_v_list) + ax_u = vec3.select(mask_list, ax_u_list) + + z = 1. + r2 * (cosθmax - 1.) + x = np.cos(phi) * np.sqrt(1. - z**2) + y = np.sin(phi) * np.sqrt(1. - z**2) + + ray_dir = ax_u*x + ax_v*y + ax_w*z + + PDF = 0. + for i in range(l): + PDF += np.where( ray_dir.dot(ax_w_list[i]) > cosθmax_list[i] , 1/((1 - cosθmax_list[i])*2*np.pi) , 0. ) + PDF = PDF/l + + return ray_dir, PDF + +def random_in_unit_spherical_cap(shape,cosθmax,normal): + + + ax_w = normal + a = vec3.where( np.abs(ax_w.x) > 0.9 , vec3(0,1,0) , vec3(1,0,0)) + ax_v = ax_w.cross(a).normalize() + ax_u = ax_w.cross(ax_v) + + phi = np.random.rand(shape)*2*np.pi + r2 = np.random.rand(shape) + + z = 1. + r2 * (cosθmax - 1.) + x = np.cos(phi) * np.sqrt(1. - z**2) + y = np.sin(phi) * np.sqrt(1. - z**2) + + + + + return ax_u*x + ax_v*y + ax_w*z \ No newline at end of file diff --git a/code/sightpy/utils/vector3.py b/code/sightpy/utils/vector3.py new file mode 100644 index 00000000..43294adb --- /dev/null +++ b/code/sightpy/utils/vector3.py @@ -0,0 +1,188 @@ +import numpy as np +import numbers + +def extract(cond, x): + if isinstance(x, numbers.Number): + return x + else: + return np.extract(cond, x) + +class vec3(): + + def __init__(self, x, y, z): + self.x = x + self.y = y + self.z = z + + def __str__(self): + # Used for debugging. This method is called when you print an instance + return "(" + str(self.x) + ", " + str(self.y) + ", " + str(self.z) + ")" + + + def __add__(self, v): + if isinstance(v, vec3): + return vec3(self.x + v.x, self.y + v.y, self.z + v.z) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(self.x + v, self.y + v, self.z + v) + def __radd__(self, v): + if isinstance(v, vec3): + return vec3(self.x + v.x, self.y + v.y, self.z + v.z) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(self.x + v, self.y + v, self.z + v) + def __sub__(self, v): + if isinstance(v, vec3): + return vec3(self.x - v.x, self.y - v.y, self.z - v.z) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(self.x - v, self.y - v, self.z - v) + def __rsub__(self, v): + if isinstance(v, vec3): + return vec3(v.x - self.x, v.y - self.y , v.z - self.z) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(v - self.x, v - self.y , v - self.z) + + def __mul__(self, v): + if isinstance(v, vec3): + return vec3(self.x * v.x , self.y * v.y , self.z * v.z ) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(self.x * v, self.y * v, self.z * v) + def __rmul__(self, v): + if isinstance(v, vec3): + return vec3(v.x *self.x , v.y * self.y, v.z * self.z ) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(v *self.x , v * self.y, v * self.z ) + def __truediv__(self, v): + if isinstance(v, vec3): + return vec3(self.x / v.x , self.y / v.y , self.z / v.z ) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(self.x / v, self.y / v, self.z / v) + + + def __rtruediv__(self, v): + if isinstance(v, vec3): + return vec3(v.x / self.x, v.y / self.y, v.z / self.z) + elif isinstance(v, numbers.Number) or isinstance(v, np.ndarray): + return vec3(v / self.x, v / self.y, v / self.z) + + + + def __abs__(self): + return vec3(np.abs(self.x), np.abs(self.y), np.abs(self.z)) + + def real(v): + return vec3(np.real(v.x), np.real(v.y), np.real(v.z)) + + def imag(v): + return vec3(np.imag(v.x), np.imag(v.y), np.imag(v.z)) + + def yzx(self): + return vec3(self.y, self.z, self.x) + def xyz(self): + return vec3(self.x, self.y, self.z) + def zxy(self): + return vec3(self.z, self.x, self.y) + def xyz(self): + return vec3(self.x, self.y, self.z) + + + def average(self): + return (self.x + self.y + self.z)/3 + + def matmul(self, matrix): + if isinstance(self.x, numbers.Number): + return array_to_vec3(np.dot(matrix,self.to_array())) + elif isinstance(self.x, np.ndarray): + return array_to_vec3(np.tensordot(matrix,self.to_array() , axes=([1,0]))) + + def change_basis(self, new_basis): + return vec3(self.dot(new_basis[0]), self.dot(new_basis[1]), self.dot(new_basis[2])) + + def __pow__(self, a): + return vec3(self.x**a, self.y**a, self.z**a) + + def dot(self, v): + return self.x*v.x + self.y*v.y + self.z*v.z + + def exp(v): + return vec3(np.exp(v.x) , np.exp(v.y) ,np.exp(v.z)) + + def sqrt(v): + return vec3(np.sqrt(v.x) , np.sqrt(v.y) ,np.sqrt(v.z)) + + def to_array(self): + return np.array([self.x , self.y , self.z]) + + def cross(self, v): + return vec3(self.y*v.z - self.z*v.y, -self.x*v.z + self.z*v.x, self.x*v.y - self.y*v.x) + + def length(self): + return np.sqrt(self.dot(self)) + + def square_length(self): + return self.dot(self) + + def normalize(self): + mag = self.length() + return self * (1.0 / np.where(mag == 0, 1, mag)) + + + + def components(self): + return (self.x, self.y, self.z) + + def extract(self, cond): + return vec3(extract(cond, self.x), + extract(cond, self.y), + extract(cond, self.z)) + + def where(cond, out_true, out_false): + return vec3(np.where(cond, out_true.x, out_false.x), + np.where(cond, out_true.y, out_false.y), + np.where(cond, out_true.z, out_false.z)) + + def select(mask_list, out_list): + out_list_x = [i.x for i in out_list] + out_list_y = [i.y for i in out_list] + out_list_z = [i.z for i in out_list] + + return vec3(np.select(mask_list, out_list_x), + np.select(mask_list, out_list_y), + np.select(mask_list, out_list_z)) + + def clip(self, min, max): + return vec3(np.clip(self.x, min, max), + np.clip(self.y, min, max), + np.clip(self.z, min, max)) + + def place(self, cond): + r = vec3(np.zeros(cond.shape), np.zeros(cond.shape), np.zeros(cond.shape)) + np.place(r.x, cond, self.x) + np.place(r.y, cond, self.y) + np.place(r.z, cond, self.z) + return r + + def repeat(self, n): + return vec3(np.repeat(self.x , n), np.repeat(self.y , n), np.repeat(self.z , n)) + + def reshape(self, *newshape): + return vec3(self.x.reshape(*newshape), self.y.reshape(*newshape), self.z.reshape(*newshape)) + + def shape(self, *newshape): + if isinstance(self.x, numbers.Number): + return 1 + elif isinstance(self.x, np.ndarray): + return self.x.shape + + def mean(self, axis): + return vec3(np.mean(self.x,axis = axis), np.mean(self.y,axis = axis), np.mean(self.z,axis = axis)) + + def __eq__(self, other): + return (self.x == other.x) & (self.y == other.y) & (self.z == other.z) + + +def array_to_vec3(array): + return vec3(array[0],array[1],array[2]) + + + +global rgb +rgb = vec3 \ No newline at end of file diff --git a/code/sightpy_test.py b/code/sightpy_test.py new file mode 100644 index 00000000..0e9c15e1 --- /dev/null +++ b/code/sightpy_test.py @@ -0,0 +1,130 @@ +from sightpy import * + + +# Set Scene + +Sc = Scene(ambient_color=rgb(0.00, 0.00, 0.00)) + + +angle = -0 + +Sc.add_Camera( + screen_width=100, + screen_height=100, + look_from=vec3(278, 278, 800), + look_at=vec3(278, 278, 0), + focal_distance=1.0, + field_of_view=40, +) + + +# define materials to use + +green_diffuse = Diffuse(diff_color=rgb(0.12, 0.45, 0.15)) +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( + Plane( + material=emissive_white, + center=vec3(213 + 130 / 2, 554, -227.0 - 105 / 2), + width=130.0, + height=105.0, + u_axis=vec3(1.0, 0.0, 0), + v_axis=vec3(0.0, 0, 1.0), + ), + importance_sampled=True, +) + + +Sc.add( + Plane( + material=white_diffuse, + center=vec3(555 / 2, 555 / 2, -555.0), + width=555.0, + height=555.0, + u_axis=vec3(0.0, 1.0, 0), + v_axis=vec3(1.0, 0, 0.0), + ) +) + +Sc.add( + Plane( + material=green_diffuse, + center=vec3(-0.0, 555 / 2, -555 / 2), + width=555.0, + height=555.0, + u_axis=vec3(0.0, 1.0, 0), + v_axis=vec3(0.0, 0, -1.0), + ) +) + +Sc.add( + Plane( + material=red_diffuse, + center=vec3(555.0, 555 / 2, -555 / 2), + width=555.0, + height=555.0, + u_axis=vec3(0.0, 1.0, 0), + v_axis=vec3(0.0, 0, -1.0), + ) +) + +Sc.add( + Plane( + material=white_diffuse, + center=vec3(555 / 2, 555, -555 / 2), + width=555.0, + height=555.0, + u_axis=vec3(1.0, 0.0, 0), + v_axis=vec3(0.0, 0, -1.0), + ) +) + +Sc.add( + Plane( + material=white_diffuse, + center=vec3(555 / 2, 0.0, -555 / 2), + width=555.0, + height=555.0, + u_axis=vec3(1.0, 0.0, 0), + v_axis=vec3(0.0, 0, -1.0), + ) +) + + +cb = Cuboid( + material=white_diffuse, + center=vec3(182.5, 165, -285 - 160 / 2), + width=165, + height=165 * 2, + length=165, + shadow=False, +) +cb.rotate(θ=15, u=vec3(0, 1, 0)) +Sc.add(cb) + + +Sc.add( + Sphere( + material=blue_glass, + center=vec3(370.5, 165 / 2, -65 - 185 / 2), + radius=165 / 2, + shadow=False, + max_ray_depth=3, + ), + importance_sampled=True, +) +# Render + +img = Sc.render(samples_per_pixel=100, progress_bar=True) +# you are going to need more than 10 samples to remove the noise. At least 1000 for a nice image. + +img.save("cornell_box.png") + +img.show()