sightpy

code/requirements.txt

@@ -6,4 +6,6 @@ matplotlib
 flake8
 black
 autopep8
-flake8-max-function-length
+flake8-max-function-length
+pillow
+progressbar

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 *
+
+
+
+
+

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()

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)

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)
+

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)
+

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)
+

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 *
+
+
+
+
+
+

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

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
+
+
+
+
+
+

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
+

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)
+            

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

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)
+            

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

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)]
+

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

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

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

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

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

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

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 <hit.surface.max_ray_depth:
+
+            """
+            if hit_orientation== UPWARDS:
+               #ray enter in the material
+            if hit_orientation== UPDOWN:
+               #ray get out of the material   
+            """
+            n1 = ray.n
+            n2 = vec3.where(hit.orientation== UPWARDS,self.n,scene.n)
+
+            n1_div_n2 =  vec3.real(n1)/vec3.real(n2) 
+            cosθi = V.dot(N)
+            sin2θt = (n1_div_n2)**2 * (1.-cosθi**2)
+
+            # compute complete fresnel term
+            cosθt = vec3.sqrt(1. - (n1/n2)**2 * (1.-cosθi**2)  )
+            r_per = (n1*cosθi - n2*cosθt)/(n1*cosθi + n2*cosθt)
+            r_par = -1.*(n1*cosθt - n2*cosθi)/(n1*cosθt + n2*cosθi) 
+            F = (np.abs(r_per)**2 + np.abs(r_par)**2)/2.
+            
+            # compute reflection
+            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
+            
+
+
+            # compute refraction 
+            # Spectrum dispersion is not implemented. 
+            # We approximate refraction direction averaging index of refraction of each wavelenght
+            n1_div_n2_aver = n1_div_n2.average()
+            sin2θt = (n1_div_n2_aver)**2 * (1.-cosθi**2)
+
+            non_TiR = (sin2θt <= 1.)
+            if np.any(non_TiR): # avoid total internal reflection
+
+                refracted_ray_dir = (ray.dir*(n1_div_n2_aver) + N*(n1_div_n2_aver * cosθi - np.sqrt(1-np.clip(sin2θt,0,1)))).normalize() 
+                nudged = hit.point - N * .000001  #nudged for refraction
+                T = 1. - F
+                refracted_color = (get_raycolor( Ray(nudged, refracted_ray_dir, ray.depth + 1, n2, ray.reflections, ray.transmissions + 1, ray.diffuse_reflections).extract(non_TiR), scene)  )*T.extract(non_TiR) 
+                color += refracted_color.place(non_TiR)
+            
+                
+            # absorption:
+            # approximation using wavelength for red = 630 nm, green 550 nm, blue 475 nm
+            color = color *vec3.exp(-2.*vec3.imag(ray.n)*2.*np.pi/vec3(630,550,475) * 1e9* hit.distance)
+        return color

code/sightpy/materials/thin_film_interference.py

@@ -0,0 +1,67 @@
+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 ..utils.image_functions import load_image
+
+class ThinFilmInterference(Material):
+    def __init__(self, thickness, noise = 0.,**kwargs):
+        super().__init__(**kwargs)
+        self.thickness = thickness
+
+        # precomputed reflectance vs cosθI (vertical axis) and thickness (horizontal axis)
+        self.thin_film_interference_reflectance = load_image("sightpy/textures/thin_film_interference_n=1.4.png")
+        self.thickness_noise = load_image("sightpy/textures/noise.png")
+        self.thickness_noise = (self.thickness_noise[:,:,0])
+        self.noise_factor = noise
+
+    def get_color(self, scene, ray, hit):
+        
+        hit.point = (ray.origin + ray.dir * hit.distance) # intersection point
+        N = hit.material.get_Normal(hit)     # normal 
+
+        # Ambient
+        color = rgb(0.,0.,0.)
+
+        V = ray.dir*-1.                            # direction to ray origin
+        
+
+        if ray.depth < hit.surface.max_ray_depth:
+
+            """
+            if hit_orientation== UPWARDS:
+               #ray enter in the material
+            if hit_orientation== UPDOWN:
+               #ray get out of the material   
+            """
+
+            cosθi = V.dot(N)
+
+            u,v = hit.get_uv()
+            thickness = self.thickness 
+            if self.noise_factor != 0.:
+                thickness += self.noise_factor*(self.thickness_noise[-((v * self.thickness_noise.shape[0]*0.5 ).astype(int)% self.thickness_noise.shape[0]) , (u   * self.thickness_noise.shape[1]*0.5).astype(int) % self.thickness_noise.shape[1]  ].T - 0.5)
+                Fim = self.thin_film_interference_reflectance[(cosθi* self.thin_film_interference_reflectance.shape[0]).astype(int),    thickness.astype(int)     ]
+            else:
+                Fim = self.thin_film_interference_reflectance[(cosθi* self.thin_film_interference_reflectance.shape[0]).astype(int),    int(thickness)    ]
+            
+            F = vec3(Fim[:,0],Fim[:,1],Fim[:,2])
+            # compute reflection
+            reflected_ray_dir = (ray.dir - N * 2. * ray.dir.dot(N)).normalize()
+
+            nudged = hit.point + N * .000001                  # M nudged to avoid itself
+            color += (scene.ambient_color + get_raycolor(Ray(nudged, reflected_ray_dir, ray.depth + 1, ray.n, ray.reflections + 1, ray.transmissions,  ray.diffuse_reflections), scene))*F
+            
+
+
+            # because the film is very thin (nm) we ignore refraction for transmitted ray.
+
+            transmitted_ray_dir = ray.dir 
+            nudged = hit.point - N * .000001  #nudged for transmitted ray
+            T = 1. - F
+            transmitted_color = get_raycolor( Ray(nudged, transmitted_ray_dir, ray.depth + 1, ray.n, ray.reflections, ray.transmissions + 1,  ray.diffuse_reflections), scene  )*T
+            color += transmitted_color
+        return color

code/sightpy/ray.py

@@ -0,0 +1,103 @@
+from .utils.constants import *
+from .utils.vector3 import vec3, extract, rgb
+import numpy as np
+from functools import reduce as reduce 
+
+
+class Ray:    
+    """Info of the ray and the media it's travelling"""
+    def __init__(self,origin, dir, depth, n, reflections, transmissions, diffuse_reflections):
+
+        self.origin  = origin   # the point where the ray comes from
+        self.dir = dir          # direction of the ray
+        self.depth = depth     # ray_depth is the number of the refrections + transmissions/refractions, starting at zero for camera rays
+        self.n = n             # ray_n is the index of refraction of the media in which the ray is travelling
+        
+        # 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)
+
+
+        self.reflections = reflections  #reflections is the number of the refrections, starting at zero for camera rays
+        self.transmissions = transmissions  #transmissions is the number of the transmissions/refractions, starting at zero for camera rays
+        self.diffuse_reflections = diffuse_reflections  #reflections is the number of the refrections, starting at zero for camera rays
+
+
+    def extract(self,hit_check):
+        return Ray(self.origin.extract(hit_check), self.dir.extract(hit_check), self.depth,  self.n.extract(hit_check), self.reflections, self.transmissions,self.diffuse_reflections)
+
+
+class Hit:  
+    """Info of the ray-surface intersection"""
+    def __init__(self, distance, orientation, material, collider,surface):
+        self.distance = distance
+        self.orientation = orientation 
+        self.material = material
+        self.collider = collider
+        self.surface = surface
+        self.u = None
+        self.v = None
+        self.N = None
+        self.point = None 
+
+    def get_uv(self):
+        if self.u is None: #this is for prevent multiple computations of u,v
+            self.u, self.v = self.collider.assigned_primitive.get_uv(self)
+        return self.u, self.v
+
+    def get_normal(self):
+        if self.N is None: #this is for prevent multiple computations of normal
+            self.N = self.collider.get_N(self)
+        return self.N
+
+
+
+def get_raycolor(ray, scene):
+
+    inters = [s.intersect(ray.origin, ray.dir) for s in scene.collider_list]
+    distances, hit_orientation = zip(*inters)
+
+    
+    # get the shortest distance collision
+    nearest = reduce(np.minimum, distances)
+    color = rgb(0., 0., 0.)
+
+    
+    for (coll, dis , orient) in zip(scene.collider_list, distances, hit_orientation):
+        hit_check = (nearest != FARAWAY) & (dis == nearest)
+
+        if np.any(hit_check):  
+
+            material = coll.assigned_primitive.material
+            hit_info = Hit(extract(hit_check,dis) , extract(hit_check,orient), material, coll, coll.assigned_primitive)
+
+            
+            
+            cc = material.get_color(scene,  ray.extract(hit_check), hit_info)
+            color += cc.place(hit_check)
+
+    return color
+    
+
+
+
+def get_distances(ray, scene): #Used for debugging ray-surface collisions. Return a grey map of objects distances.
+
+    inters = [s.intersect(ray.origin, ray.dir) for s in scene.collider_list]
+    distances, hit_orientation = zip(*inters)
+    # get the shortest distance collision
+    nearest = reduce(np.minimum, distances)
+    
+    max_r_distance = 10
+    r_distance = np.where(nearest <= max_r_distance, nearest, max_r_distance)
+    norm_r_distance = r_distance/max_r_distance
+    return rgb(norm_r_distance, norm_r_distance, norm_r_distance)
+
+
+    
+    
+

code/sightpy/scene.py

@@ -0,0 +1,114 @@
+from PIL import Image
+import numpy as np
+import time
+from .utils import colour_functions as cf
+from .camera import Camera
+from .utils.constants import *
+from .utils.vector3 import vec3, rgb
+from .ray import Ray, get_raycolor, get_distances
+from . import lights
+from .backgrounds.skybox import SkyBox
+from .backgrounds.panorama import Panorama
+
+
+class Scene():
+    def __init__(self, ambient_color = rgb(0.01, 0.01, 0.01), n = vec3(1.0,1.0,1.0)) :
+        # n = index of refraction (by default index of refraction of air n = 1.)
+        
+        self.scene_primitives = []
+        self.collider_list = []
+        self.shadowed_collider_list = []
+        self.Light_list = []
+        self.importance_sampled_list = []
+        self.ambient_color = ambient_color
+        self.n = n
+        self.importance_sampled_list = []
+    def add_Camera(self, look_from, look_at, **kwargs):
+        self.camera = Camera(look_from, look_at, **kwargs)
+
+
+    def add_PointLight(self, pos, color):
+        self.Light_list += [lights.PointLight(pos, color)]
+        
+    def add_DirectionalLight(self, Ldir, color):
+        self.Light_list += [lights.DirectionalLight(Ldir.normalize() , color)]  
+
+    def add(self,primitive, importance_sampled = False):
+        self.scene_primitives += [primitive]
+        self.collider_list += primitive.collider_list
+
+        if importance_sampled == True:
+            self.importance_sampled_list += [primitive]
+
+        if primitive.shadow == True:
+            self.shadowed_collider_list += primitive.collider_list
+            
+        
+    def add_Background(self, img, light_intensity = 0.0, blur =0.0 , spherical = False):
+
+        primitive = None
+        if spherical == False:
+            primitive = SkyBox(img, light_intensity = light_intensity, blur = blur)
+        else:
+            primitive = Panorama(img, light_intensity = light_intensity, blur = blur)
+
+        self.scene_primitives += [primitive]        
+        self.collider_list += primitive.collider_list
+
+        
+    def render(self, samples_per_pixel, progress_bar = False):
+
+        print ("Rendering...")
+
+        t0 = time.time()
+        color_RGBlinear = rgb(0.,0.,0.)
+
+        if progress_bar == True:
+
+
+            try:
+                import progressbar
+
+            except ModuleNotFoundError:
+                 print("progressbar module is required. \nRun: pip install progressbar")
+            
+            bar = progressbar.ProgressBar()
+            for i in bar(range(samples_per_pixel)):
+                color_RGBlinear += get_raycolor(self.camera.get_ray(self.n), scene = self)
+                bar.update(i)
+        else:
+
+
+            for i in range(samples_per_pixel):
+                color_RGBlinear += get_raycolor(self.camera.get_ray(self.n), scene = self)
+                
+
+
+        #average samples per pixel (antialiasing)
+        color_RGBlinear = color_RGBlinear/samples_per_pixel
+        #gamma correction
+        color = cf.sRGB_linear_to_sRGB(color_RGBlinear.to_array())
+        
+        print ("Render Took", time.time() - t0)
+
+        img_RGB = []
+        for c in color:
+            # average ray colors that fall in the same pixel. (antialiasing) 
+            img_RGB += [Image.fromarray((255 * np.clip(c, 0, 1).reshape((self.camera.screen_height, self.camera.screen_width))).astype(np.uint8), "L") ]
+
+        return Image.merge("RGB", img_RGB)
+
+
+    def get_distances(self): #Used for debugging ray-primitive collisions. Return a grey map of objects distances.
+
+        print ("Rendering...")
+        t0 = time.time()
+        color_RGBlinear = get_distances( self.camera.get_ray(self.n), scene = self)
+        #gamma correction
+        color = color_RGBlinear.to_array()
+        
+        print ("Render Took", time.time() - t0)
+
+
+        img_RGB = [Image.fromarray((255 * np.clip(c, 0, 1).reshape((self.camera.screen_height, self.camera.screen_width))).astype(np.uint8), "L") for c in color]
+        return Image.merge("RGB", img_RGB)

code/sightpy/textures/__init__.py

@@ -0,0 +1 @@
+from .texture import texture, image, solid_color

code/sightpy/textures/texture.py

@@ -0,0 +1,38 @@
+from ..utils.constants import *
+from ..utils.vector3 import vec3, rgb
+from ..ray import Ray, get_raycolor
+from ..utils.image_functions import load_image, load_image_as_linear_sRGB
+import numpy as np
+from abc import abstractmethod 
+
+class texture():
+
+    @abstractmethod   
+    def __init__(self):
+    	pass
+
+    @abstractmethod  
+    def get_color(self, hit):
+    	pass
+
+
+class solid_color(texture):
+	 
+    def __init__(self,color):
+    	self.color = color
+
+    def get_color(self, hit):
+    	return self.color
+
+
+class image(texture):
+	 
+    def __init__(self,img, repeat = 1.0):
+    	self.img = load_image_as_linear_sRGB("sightpy/textures/" + img)
+    	self.repeat = repeat
+
+    def get_color(self, hit):
+        u,v = hit.get_uv()
+        im = self.img[-((v * self.img.shape[0]*self.repeat ).astype(int)% self.img.shape[0]) , (u   * self.img.shape[1]*self.repeat).astype(int) % self.img.shape[1]  ].T
+        color = vec3(im[0],im[1],im[2])
+        return color

code/sightpy/utils/colour_functions.py

@@ -0,0 +1,22 @@
+import numpy as np
+
+def sRGB_linear_to_sRGB(rgb_linear):
+
+    '''sRGB standard for gamma inverse correction.'''
+    rgb = np.where( rgb_linear <= 0.00304,  12.92 * rgb_linear,  1.055 * np.power(rgb_linear, 1.0/2.4) - 0.055)
+    
+    # clip intensity if needed (rgb values > 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
+

code/sightpy/utils/constants.py

@@ -0,0 +1,5 @@
+UPWARDS = 1
+UPDOWN = -1         
+FARAWAY = 1.0e39  
+SKYBOX_DISTANCE = 1.0e6   
+

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
+
+
+
+
+

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

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

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()