progress

exam-2/exam2.md

@@ -15,7 +15,7 @@ author: |
 \newcommand{\now}[1]{\textcolor{blue}{#1}}
 \newcommand{\todo}[0]{\textcolor{red}{\textbf{TODO}}}
 
-[ 2 3 8 9 ]
+[ 3 8 9 ]
 
 ## Reflection and Refraction
 
@@ -95,11 +95,11 @@ author: |
    $\frac{v_1}{\|v_1\|_2}$ as our unit vector for $v_5$.
 
    \begin{align}
-   v_4 &= (s - p) + \tan \theta_t \times \|s - p\|_2 \times \frac{v_1}{\|v_1\|_2} \\
+   v_4 &= (s - p) + \tan \theta_t \times \|s - p\|\_2 \times \frac{v_1}{\|v_1\|\_2} \\
-   v_4 &= (\r{ (2, 2, 10) - (1, 4, 8) }) + \tan \theta_t \times \|s - p\|_2 \times \frac{v_1}{\|v_1\|_2} \\
+   v_4 &= (\r{ (2, 2, 10) - (1, 4, 8) }) + \tan \theta_t \times \|s - p\|\_2 \times \frac{v_1}{\|v_1\|\_2} \\
-   v_4 &= \r{ (1, -2, 2) } + \tan \theta_t \times \|\r{(1, -2, 2)}\|_2 \times \frac{v_1}{\|v_1\|_2} \\
+   v_4 &= \r{ (1, -2, 2) } + \tan \theta_t \times \|\r{(1, -2, 2)}\|\_2 \times \frac{v_1}{\|v_1\|\_2} \\
-   v_4 &= (1, -2, 2) + \tan \theta_t \times \r{3} \times \frac{v_1}{\|v_1\|_2} \\
+   v_4 &= (1, -2, 2) + \tan \theta_t \times \r{3} \times \frac{v_1}{\|v_1\|\_2} \\
-   v_4 &= (1, -2, 2) + \tan \theta_t \times 3 \times \frac{\r{(0, 6, 6)}}{\|\r{(0, 6, 6)}\|_2} \\
+   v_4 &= (1, -2, 2) + \tan \theta_t \times 3 \times \frac{\r{(0, 6, 6)}}{\|\r{(0, 6, 6)}\|\_2} \\
    v_4 &= (1, -2, 2) + \tan \theta_t \times 3 \times \r{\left(0, \frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}\right)} \\
    v_4 &= (1, -2, 2) + \tan \theta_t \times \left(0, \r{\frac{3}{2}\sqrt{2}}, \r{\frac{3}{2}\sqrt{2}} \right) \\
    v_4 &= (1, -2, 2) + \r{\frac{4}{7}\sqrt{2}} \times \left(0, \frac{3}{2}\sqrt{2}, \frac{3}{2}\sqrt{2} \right) \\
@@ -120,46 +120,59 @@ author: |
    at the location $p = (4, 4, 7)$, with a direction of flight $w = (2, 1, -2)$
    and the wings aligned with the direction $d = (-2, 2, -1)$.}
 
-   The order we want is (1) do all the rotations, and then (2) translate to the
+   There's 2 discrete transformations going on here:
-   spot we want. The rotation is done in multiple steps:
 
-   - First we want to make sure the nose of the plane points in the correct
+   - First, we must rotate the plane. Since we are given orthogonal directions
-     direction in the $xz$ plane (rotating around the $y$ axis). The desired
+     for the flight and wings, we can just come up with another vector for the
-     resulting direction is $(2, y, -2)$, so that means for a nose currently facing
+     tail direction by doing:
-     the $+z$ direction which is $(0, y, 1)$, we want to rotate around by around
+
-     $-\frac{3}{4}\pi$. The transformation matrix is:
+     $y' = (2, 1, -2) \times (-2, 2, -1) = (1, 2, 2)$
+
+     Now we can construct a 3D rotation matrix:
 
      $$
-     M_1 = \begin{bmatrix}
+     R =
-        1 & 0 & 0 & 0 \\
+     \begin{bmatrix}
-        0 & 1 & 0 & 0 \\
+       d_x & y_x & w_x & 0 \\
-        0 & 0 & 1 & 0 \\
+       d_y & y_y & w_y & 0 \\
-        0 & 0 & 0 & 1 \\
+       d_z & y_z & w_z & 0 \\
+       0 & 0 & 0 & 1
+     \end{bmatrix}
+     =
+     \begin{bmatrix}
+       -2 & 1 & 2 & 0 \\
+       2 & 2 & 1 & 0 \\
+       -1 & 2 & -2 & 0 \\
+       0 & 0 & 0 & 1
      \end{bmatrix}
      $$
 
-   - Then we want to rotate the plane vertically, so it's pointing in the right
+   - Now we just need to translate this to the position (4, 4, 7). This is easy
-     direction.
+     with:
+
+     $$
+     T =
+     \begin{bmatrix}
+       1 & 0 & 0 & 4 \\
+       0 & 1 & 0 & 4 \\
+       0 & 0 & 1 & 7 \\
+       0 & 0 & 0 & 1 \\
+     \end{bmatrix}
+     $$
+
+   We can just compose these two matrices together (by doing the rotation first,
+   then translating!)
 
    $$
+   TR =
    \begin{bmatrix}
-     1 & 0 & 0 & x \\
+     -2 & 1 & 2 & 4 \\
-     0 & 1 & 0 & y \\
+     2 & 2 & 1 & 4 \\
-     0 & 0 & 1 & z \\
+     -1 & 2 & -2 & 7 \\
-     0 & 0 & 0 & 1 \\
+     0 & 0 & 0 & 1
-   \end{bmatrix}
-   =
-   \begin{bmatrix}
-     1 & 0 & 0 & 4 \\
-     0 & 1 & 0 & 4 \\
-     0 & 0 & 1 & 7 \\
-     0 & 0 & 0 & 1 \\
    \end{bmatrix}
    $$
 
-   Since the direction of flight was originally $(0, 0, 1)$, we have to
-   transform it to $(2, 1, -2)$.
-
 3. \c{Consider the earth model shown below, which is defined in object
    coordinates with its center at $(0, 0, 0)$, the vertical axis through the
    north pole aligned with the direction $(0, 1, 0)$, and a horizontal plane
@@ -680,70 +693,222 @@ author: |
 
    a. \c{(2 points) What are the entries in $P$?}
 
-   The left / right values are found by using the tangent of the field-of-view
+   In this case:
-   triangle: $\tan(60^\circ) = \frac{\textrm{right}}{0.5}$, so $\textrm{right} =
+
-   \tan(60^\circ) \times 0.5 = \boxed{\frac{\sqrt{3}}{2}}$. The same goes for the
+   - near = 0.5
-   vertical, which also yields $\frac{\sqrt{3}}{2}$.
+   - far = 20
+   - left = bottom = $-\tan 30^\circ \times 0.5 = -\frac{\sqrt{3}}{6}$
+   - right = top = $\tan 30^\circ \times 0.5 = \frac{\sqrt{3}}{6}$
+   - right - left = top - bottom = $\frac{\sqrt{3}}{3} = \frac{1}{\sqrt{3}}$
 
    $$
    \begin{bmatrix}
-     \frac{2\times near}{right - left} & 0 & \frac{right + left}{right - left} & 0 \\
+     \frac{2\cdot near}{right-left} & 0 & \frac{right+left}{right-left} & 0 \\
-     0 & \frac{2\times near}{top - bottom} & \frac{top + bottom}{top - bottom} & 0 \\
+     0 & \frac{2\cdot near}{top-bottom} & \frac{top+bottom}{top-bottom} & 0 \\
-     0 & 0 & -\frac{far + near}{far - near} & -\frac{2\times far\times near}{far - near} \\
+     0 & 0 & -\frac{far+near}{far-near} & -\frac{2\cdot far\cdot near}{far-near} \\
-     0 & 0 & -1 & 0
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   =
+   \begin{bmatrix}
+     \sqrt{3} & 0 & 0 & 0 \\
+     0 & \sqrt{3} & 0 & 0 \\
+     0 & 0 & -\frac{20.5}{19.5} & \frac{-20}{19.5} \\
+     0 & 0 & -1 & 0 \\
    \end{bmatrix}
    $$
 
-   $$
-   = \begin{bmatrix}
-     \frac{2\times 0.5}{\frac{\sqrt{3}}{2} - (-\frac{\sqrt{3}}{2})} & 0 & \frac{\frac{\sqrt{3}}{2} + (-\frac{\sqrt{3}}{2})}{\frac{\sqrt{3}}{2} - (-\frac{\sqrt{3}}{2})} & 0 \\
-     0 & \frac{2\times 0.5}{\frac{\sqrt{3}}{2} - (-\frac{\sqrt{3}}{2})} & \frac{\frac{\sqrt{3}}{2} + (-\frac{\sqrt{3}}{2})}{\frac{\sqrt{3}}{2} - (-\frac{\sqrt{3}}{2})} & 0 \\
-     0 & 0 & -\frac{20 + 0.5}{20 - 0.5} & -\frac{2\times 20\times 0.5}{20 - 0.5} \\
-     0 & 0 & -1 & 0
-   \end{bmatrix}
-   $$
-
-   $$
-   = \boxed{\begin{bmatrix}
-     \frac{1}{\sqrt{3}} & 0 & 0 & 0 \\
-     0 & \frac{1}{\sqrt{3}} & 0 & 0 \\
-     0 & 0 & -\frac{41}{39} & -\frac{40}{39} \\
-     0 & 0 & -1 & 0
-   \end{bmatrix}}
-   $$
-
-   \todo the numbers are wrong lmao
-
    b. \c{(3 points) How should be matrix $P$ be re-defined if the viewing window
    is re-sized to be twice as tall as it is wide?}
 
+   Only the second row changes, since it is the only thing that references
+   bottom or top. top + bottom is still 0, so the non-diagonal cell doesn't
+   change, so only top - bottom gets doubled. The result is:
+
+   $$
+   \begin{bmatrix}
+     \sqrt{3} & 0 & 0 & 0 \\
+     0 & \r{2}\sqrt{3} & 0 & 0 \\
+     0 & 0 & -\frac{20.5}{19.5} & \frac{-20}{19.5} \\
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   $$
+
    c. \c{(3 points) What are the new horizontal and vertical fields of view
    after this change has been made?}
 
+   The horizontal one doesn't change, so only the vertical one does. Since the
+   height has increased relative to the near value, the FOV increases to
+   $120^\circ$.
+
+   \c{When the viewing frustum is known to be symmetric, we will have $left =
+   -right$ and $bottom = -top$. In that case, an alternative definition can be
+   used for the perspective projection matrix where instead of defining
+   parameters left, right, top, bottom, the programmer instead specifies a
+   vertical field of view angle and the aspect ratio of the viewing frustum.}
+
+   d. \c{(1 point) What are the entries in $P_{alt}$ when the viewing frustum is
+   defined by: a near clipping plane located 0.5 units in front of the camera, a
+   far clipping plane located 20 units from the front of the camera, a
+   $60^\circ$ vertical field of view, and a square aspect ratio?}
+
+   Using a 1:1 aspect ratio (since the field of views are the same)
+
+   - $\cot(30^\circ) = \sqrt{3}$
+
+   $$
+   \begin{bmatrix}
+     \cot \left( \frac{\theta_v}{2} \right) & 0 & 0 & 0 \\
+     0 & \cot \left( \frac{\theta_v}{2} \right) & 0 & 0 \\
+     0 & 0 & -\frac{far + near}{far - near} & \frac{-2 \cdot far \cdot near}{far - near} \\
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   =
+   \begin{bmatrix}
+     \sqrt{3} & 0 & 0 & 0 \\
+     0 & \sqrt{3} & 0 & 0 \\
+     0 & 0 & -\frac{20.5}{19.5} & \frac{-20}{19.5} \\
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   $$
+
+   e. \c{(1 points) Suppose the viewing window is re-sized to be twice as wide as
+   it is tall. How might you re-define the entries in $P_{alt}$?}
+
+   This means the aspect ratio is $\frac{1}{2}$, and the aspect ratio is applied
+   to the first entry in the matrix.
+
+   $$
+   \begin{bmatrix}
+     \r{2}\sqrt{3} & 0 & 0 & 0 \\
+     0 & \sqrt{3} & 0 & 0 \\
+     0 & 0 & -\frac{20.5}{19.5} & \frac{-20}{19.5} \\
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   $$
+
+   f. \c{(2 points) What would the new horizontal and vertical fields of view be
+   after this change has been made? How would the image contents differ from
+   when the window was square?}
+
+   The horizontal FOV has doubled, so it goes to $120^\circ$, but the vertical
+   doesn't change.
+
+   g. \c{(1 points) Suppose the viewing window is re-sized to be twice as tall
+   as it is wide. How might you re-define the entries in $P_{alt}$?}
+
+   $$
+   \begin{bmatrix}
+     \r{\frac{1}{2}}\sqrt{3} & 0 & 0 & 0 \\
+     0 & \sqrt{3} & 0 & 0 \\
+     0 & 0 & -\frac{20.5}{19.5} & \frac{-20}{19.5} \\
+     0 & 0 & -1 & 0 \\
+   \end{bmatrix}
+   $$
+
+   So in this case, the $\theta_v$ would be changed, and then we scale the
+   $\theta_h$ part using the aspect ratio so it stays consistent.
+
+   h. \c{(2 points) What would the new horizontal and vertical fields of view be
+   after this change has been made? How would the image contents differ from
+   when the window was square?}
+
+   Since the horizontal field of view hasn't changed, it remains $60^\circ$. But
+   the vertical field of view has doubled, which makes it $120^\circ$. More of
+   the vertical view would be available than if it was just in a square.
+
+   i. \c{(1 points) Suppose you wanted the user to be able to see more of the
+   scene in the vertical direction as the window is made taller. How would you
+   need to adjust $P_{alt}$ to achieve that result?}
+
+   As long as the FOV is increased, more of the scene is available. If for
+   example, when we changed the height of the window above, the _horizontal_
+   field of view changed, then we would've _reduced_ the amount of visibility in
+   the scene.
+
 ## Clipping
 
 9. \c{Consider the triangle whose vertex positions, after the viewport
    transformation, lie in the centers of the pixels: $p_0 = (3, 3), p_1 = (9,
    5), p_2 = (11, 11)$.}
 
-   Starting at $p_0$, the three vectors are:
-
-   - $v_0 = p_1 - p_0 = (9 - 3, 5 - 3) = (6, 2)$
-   - $v_1 = p_2 - p_1 = (11 - 9, 11 - 5) = (2, 6)$
-   - $v_2 = p_0 - p_2 = (3 - 11, 3 - 11) = (-8, -8)$
-
-   The first edge vector $e$ would be $(6, 2)$, and the edge normal would be
-   that rotated by $90^\circ$.
-
    a. \c{(6 points) Define the edge equations and tests that would be applied,
    during the rasterization process, to each pixel $(x, y)$ within the bounding
    rectangle $3 \le x \le 11, 3 \le y \le 11$ to determine if that pixel is
    inside the triangle or not.}
 
+   So for each edge, we have to test if the point is on the inside half of the
+   line that divides the plane. Here are the (a, b, c) values for each of the
+   edges:
+
+   - Edge 1: (-2, 6, -12)
+   - Edge 2: (-6, 2, 44)
+   - Edge 3: (8, -8, 0)
+
+   We just shove these numbers into $e(x, y) = ax + by + c$ for each point to
+   determine if it lies inside the triangle or not. I've written this Python
+   script to do the detection:
+
+   ```py
+   for (i, (p0_, p1_)) in enumerate([(p0, p1), (p1, p2), (p2, p0)]):
+     a= -(p1_[1] - p0_[1])
+     b = (p1_[0] - p0_[0])
+     c = (p1_[1] - p0_[1]) * p0_[0] - (p1_[0] - p0_[0]) * p0_[1]
+
+     for x, y in itertools.product(range(3, 12), range(3, 12)):
+       if (x, y) not in statuses: statuses[x, y] = [None, None, None]
+       e = a * x + b * y + c
+       statuses[x, y][i] = e >= 0
+   ```
+
+   I then plotted the various numbers to see if they match:
+
+   ![](9a.jpg){width=40%}
+
+   The 3 digit number corresponds to 1 if it's "inside" and 0 if it's not
+   "inside" for each of the 3 edges. The first digit corresponds to the top
+   horizontal edge, the second digit corresponds to the right most edge, and the
+   last digit corresponds to the long diagonal. When all three are 1, the pixel
+   is officially "inside" the triangle for sure.
+
+   There is also edge detection to see if the edge pixels belong to the left or
+   the top edges. I didn't implement that here but I talk about it below in the
+   second part b.
+
    b. \c{(3 points) Consider the three pixels $p_4 = (6, 4), p_5 = (7, 7)$, and
    $p_6 = (10, 8)$. Which of these would be considered to lie inside the
    triangle, according to the methods taught in class?}
 
+   For these three pixels, we can start with $p_4$ and define $a$ and $b$ using
+   it (going to $p_6$ first to remain in counter-clockwise order).
+
+   Then we use the checks to determine if the $a$ and $b$ values satisfy the
+   conditions for being left or top edges:
+
+   ```py
+   p4 = (6, 4)
+   p5 = (7, 7)
+   p6 = (10, 8)
+
+   for (i, (p0_, p1_)) in enumerate([(p4, p6), (p6, p5), (p5, p4)]):
+     a= -(p1_[1] - p0_[1])
+     b = (p1_[0] - p0_[0])
+     c = (p1_[1] - p0_[1]) * p0_[0] - (p1_[0] - p0_[0]) * p0_[1]
+
+     print(a, b, c, end=" ")
+
+     if a == 0 and b < 0: print("top")
+     elif a > 0: print("left")
+     else: print()
+   ```
+
+   This tells us that the $p_6 \rightarrow p_5$ and the $p_5 \rightarrow p_4$
+   edges are both left edges. If you graph this on the grid, this is accurate.
+   This means for those particular edges, the points that lie exactly on the
+   edge will be considered "inside" and for others, it will not.
+
+   Edge detection can be done by subtracting the point from the normal and
+   seeing if the resulting vector is normal or not.
+
 10. \c{When a model contains many triangles that form a smoothly curving surface
     patch, it can be inefficient to separately represent each triangle in the
     patch independently as a set of three vertices because memory is wasted when

exam-2/exam2.py

@@ -1,14 +1,37 @@
 import itertools
 import numpy as np
 import math
-from sympy import N, Number, Rational, init_printing, latex, simplify, Expr
+from sympy import N, Matrix, Number, Rational, init_printing, latex, simplify, Expr
 from sympy.vector import CoordSys3D, Vector
+from PIL import Image, ImageDraw
 
 init_printing()
 
 import sympy
 
+C = CoordSys3D('C')
 unit = lambda v: v/np.linalg.norm(v)
+vec = lambda a, b, c: a * C.i + b * C.j + c * C.k
+
+def ap(matrix, vector):
+  vector_ = np.r_[vector, [1]]
+  trans_ = matrix @ vector_
+  trans = trans_[:3]
+  return trans
+
+def ap2(matrix, vector):
+  c = vector.components
+  vector_ = vector.to_matrix(C).col_join(Matrix([[1]]))
+  trans_ = matrix @ vector_
+  return trans_[0] * C.i + trans_[1] * C.j + trans_[2] * C.k
+
+def pv(vector):
+  c = vector.components
+  x = c.get(C.i, 0)
+  y = c.get(C.j, 0)
+  z = c.get(C.k, 0)
+  return (x, y, z)
+
 
 def perspective_matrix(vfov, width, height, left, right, bottom, top, near, far):
   aspect = width / height
@@ -62,7 +85,6 @@ def print_bmatrix(arr):
     print("\\\\")
 
 def problem_1():
-  C = CoordSys3D('C')
   p = 1 * C.i + 4 * C.j + 8 * C.k
   e = 0 * C.i + 0 * C.j + 0 * C.k
   s = 2 * C.i + 2 * C.j + 10 * C.k
@@ -302,10 +324,9 @@ def problem_6():
     near = min(map(lambda p: p[2], points))
     far = max(map(lambda p: p[2], points))
 
+    M_this = M(left, right, bottom, top, near, far)
     for point in points:
-      point_ = np.r_[point, [1]]
+      trans = ap(M_this, point)
-      trans_ = M(left, right, bottom, top, near, far) @ point_
-      trans = trans_[:3]
       v = np.vectorize(lambda x: float(x.evalf()))
       l = np.vectorize(lambda x: latex(simplify(x)))
       point = l(point)
@@ -330,12 +351,139 @@ def problem_6():
   calculate(oblique_transform, points)
 
 def problem_9():
+  p0 = (3, 3)
+  p1 = (9, 5)
+  p2 = (11, 11)
+
+  statuses = {}
+
+  for (i, (p0_, p1_)) in enumerate([(p0, p1), (p1, p2), (p2, p0)]):
+    a= -(p1_[1] - p0_[1])
+    b = (p1_[0] - p0_[0])
+    c = (p1_[1] - p0_[1]) * p0_[0] - (p1_[0] - p0_[0]) * p0_[1]
+
+    for x, y in itertools.product(range(3, 12), range(3, 12)):
+      if (x, y) not in statuses: statuses[x, y] = [None, None, None]
+      e = a * x + b * y + c
+      statuses[x, y][i] = e >= 0
+
+  CELL_SIZE = 30
+  im = Image.new("RGB", (9 * CELL_SIZE, 9 * CELL_SIZE))
+  draw = ImageDraw.Draw(im)
+  in_color = (180, 255, 180)
+  out_color = (255, 180, 180)
+  for (x, y), status in statuses.items():
+    color = in_color if all(status) else out_color
+    sx, sy = x - 3, y - 3
+    ex, ey = sx + 1, sy + 1
+    draw.rectangle([
+      (sx * CELL_SIZE, sy * CELL_SIZE),
+      (ex * CELL_SIZE, ey * CELL_SIZE),
+    ], color)
+    text = "".join(map(lambda s: "1" if s else "0", status))
+    _, _, w, h = draw.textbbox((0, 0), text)
+    draw.text(
+      (sx * CELL_SIZE + (CELL_SIZE - w) / 2.0, sy * CELL_SIZE + (CELL_SIZE - h) / 2.0),
+      text,
+      "black"
+    )
+
+  im.save("9a.jpg")
+
+  p4 = (6, 4)
+  p5 = (7, 7)
+  p6 = (10, 8)
+
+  for (i, (p0_, p1_)) in enumerate([(p4, p6), (p6, p5), (p5, p4)]):
+    a= -(p1_[1] - p0_[1])
+    b = (p1_[0] - p0_[0])
+    c = (p1_[1] - p0_[1]) * p0_[0] - (p1_[0] - p0_[0]) * p0_[1]
+
+    print(a, b, c, end=" ")
+
+    if a == 0 and b < 0: print("top")
+    elif a > 0: print("left")
+    else: print()
+
   pass
 
-print("\nPROBLEM 8 -------------------------"); problem_8()
+def problem_2():
-print("\nPROBLEM 5 -------------------------"); problem_5()
+  y_axis_angle = 3 * sympy.pi / 4
+
+  cos_t = -2
+  sin_t = 2
+  step_1 = Matrix([
+    [cos_t, 0, sin_t, 0],
+    [0, 1, 0, 0],
+    [-sin_t, 0, cos_t, 0],
+    [0, 0, 0, 1],
+  ])
+
+  sqrt2 = sympy.sqrt(2)
+  cos_t = 2 * sqrt2
+  sin_t = 1
+  step_2 = Matrix([
+    [1, 0, 0, 0],
+    [0, cos_t, -sin_t, 0],
+    [0, sin_t, cos_t, 0],
+    [0, 0, 0, 1],
+  ])
+
+  up_dir = vec(0, 1, 0)
+  nose_dir = vec(0, 0, 1)
+  left_wing_dir = vec(1, 0, 0)
+  right_wing_dir = vec(+1, 0, 0)
+  
+  def apply(m):
+    print("- nose (z):", pv(ap2(m, nose_dir)))
+    print("- up (y):", pv(ap2(m, up_dir)))
+    print("- leftwing (+x):", pv(ap2(m, left_wing_dir)))
+    print("- rightwing (-x):", pv(ap2(m, right_wing_dir)))
+
+  print("step 1")
+  apply(step_1)
+  print()
+
+  print("step 2")
+  apply(step_2)
+  print()
+
+  print("step 2 @ step 1")
+  apply(step_2 @ step_1)
+  print()
+
+  print("step 1 @ step 2")
+  apply(step_1 @ step_2)
+  print()
+
+  print("SHIET")
+  print(vec(2, 1, -2).cross(vec(-2, 2, -1)))
+  R = Matrix([
+    [-2, 2, -1, 0],
+    [1, 2, 2, 0],
+    [2, 1, -2, 0],
+    [0, 0, 0, 1],
+  ]).transpose()
+  print(R)
+  apply(R)
+  print()
+
+  T = Matrix([
+    [1, 0, 0, 4],
+    [0, 1, 0, 4],
+    [0, 0, 1, 7],
+    [0, 0, 0, 1],
+  ])
+
+  print("T @ R")
+  print(T @ R)
+
+# print("\nPROBLEM 8 -------------------------"); problem_8()
+# print("\nPROBLEM 5 -------------------------"); problem_5()
+# print("\nPROBLEM 9 -------------------------"); problem_9()
+# print("\nPROBLEM 7 -------------------------"); problem_7()
+# print("\nPROBLEM 4 -------------------------"); problem_4()
+# print("\nPROBLEM 6 -------------------------"); problem_6()
+# print("\nPROBLEM 1 -------------------------"); problem_1()
+print("\nPROBLEM 2 -------------------------"); problem_2()
 print("\nPROBLEM 9 -------------------------"); problem_9()
-print("\nPROBLEM 7 -------------------------"); problem_7()
-print("\nPROBLEM 4 -------------------------"); problem_4()
-print("\nPROBLEM 6 -------------------------"); problem_6()
-print("\nPROBLEM 1 -------------------------"); problem_1()

flake.nix

@@ -43,7 +43,8 @@
             zip
             zathura
 
-            (python310.withPackages (p: with p; [ ipython numpy scipy sympy ]))
+            (python310.withPackages
+              (p: with p; [ ipython numpy scipy sympy pillow ]))
           ]) ++ (with toolchain; [
             cargo
             rustc