Untitled

from PIL import Image
import numpy as np
import math
import colorsys

#Define a vector class
class Vector:
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z
    
    def __mul__(self, a):
        return Vector(a*self.x, a*self.y, a*self.z)
    
    def __truediv__(self, a):
        return Vector(self.x/a, self.y/a, self.z/a)
    
    def __add__(self, other):
        return Vector(self.x + other.x, self.y + other.y, self.z + other.z)
        
    def __sub__(self, other):
        return Vector(self.x - other.x, self.y - other.y, self.z - other.z)
    
    def dot_product(self, other):
        return (self.x*other.x + self.y*other.y + self.z*other.z)*1.0
    
    def abs(self):
        return math.sqrt(self.dot_product(self))
    
    def normalize(self):
        return self / self.abs()
    
    def cos_angle(self, other):
        return self.dot_product(other)/(self.abs() * other.abs())
    
    def sin_square_angle(self, other):
        cos = self.cos_angle(other)
        
        return 1 - cos**2
    
    def print(self):
        print((self.x, self.y, self.z))

#Define a sphere class
class Sphere:
    def __init__(self, radio, x = 0, y = 0, z = 0):
        self.center = Vector(x, y, z)
        self.radio = radio
    
    def normal_vector(self, point):
        v = point - self.center
        return v.normalize()
    
    def ray_casting(self, pixel):
        v = (pixel - camera).normalize()
        position = camera - self.center
        
        delta_over4 = (v.dot_product(position) ** 2) - (position.dot_product(position) - (self.radio ** 2))
        #delta_over4 = self.radio**2 - position.dot_product(position)*v.sin_square_angle(position)
        
        if delta_over4 < 0:
            #print("Negative")
            return False
        elif delta_over4 == 0:
            #print("Zero")
            t = - v.dot_product(position)
    
            if t >= 0:
                return camera + v*t
            else:
                return False
        else:
            #print("Positive")
            t1 = - (v.dot_product(position) + math.sqrt(delta_over4))
            t2 = math.sqrt(delta_over4) - v.dot_product(position) 
            
            maximal = max(t1, t2)
            minimal = min(t1, t2)
            
            if minimal >= 0:
                return camera + v*minimal
            elif maximal >= 0:
                return camera + v*maximal
            else:
                return False

class Color:
    def __init__(self, r, g, b):
        self.r, self.g, self.b = r, g, b
        
    def __add__(self, other):
        return Color(self.r + other.r, self.g + other.g, self.b + other.b)
    
    def to_tuple(self):
        return (int(self.r), int(self.g), int(self.b))
        
    def intensity(self, k):
        hls_color = colorsys.rgb_to_hls(self.r, self.g, self.g)
        (r, g, b) = colorsys.hls_to_rgb(hls_color[0], k*hls_color[1], hls_color[2])
        
        return Color(r, g, b)
        
    def intensity_dif(self, kr, kg, kb):
        hls_color = [0, 0, 0]
        hls_color[0] = colorsys.rgb_to_hls(self.r, 0, 0)
        hls_color[1] = colorsys.rgb_to_hls(0, self.g, 0)
        hls_color[2] = colorsys.rgb_to_hls(0, 0, self.b)     
        
        color = [0, 0, 0]
        color[0] = colorsys.hls_to_rgb(hls_color[0][0], kr*hls_color[0][1], hls_color[0][2])
        color[1] = colorsys.hls_to_rgb(hls_color[1][0], kg*hls_color[1][1], hls_color[1][2])
        color[2] = colorsys.hls_to_rgb(hls_color[2][0], kb*hls_color[2][1], hls_color[2][2])
        
        return Color(color[0][0], color[0][1], color[0][2]) + Color(color[1][0], color[1][1], color[1][2]) + Color(color[2][0], color[2][1], color[2][2])

#Image start
width = 400
height = 400
image = Image.new(mode="RGB", size=(width, height), color=(0,0,0))

#Setup parameters
camera = Vector(0, 0, -500)
plan = Vector(-width/2, -height/2, -200)
light_source = Vector(1000, 0, 0)

#Color parametes
light_color = Color(255, 255, 255)
light_intensity = 1
env_color = Color(255, 255, 255)
env_intensity = 0.1

#Constant terms
n = 20
kd = (1, 1, 1)
ks = 1

#Enviroment
env = env_color.intensity(env_intensity)
env = env.intensity_dif(kd[0], kd[1], kd[2])

sphere = Sphere(200)
for i in range(width):
    for j in range(height):
        pixel = plan + Vector(i, j, 0)
        spectral = (0, 0, 0)
        
        point = sphere.ray_casting(pixel)
        if point != False:
            #Calculate parameters
            normal = sphere.normal_vector(point)
            light = (light_source - point).normalize()
            reflection = normal*(2*light.dot_product(normal)) - light
            visual = (camera - point).normalize()
            
            #Calculate dot products
            N_L = normal.dot_product(light)
            V_R = visual.dot_product(reflection)
            
            #Phong components
            difuse = (kd[0]*N_L, kd[1]*N_L, kd[2]*N_L)
            specular = ks*(visual.dot_product(reflection) ** n)
            
            if (N_L >= 0):
                if (V_R >= 0):
                    spectral = (difuse[0] + specular, difuse[1] + specular, difuse[2] + specular)
                else:
                    spectral = difuse
        
            spectral = (light_color.intensity_dif(spectral[0], spectral[1], spectral[2]) + env).to_tuple()
            
        image.putpixel([i, j], spectral)

image.show()