makeuoft.py

# -*- coding: utf-8 -*-
"""MakeUofT.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1SxMdm3sDstckEw5Ys8vUI2ZFxM6NqdUQ

### Import Modules
"""

!pip install open3d
!pip install trimesh

import numpy as np
import math
import open3d as o3d
import trimesh

"""### Useless Stuff"""

# Define the vertices of the cube
vertices = [
    (-0.5, -0.5, 0.5),
    (0.5, -0.5, 0.5),
    (0.5, 0.5, 0.5),
    (-0.5, 0.5, 0.5),
    (-0.5, -0.5, -0.5),
    (0.5, -0.5, -0.5),
    (0.5, 0.5, -0.5),
    (-0.5, 0.5, -0.5)
]

# Define the faces of the cube
faces = [
    (1, 2, 3, 4),
    (5, 8, 7, 6),
    (1, 5, 6, 2),
    (2, 6, 7, 3),
    (3, 7, 8, 4),
    (5, 1, 4, 8)
]

# Open the .obj file for writing
with open('cube.obj', 'w') as f:
    # Write the vertices to the file
    for v in vertices:
        f.write('v {} {} {}\n'.format(v[0], v[1], v[2]))

    # Write the faces to the file
    for face in faces:
        f.write('f {} {} {} {}\n'.format(face[0], face[1], face[2], face[3]))

mesh = o3d.io.read_triangle_mesh("cube.obj")
# print(np.asarray(mesh.triangles))

def generateV(phit, normal, cos_theta_i, startdB, distanceFromOrigin, num_samples):
    """
    Generates an array of rays using a cosine-weighted hemisphere sampling technique
    """
    V = []
    for i in range(num_samples):
        r = sampleHemisphere(normal)
        if r.dot(normal) < 0:
            r = -r
        # Calculate the reflection of the sample vector around the incident vector
        sample_reflection = r - 2 * r.dot(normal) * normal
        # Calculate the cosine factor for the specular term
        cos_theta_r = sample_reflection.dot(-cos_theta_i)
        # Calculate the dB level for the ray
        ray_dB = startdB * (cos_theta_i + cos_theta_r) / (4 * math.pi * distanceFromOrigin**2)
        # Create a new ray and add it to the list
        ray = Ray(phit, r, ray_dB)
        V.append(ray)
    return V

"""### Class Definitions

### Helper Functions
"""

tolerance=math.cos(10 * math.pi / 180.0)

def line_plane_intersection(ray: Ray, triangle: Triangle):
    line_direction = normalize(ray.direction)
    plane_point, plane_normal = triangle.vertices[0], normalize(triangle.normal)
    
    dot_product = np.dot(line_direction, plane_normal)
    if abs(dot_product) < tolerance:
        return None
    
    unit_vector = normalize(plane_point - ray.origin)
    distance = np.dot(unit_vector, plane_normal) / dot_product
    intersection_point = ray.origin + distance * line_direction
    
    return intersection_point

origin = np.array([0, 0, 0])
line_direction = np.array([1, 1, 1])
plane_point, plane_normal = np.array([2, 0, 0]), np.array([3, 4, 5])
# plane_point, plane_normal = np.array([2, 0, 0]), np.array([-3, -4, -5])

dot_product = np.dot(normalize(line_direction), normalize(plane_normal))
print(math.acos(dot_product) * 180.0 / math.pi)
if abs(dot_product) < tolerance: print("Whoops")
else:
    unit_vector = np.array(plane_point - origin)
    distance = np.dot(unit_vector, plane_normal) / dot_product
    intersection_point = origin + distance * line_direction

print(intersection_point)

11.536959032815469 + 168.46304096718453

def reflect_ray(incident_ray: Ray, surface_normal, phit, coefficient):
  
    v1 = normalize(incident_ray.direction)
    v2 = normalize(surface_normal)
    dot_product = np.dot(v1, v2)
    sub_term = normalize(2 * dot_product * v2)
    reflection = normalize(v1 - sub_term)
    reflectedRay = Ray(phit, reflection, (1-coefficient)*incident_ray.soundlevel)
    
    return reflectedRay

def line_cube_intersection(start, direction, cube_min, cube_max):
    tmin = -np.inf
    tmax = np.inf

    for i in range(3):
        if direction[i] != 0:
            t1 = (cube_min[i] - start[i]) / direction[i]
            t2 = (cube_max[i] - start[i]) / direction[i]

            if t1 > t2:
                t1, t2 = t2, t1

            if t1 > tmin:
                tmin = t1

            if t2 < tmax:
                tmax = t2

            if tmax < 0:
                return None

        # else: return None

    if tmin > tmax:
        return None

    intersection = start + tmin * direction
    return intersection

def normalize(v):
    magnitude = math.sqrt(v[0]**2+v[1]**2+v[2]**2)
    if magnitude == 0:
      return np.array([0, 0, 0])
    return np.array(v)/magnitude

"""### Utility Classes"""

class Room:
    def __init__(self, path):
        self.path=path
        self.volume=0
        self.vertices=list()
        self.triangles=np.array([None])

    def compute_details(self, point):
        mesh = trimesh.load(self.path)
        self.volume = mesh.volume
        triangles = np.array(mesh.triangles)
        new_triangles = np.array([None]*len(triangles))
        
        # for i, triangle in enumerate(triangles):
        #     v1, v2, v3 = triangle[0], triangle[1], triangle[2]
        #     self.vertices.append(v1-point)
        #     self.vertices.append(v2-point)
        #     self.vertices.append(v3-point)

        # self.vertices = np.array(self.vertices)

        for i, triangle in enumerate(triangles):
            new_t = [None]*3
            for j in range(3):
                new_v = triangle[j] - point
                new_t[j] = new_v
            # new_triangle = new_triangle.compute_normal()
            new_triangle = Triangle(i, new_t)
            new_triangle = new_triangle.compute_normal()
            new_triangles[i] = new_triangle
        self.triangles = new_triangles
        
        return self

class Triangle:
    def __init__(self, id, vertices):
        self.id=id
        self.vertices = np.array(vertices)
        self.normal = None

    def compute_normal(self):
        v1, v2, v3 = self.vertices
        vector1 = v2 - v1
        vector2 = v3 - v1
        normal = normalize(np.cross(vector1, vector2))
        self.normal = normal

        return self

    def is_point_on_plane(self, point):
        point_vector = point - self.vertices[0]
        unit_vector = normalize(point_vector)
        unit_normal = self.normal
        dot_product = abs(np.dot(unit_vector, unit_normal))
        angle = math.acos(dot_product)

        return angle >= math.pi * 88.5 / 180.0

class Source:
    def __init__(self, name, origin, numrays, soundlevel):
        self.name=name
        self.origin=np.array(origin)
        self.numrays=numrays
        self.soundlevel=soundlevel

    def generate_rays(self):
        points = []
        for i in range(self.numrays):
            for j in range(self.numrays):
                u = (i + 0.5) / self.numrays
                v = (j + 0.5) / self.numrays
                theta = 2 * math.pi * u
                phi = math.acos(2 * v - 1)
                x = math.sin(phi) * math.cos(theta)
                y = math.sin(phi) * math.sin(theta)
                z = math.cos(phi)
                direction = normalize(np.array([x, y, z]))
                points.append(Ray(self.origin, direction, self.soundlevel))
        return points

class Ray:
    def __init__(self, origin, direction, soundlevel):
        self.origin=origin
        self.direction=normalize(direction)
        self.soundlevel=soundlevel

class Listener:
    def __init__(self, position, unit):
        self.position=np.array(position)
        self.unit=unit
        self.min_coord = self.position - 0.5 * self.unit
        self.max_coord = self.position + 0.5 * self.unit

num_rays = 50
centroid = np.array([-0.0121977, 2.06714929, 1.4618923 ])
src = Source("", np.array([0, 0, 0]), num_rays, 100)
incident_rays = src.generate_rays()
num_reflections=15

room = Room("model.obj")
room = room.compute_details(centroid)

triangles = room.triangles
num_triangles = len(triangles)
coefficient = 0.05
#listener is with respect to 0 0 0
listener=Listener([2, 2, 2], 1)

reflectionData=np.array([None]*num_rays)
for i in range(num_rays):
    # print("For Ray", i+1, ":")
    reflections = np.array([(None, -2)]*(num_reflections+1))
    incident_ray = incident_rays[i]
    reflections[0] = (incident_ray, -1)

    for j in range(num_reflections):
        # print("Reflection", j, "For Ray", i+1)
        incoming_ray = reflections[j][0]

        reflected_ray = None
        for k in range(num_triangles):
            triangle = triangles[k]
            if triangle.is_point_on_plane(incoming_ray.origin): continue

            point_of_incidence = line_plane_intersection(incoming_ray, triangle)
            # print(triangle.is_point_on_plane(point_of_incidence))
            if point_of_incidence is None:
                # print("Oops")
                continue

            reflected_ray = reflect_ray(incoming_ray, triangle.normal, point_of_incidence, coefficient)
            reflections[j+1] = (reflected_ray, triangle.id)
            # break
        reflectionReachedListener = line_cube_intersection(incoming_ray.origin, incoming_ray.direction, listener.min_coord, listener.max_coord)
        if reflectionReachedListener is not None:
            break;

    reflectionData[i]=reflections

for i, ray in enumerate(reflectionData):
    print("Ray", i+1, ":")
    reflections = reflectionData[i]
    for j in range(len(reflections)):
        # if j <=100: continue
        reflected_ray = reflections[j]
        if reflected_ray[1]==-2: continue
        print("Reflected Ray", j,":", reflected_ray[0].direction, "From Triangle", reflected_ray[1])

#for each ray there is going to be a list of reflections
list_of_number_of_reflections = [0]*num_rays
for i, ray in enumerate(reflectionData):
    actual_reflections = list()
    for reflection in ray:
        reflected_ray = reflection[0]
        if reflected_ray is None: break
        actual_reflections.append(reflected_ray)
    list_of_number_of_reflections[i] = len(actual_reflections)

shortest_hit = min(list_of_number_of_reflections) - 1
output_amplitude = reflectionData[0][0][0].soundlevel * math.pow(0.95, shortest_hit)

print(output_amplitude)

#reverb time