import itertools import numpy as np from sympy import * import math unit = lambda v: v/np.linalg.norm(v) def perspective_matrix(vfov, width, height, left, right, bottom, top, near, far): aspect = width / height return np.array([ [1.0 / math.tan(vfov / 2.0) / aspect, 0, 0, 0], [0, 1.0 / math.tan(vfov / 2.0), 0, 0], [0, 0, -(far + near) / (far - near), -2.0 * far * near / (far - near)], [0, 0, -1, 0] ]) # return np.array([ # [2.0 * near / (right - left), 0, (right + left) / (right - left), 0], # [0, 2.0 * near / (top - bottom), (top + bottom) / (top - bottom), 0], # [0, 0, -(far + near) / (far - near), -(2.0 * far * near) / (far - near)], # [0, 0, -1, 0], # ]) def view_matrix(camera_pos, view_dir, up_dir): n = unit(-view_dir) u = unit(np.cross(up_dir, n)) v = np.cross(n, u) return np.array([ [u[0], u[1], u[2], -np.dot(camera_pos, u)], [v[0], v[1], v[2], -np.dot(camera_pos, v)], [n[0], n[1], n[2], -np.dot(camera_pos, n)], [0, 0, 0, 1], ]) def print_trans(before, after): def style_vec(v): start = "$[\\begin{matrix}" mid = str(v[0]) + " & " + str(v[1]) + " & " + str(v[2]) end = "\\end{matrix}]$" return f"{start}{mid}{end}" return style_vec(before) + " $\\rightarrow$ " + style_vec(after) def compute_view(near, vfov, hfov): width = 2.0 * near * math.tan(hfov / 2.0) height = 2.0 * near * math.tan(hfov / 2.0) left = -width / 2.0 right = width / 2.0 bottom = -height / 2.0 top = height / 2.0 return width, height, left, right, bottom, top def print_bmatrix(arr): for row in arr: for j, col in enumerate(row): end = " " if j == len(row) - 1 else " & " print(col, end=end) print("\\\\") def problem_1(): p = np.array([1, 4, 8]) e = np.array([0, 0, 0]) s = np.array([2, 2, 10]) i = p - e print("incoming", i) print("|I| =", np.linalg.norm(i)) n = s - p print("normal", n) n_norm = unit(n) print("normal_norm", n_norm) cos_theta_i = np.dot(-i, n) / (np.linalg.norm(i) * np.linalg.norm(n)) print("part a = cos^{-1} of ", cos_theta_i) print(np.arccos(cos_theta_i)) proj = n_norm * np.linalg.norm(i) * cos_theta_i print("proj", proj) p_ = p + proj print("proj point", p_) v2 = p_ - e print("v2", v2) sin_theta_i = np.sin(np.arccos(cos_theta_i)) print("sin theta_i =", sin_theta_i) theta_t = np.arcsin(1.0 / 1.5 * sin_theta_i) print("approx answer for part d", theta_t) uin = unit(n_norm) print("unit inv normal", uin) uperp = unit(v2) print("uperp", uperp) answer_1e = uin + math.tan(theta_t) * uperp - p print("1e answer", unit(answer_1e)) def problem_4(): camera_pos = np.array([2, 3, 5]) view_dir = np.array([1, -1, -1]) up_dir = np.array([0, 1, 0]) V = view_matrix(camera_pos, view_dir, up_dir) print(V) f = np.vectorize(lambda c: (c * c)) print(f(V)) def build_translation_matrix(vec): return np.array([ [1, 0, 0, vec[0]], [0, 1, 0, vec[1]], [0, 0, 1, vec[2]], [0, 0, 0, 1], ]) def problem_5(): b1 = build_translation_matrix(np.array([0, 0, 5])) theta = math.radians(-90) sin_theta = round( math.sin(theta), 5) cos_theta = round(math.cos(theta), 5) b2 = np.array([ [cos_theta, 0, sin_theta, 0], [0, 1, 0, 0], [-sin_theta, 0, cos_theta, 0], [0, 0, 0, 1], ]) b3 = build_translation_matrix(np.array([0, 0, -5])) print("b1", b1) print("b2", b2) print("b3", b3) M = b3 @ b2 @ b1 print("M", M) ex1 = np.array([1, 1, -4, 1]) print("ex1", ex1, M @ ex1) up = np.array([0, 1, 0]) n = np.array([1, 0, 0]) u = unit(np.cross(up, n)) v = np.cross(n, u) print(f"{up = }, {n = }, {u = }, {v = }") eye = np.array([5, 0, -5]) dx = -(np.dot(eye, u)) dy = -(np.dot(eye, v)) dz = -(np.dot(eye, n)) print(f"{dx = }, {dy = }, {dz = }") def problem_8(): near = 0.5 far = 20 print("part 8a") vfov = hfov = math.radians(60) width, height, left, right, bottom, top = compute_view(near, vfov, hfov) print(perspective_matrix(vfov, width, height, left, right, bottom, top, near, far)) print() def problem_7(): def solve(camera_pos, angle): angle_radians = math.radians(angle) near = 1 far = 10 view_dir = np.array([0, 0, -1]) up_dir = np.array([0, 1, 0]) width, height, left, right, bottom, top = compute_view(near, angle_radians, angle_radians) print("faces of the viewing frustum", left, right, bottom, top) P = perspective_matrix(angle_radians, width, height, left, right, bottom, top, near, far) V = view_matrix(camera_pos, view_dir, up_dir) print("P") print_bmatrix(np.around(P, 4)) print("V") print_bmatrix(np.around(V, 4)) m = P @ V points = [ np.array([0.5, 0.5, -4]), np.array([0.5, -0.5, -4]), np.array([-0.5, -0.5, -4]), np.array([-0.5, 0.5, -4]), np.array([0.5, 0.5, -5]), np.array([0.5, -0.5, -5]), np.array([-0.5, -0.5, -5]), np.array([-0.5, 0.5, -5]), np.array([1, 1.5, -7]), np.array([1, -0.5, -7]), np.array([-1, -0.5, -7]), np.array([-1, 1.5, -7]), np.array([1, 1.5, -9]), np.array([1, -0.5, -9]), np.array([-1, -0.5, -9]), np.array([-1, 1.5, -9]), ] for point in points: point_ = np.r_[point, [1]] trans = m @ point_ def style_vec(v): start = "$[\\begin{matrix}" mid = str(v[0]) + " & " + str(v[1]) + " & " + str(v[2]) end = "\\end{matrix}]$" return f"{start}{mid}{end}" print("-", style_vec(point), "$\\rightarrow$", style_vec(np.around(trans[:3], 4))) print("Part A") camera_pos = np.array([0, 0, 0]) angle = 90 solve(camera_pos, angle) print("Part B") camera_pos = np.array([0, 0, -2]) angle = 90 solve(camera_pos, angle) print("Part C") camera_pos = np.array([0, 0, 0]) angle = 45 solve(camera_pos, angle) def problem_6(): def calculate(points): left = min(map(lambda p: p[0], points)) right = max(map(lambda p: p[0], points)) bottom = min(map(lambda p: p[1], points)) top = max(map(lambda p: p[1], points)) near = min(map(lambda p: p[2], points)) far = max(map(lambda p: p[2], points)) step_1 = np.array([ [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1], ]) step_2 = np.array([ [1, 0, 0, -(left + right) / 2.0], [0, 1, 0, -(bottom + top) / 2.0], [0, 0, 1, -(near + far) / 2.0], [0, 0, 0, 1], ]) step_3 = np.array([ [2.0 / (right - left), 0, 0, 0], [0, 2.0 / (top - bottom), 0, 0], [0, 0, 2.0 / (far - near), 0], [0, 0, 0, 1], ]) M_ortho = step_3 @ step_2 @ step_1 for point in points: point_ = np.r_[point, [1]] trans_ = M_ortho @ point_ trans = trans_[:3] print("-", print_trans(np.around(point, 3), np.around(trans, 3))) sqrt3 = math.sqrt(3) width = 2.0 * sqrt3 cube_center = np.array([0, 0, -3.0 * sqrt3]) print("HELLOSU") points = [] for (dx, dy, dz) in itertools.product([-1, 1], [-1, 1], [-1, 1]): point = cube_center + np.array([dx * width / 2, dy * width / 2, dz * width / 2]) points.append(point) calculate(points) def problem_9(): pass print("\nPROBLEM 8 -------------------------"); problem_8() print("\nPROBLEM 1 -------------------------"); problem_1() print("\nPROBLEM 5 -------------------------"); problem_5() print("\nPROBLEM 9 -------------------------"); problem_9() print("\nPROBLEM 7 -------------------------"); problem_7() print("\nPROBLEM 6 -------------------------"); problem_6() print("\nPROBLEM 4 -------------------------"); problem_4()