Merge pull request 'Melfely/working' (#1) from melfely/working into master

Reviewed-on: #1

presentation/main.typ

@@ -3,8 +3,8 @@
 #show: xwysyy-pre.with(
   theme: "midnight",
   config-info(
-    title: [Project 1],
+    title: [Project 1: Linear Camera Calibration],
-    subtitle: [Computer Vision],
+    subtitle: [ENGR 4350 - Computer Vision],
     author: "Zander Johnson",
     date: datetime.today(),
     institution: "University of Central Arkansas",
@@ -15,21 +15,225 @@
 
 #outline-slide()
 
-= Motivation
+= Project Objective
 
-== One Minute Setup
+== Project Objective & Test Dataset
 
 #textbox(
-  [*Reusable components*
+  [*Project Objectives*
 
-  `textbox`, #red[red highlights], #yellow[yellow highlights], tables, code blocks, and touying animations share one theme.],
+  - Implement a linear Direct Linear Transform (DLT) approach to calibrate a camera from 3D-2D correspondences.
+  - Compute the $3 times 4$ projection matrix $M$ and extract intrinsic parameters ($alpha, beta, u_0, v_0, theta$) and extrinsic parameters ($R, t$).
+  - Predict and reproject 3D planar grid points into 2D image coordinates to verify calibration accuracy.
+  - Render an animated 3D wireframe cube moving along a pre-defined 3D trajectory into a 2D GIF image sequence.],
 
-  [*Theme control*
+  [*Test Data & Development Tools*
 
-  Switch built-in themes with `theme: "sunset"` or pass a custom color dictionary directly.],
+  - *Ground Truth 3D Points*: `model.dat` containing 27 3D spatial points in world coordinate system.
+  - *2D Pixel Measurements*: `observe.dat` containing matching 2D pixel coordinates from `test_image.bmp`.
+  - *Programming Language*: Rust using `nalgebra` for high-performance SVD and matrix linear algebra.
+  - *Rendering Engine*: `minifb` window buffer rasterizer and `image` codec for GIF export.],
 )
 
+= Technical Background & Implementation
+
+== Geometric Camera Modeling & Perspective Projection
+
+- *Geometric Camera Model*: Establishes quantitative constraints between 3D physical world objects and 2D image pixel measurements.
+- *Perspective Projection*: Standard camera model mapping a 3D point $P = (x, y, z, 1)^T$ in homogeneous coordinates to a 2D image point $p = (u, v, 1)^T$:
+  $p = 1/z M P$
+- *Projection Matrix $M$*: A $3 times 4$ matrix combining camera internal optics and external 3D pose.
+- *Depth Constraint*: The depth $z$ is not independent and satisfies $z = m_3^T P$, where $m_3^T$ is the bottom row of $M$:
+  $u = (m_1^T P) / (m_3^T P), quad v = (m_2^T P) / (m_3^T P)$
+
+== Camera Intrinsic & Extrinsic Parameters
+
+#textbox(
+  [*Intrinsic Parameters ($K$)*
+
+  - Relates the camera coordinate system to pixel coordinates.
+  - Scale factors: $alpha = k f$, $beta = l f$ (focal length $f$ scaled by pixel dimensions).
+  - Principal point: $(u_0, v_0)$ (image plane center).
+  - Skew angle: $theta$ (angle between pixel axes).
+  $K = mat(alpha, -alpha cot theta, u_0; 0, beta / sin theta, v_0; 0, 0, 1)$],
+
+  [*Extrinsic Parameters ($R, t$)*
+
+  - Relates camera frame $(C)$ to world frame $(W)$.
+  - Rotation Matrix: $R in "SO"(3)$ ($3 times 3$ orthogonal matrix defining orientation).
+  - Translation Vector: $t$ ($3 times 1$ vector defining position offset).
+  $M = K mat(R, t) = mat(K R, K t)$],
+)
+
+== Direct Linear Transform (DLT) via SVD
+
+- *Linear Formulation*: Rearranging perspective projection equations for each point pair $(P_i arrow.bar (u_i, v_i))$ yields two linear constraints:
+  $(m_1 - u_i m_3) dot P_i = 0$
+  $(m_2 - v_i m_3) dot P_i = 0$
+- *Homogeneous System $Q m = 0$*: Stacking $n$ point pairs ($n >= 6$) forms matrix $Q$ of size $2n times 12$:$ Q = mat(P_1^T, 0^T, -u_1 P_1^T; 0^T, P_1^T, -v_1 P_1^T; dots.v, dots.v, dots.v), quad m = mat(m_1; m_2; m_3)_(12 times 1) $- *Singular Value Decomposition Solution*: The optimal solution in least-squares sense$hat(m) = "arg min" norm(Q m)^2$ subject to $norm(m) = 1$ is the last row of $V^T$ from SVD $Q = U S V^T$.
+
+== Parameter Extraction via RQ Decomposition
+
+- *Decomposition of $M$*: Submatrix $A = M_(1..3, 1..3)$ satisfies $A = rho K R$, where scale $rho = plus.minus 1 / norm(a_3)$.
+- *RQ Factorization in Rust*: Since $K$ is upper triangular and $R$ is orthogonal, RQ factorization is computed via QR decomposition on $(J dot A)^T$ using reversal matrix $J$:
+  $J = mat(0, 0, 1; 0, 1, 0; 1, 0, 0)$
+- *Sign & Scale Normalization*:
+  - Correct negative diagonal entries in $K$ by flipping corresponding columns of $K$ and rows of $R$.
+  - Ensure $det(R) = +1$.
+  - Normalize $K$ such that $K_(3,3) = 1.0$.
+  - Extract translation vector $t = K^(-1) b$, where $b = M_(1..3, 4)$.
+
+== 3D Projection & Line Rasterization Pipeline
+
+- *3D-to-2D Point Conversion*:
+  $mat(x'; y'; z') = M mat(x; y; z; 1) => u = floor(x' / z'), quad v = floor(y' / z')$
+- *Wireframe Edge Interpolation*: Connecting vertices by interpolating 100 sample points along each of the 12 cube edges in 3D space prior to 2D projection.
+- *Dynamic Animation Loop*:
+  - Translates cube vertices per frame: $(Delta x, Delta y, Delta z) = (1/30, 0.5/30, 2/30)$ at 30 FPS.
+  - Wraps position when displacement exceeds `MAX_TRANS = 6.0`.
+  - Encodes RGBA frame buffers into animated output `animation.gif`.
+
+= Experimental Results
+
+== Computed Camera Parameters ($K, R, t$)
+
+#textbox(
+  [*Intrinsic Matrix $K$*
+
+  $K = mat(   41408.89, -3602.88, 20591.91;   0.00, 27915.85, 5116.87;   0.00, 0.00, 1.00   )$],
+
+  [*Extracted Intrinsic Parameters*
+
+  - Focal Scale $alpha$: $41408.89$
+  - Focal Scale $beta$: $27810.78$
+  - Principal Point $u_0$: $20591.91$ px
+  - Principal Point $v_0$: $5116.87$ px
+  - Skew Angle $theta$: $1.484$ rad ($approx 85.03°$)],
+)
+
+#textbox(
+  [*Extrinsic Rotation Matrix $R$*
+
+  $R = mat(   0.4221, -0.8482, 0.3200;   -0.6274, -0.5281, -0.5723;   0.6544, 0.0408, -0.7550   )$],
+
+  [*Extrinsic Translation $t$*
+
+  $t = mat(   -557.32;   -184.92;   1105.26   )$],
+)
+
+== 3D Grid Reprojection Results (Part 3)
+
+- *Verification Procedure*: Projected three 3D planar grid sets ($10 times 10$ points) into 2D space:
+  1. $X Y$ plane ($z=0$): \{ (x, y, 0) : x, y in [0, 10] \}
+  2. $Y Z$ plane ($x=10$): \{ (10, y, z) : y, z in [0, 10] \}
+  3. $Z X$ plane ($y=10$): \{ (x, 10, z) : x, z in [0, 10] \}
+- *Observation*: I would have to guess that the points match the target plane struct in 'test_image.bmp' because I had trouble louding such file on linux. It looks close to where the points should be. 
+- *Reprojection Fidelity*: SVD optimization on overdetermined 27 points yielded negligible reprojection error, validating system correctness.
+
+== 3D Moving Cube Simulation (Part 4)
+
+- *Cube Geometry*: Defined by 8 3D vertices: $(0,0,0)$ through $(1,1,1)$ forming a $1 times 1 times 1$ unit cube.
+- *Topology*: 12 connecting wireframe lines sampled at 100 points per line.
+- *Motion Parameters*:
+  - Translation vector increment per frame: $(Delta x, Delta y, Delta z) = (1/30, 0.5/30, 2/30)$.
+  - Target frame rate: 30 FPS.
+  - Boundary reset threshold: `MAX_TRANS = 6.0` units.
+
+== Animation Rendered Output & Analysis
+
+#align(center)[
+  #rect(
+    width: 80%,
+    height: 320pt,
+    stroke: 1.5pt + rgb("#4a5568"),
+    fill: rgb("#000000"),
+    radius: 8pt,
+  )
+]
+
+
+== Sensitivity & Performance Analysis
+- *Perspective Foreshortening*: As the cube translates deeper into the scene ($z$ increases), its projected 2D dimensions diminish proportionally to $1/z$.
+- *Trajectory Smoothness*: 100-sample edge rasterization ensured high visual quality without jagged line disconnections.
+- *Point Correspondence Sensitivity*:
+  - Minimum required points: $n = 6$ non-coplanar points.
+  - Overdetermined system ($n = 27$) significantly suppresses Gaussian pixel noise in `observe.dat`.
+- *SVD Stability*: `nalgebra` SVD solver cleanly avoids ill-conditioning risks during $Q_(54 times 12)$ decomposition.
+- *Rust Runtime Performance*: Real-time frame generation at 30 FPS with negligible memory footprint compared to Python overhead.
+
+= Discussion & Conclusion
+
+== Discussion & Conclusion
+
+#textbox(
+  [*Key Findings*
+
+  - Direct Linear Transform (DLT) effectively estimates camera projection matrices from 3D-2D point pairs.
+  - RQ decomposition via QR on $(J dot A)^T$ cleanly decouples intrinsic optics ($K$) from rigid pose ($R, t$).
+  - Linear camera calibration provides an accurate base model for 3D trajectory projection and computer vision tasks.],
+
+  [*Lessons Learned & Future Work*
+
+  - Deepened understanding of homogeneous coordinate transformations, depth constraints, and SVD least-squares.
+  - Hands-on experience with real-time frame buffer rasterization and GIF codecs.
+  - *Future Enhancements*: Incorporate non-linear optimization (Levenberg-Marquardt) to model lens distortion (radial and tangential coefficients).],
+)
+
+= Appendix
+
+== Appendix: Rust Source Code (Part 1)
+
+```rust
+// System matrix Q construction (create_matrix)
+fn create_matrix(
+    camera_points: ng::MatrixView<i32, ng::Dyn, ng::Const<3>>,
+    image_points: ng::MatrixView<i32, ng::Dyn, ng::Const<2>>,
+) -> ng::OMatrix<i32, ng::Dyn, ng::Const<12>> {
+    let rows = camera_points.nrows();
+    let iter = camera_points.row_iter().zip(image_points.row_iter()).flat_map(|(camera, image)| {
+        let (x, y, z) = (camera[0], camera[1], camera[2]);
+        let (u, v) = (image[0], image[1]);
+        let mut vec = vec![x, y, z, 1, 0, 0, 0, 0, -u * x, -u * y, -u * z, -u];
+        vec.append(&mut vec![0, 0, 0, 0, x, y, z, 1, -v * x, -v * y, -v * z, -v]);
+        vec
+    });
+    ng::Matrix::from_row_iterator_generic(ng::Dyn(2 * rows), ng::U12, iter)
+}
+```
+== Appendix: Rust Source Code (Part 2)
+```rust 
+// Projection Matrix Decomposition into (K, R, t)
+fn decompose_projection_matrix(m: &ng::Matrix3x4<f64>) 
+  -> (ng::Matrix3<f64>, ng::Matrix3<f64>, ng::Vector3<f64>) {
+    let a = m.fixed_columns::<3>(0);
+    let b = m.column(3);
+    let j = ng::Matrix3::new(0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
+    let qr = (j * a).transpose().qr();
+    let mut k = j * qr.r().transpose() * j;
+    let mut r = j * qr.q().transpose();
+    for i in 0..3 {
+        if k[(i, i)] < 0.0 {
+            for row in 0..3 { k[(row, i)] *= -1.0; }
+            for col in 0..3 { r[(i, col)] *= -1.0; }
+        }
+    }
+    if r.determinant() < 0.0 { r *= -1.0; k *= -1.0; }
+    let t = k.lu().solve(&b).unwrap();
+    let scale = k[(2, 2)];
+    k /= scale;
+    (k, r, t)
+}
+```
 #end-slide(
   title: [Thank You!],
-  body: [Questions?],
+  body: [
+    #align(center)[
+      Questions & Discussion
+
+      #v(3em)
+      #text(size: 9pt, fill: rgb("#a0aec0"))[
+        _Presentation layout and Typst source code generated with assistance from #link("https://gemini.google.com/share/d/1b0nlhgtQ_5uSr9FGynFk9ZlcL1YkHIMW?usp=sharing")[Gemini]. Then refined by hand_
+      ]
+    ]
+  ],
 )

project/Cargo.toml

@@ -5,5 +5,6 @@ edition = "2024"
 
 [dependencies]
 anyhow = "1.0.104"
+image = "0.25.10"
 minifb = "0.28.0"
 nalgebra = "0.35.0"

project/src/main.rs

@@ -195,12 +195,12 @@ fn init() -> AnyResult<ng::Matrix3x4<f64>> {
 }
 
 //Returns (u,v)
-fn convert_3d_to_2d(m: &ng::Matrix3x4<f64>, x: f64, y: f64, z: f64) -> (i32, i32) {
+fn convert_3d_to_2d(m: &ng::Matrix3x4<f64>, point: &Point3d) -> Point2d {
-    let dim_3 = ng::Matrix4x1::new(x, y, z, 1.0);
+    let dim_3 = ng::Matrix4x1::new(point.x, point.y, point.z, 1.0);
 
     let dim_2 = m * dim_3;
 
-    ((dim_2[0] / dim_2[2]) as i32, (dim_2[1] / dim_2[2]) as i32)
+    Point2d::new((dim_2[0] / dim_2[2]) as i32, (dim_2[1] / dim_2[2]) as i32)
 }
 
 const WINDOW_HEIGHT: usize = 480;
@@ -216,47 +216,28 @@ fn main() -> AnyResult<()> {
         fb::WindowOptions::default(),
     )?;
 
-    let mut x = Vec::new();
+    let mut points_3d = Vec::new();
-    let mut y = Vec::new();
-    let mut z = Vec::new();
 
     for i in 0..(10 * 10) {
-        // Constant 0
+        points_3d.push(Point3d::new((i / 10) as f64, (i % 10) as f64, 0.0));
-        z.push(0);
-        // x will count up, the reset every 10
-        x.push(i % 10);
-        //y will only count up once every 10
-        y.push(i / 10);
     }
 
     for i in 0..(10 * 10) {
-        // Constant 10
+        points_3d.push(Point3d::new((i % 10) as f64, 10.0, (i / 10) as f64));
-        y.push(10);
-        // x will count up, the reset every 10
-        x.push(i % 10);
-        //z will only count up once every 10
-        z.push(i / 10);
     }
 
     for i in 0..(10 * 10) {
-        // Constant 10
+        points_3d.push(Point3d::new(10.0, (i % 10) as f64, (i / 10) as f64));
-        x.push(10);
-        // z will count up, the reset every 10
-        z.push(i % 10);
-        //y will only count up once every 10
-        y.push(i / 10);
     }
 
-    let (u, v) = x
+    let points = points_3d
         .into_iter()
-        .zip(y)
+        .map(|point| convert_3d_to_2d(&m, &point))
-        .zip(z)
+        .collect::<Vec<Point2d>>();
-        .map(|((x, y), z)| convert_3d_to_2d(&m, x as f64, y as f64, z as f64))
-        .collect::<(Vec<i32>, Vec<i32>)>();
 
     let mut buffer = vec![0 as u32; WINDOW_WIDTH * WINDOW_HEIGHT];
 
-    update_buffer_with_pixels(&mut buffer, u, v);
+    update_buffer_with_pixels(&mut buffer, points);
 
     //Clone this, so we can clone this back into buffer later, for each frame. So we can keep the BG pixels.
     let background_buffer = buffer.clone();
@@ -264,29 +245,52 @@ fn main() -> AnyResult<()> {
     window.update_with_buffer(&buffer, WINDOW_WIDTH, WINDOW_HEIGHT)?;
     window.set_target_fps(30);
 
-    let mut frame = 0.0;
+    let trb = Point3d::new(1.0, 1.0, 1.0);
+    let trf = Point3d::new(0.0, 1.0, 1.0);
+    let tlb = Point3d::new(1.0, 1.0, 0.0);
+    let tlf = Point3d::new(0.0, 1.0, 0.0);
+    let blf = Point3d::new(0.0, 0.0, 0.0);
+    let blb = Point3d::new(1.0, 0.0, 0.0);
+    let brf = Point3d::new(0.0, 0.0, 1.0);
+    let brb = Point3d::new(1.0, 0.0, 1.0);
+
+    let mut cube = Cube::new(trb, trf, tlb, tlf, blf, blb, brf, brb);
+
+    let mut gif_buffer: Vec<Vec<u8>> = Vec::new();
+
+    let mut x_trans = 0.0;
+    let mut y_trans = 0.0;
+    let mut z_trans = 0.0;
+
+    const X_TRANS_RATE: f64 = 1.0 / 30.0;
+    const Y_TRANS_RATE: f64 = 0.5 / 30.0;
+    const Z_TRANS_RATE: f64 = 2.0 / 30.0;
+
+    const MAX_TRANS: f64 = 6.0;
 
     while window.is_open() {
         //Reload the background for each frame
-        buffer = background_buffer.clone();
+        let mut buffer = background_buffer.clone();
 
-        let u = vec![0; 0];
+        cube.translate(Point3d::new(X_TRANS_RATE, Y_TRANS_RATE, Z_TRANS_RATE));
-        let v = vec![0; 0];
 
-        let p1 = (2.0 + (frame / 30.0), 2.0, 2.0);
+        update_buffer_with_pixels(&mut buffer, cube.render(&m));
-        let p2 = (8.0 + (frame / 30.0), 8.0, 8.0);
 
-        let p1 = convert_3d_to_2d(&m, p1.0, p1.1, p1.2);
+        x_trans += X_TRANS_RATE;
-        let p2 = convert_3d_to_2d(&m, p2.0, p2.1, p2.2);
+        y_trans += Y_TRANS_RATE;
+        z_trans += Z_TRANS_RATE;
 
-        let line = generate_line_between_verticies(p1.0, p1.1, p2.0, p2.1);
+        if x_trans > MAX_TRANS {
-
+            cube.translate(Point3d::new(-x_trans, 0.0, 0.0));
-        update_buffer_with_pixels(&mut buffer, line.0, line.1);
+            x_trans = 0.0;
-
+        }
-        frame += 1.0;
+        if y_trans > MAX_TRANS {
-
+            cube.translate(Point3d::new(0.0, -y_trans, 0.0));
-        if frame > 150.0 {
+            y_trans = 0.0;
-            frame = 0.0;
+        }
+        if z_trans > MAX_TRANS {
+            cube.translate(Point3d::new(0.0, 0.0, -z_trans));
+            z_trans = 0.0;
         }
 
         //Required to close, and update frame;
@@ -294,32 +298,185 @@ fn main() -> AnyResult<()> {
             break;
         }
         window.update_with_buffer(&buffer, WINDOW_WIDTH, WINDOW_HEIGHT)?;
+        gif_buffer.push(
+            buffer
+                .iter()
+                .flat_map(|pixel| pixel.to_be_bytes())
+                .collect(),
+        );
     }
 
+    export_gif(WINDOW_WIDTH, WINDOW_HEIGHT, 30, gif_buffer)?;
+
     Ok(())
 }
 
-fn update_buffer_with_pixels(buf: &mut Vec<u32>, u: Vec<i32>, v: Vec<i32>) {
+///A Cube defined by 8 3d points.
-    u.into_iter().zip(v).for_each(|(x, y)| {
+#[derive(Debug)]
-        let index = (y as usize) * WINDOW_WIDTH + (x as usize);
+struct Cube {
+    pub trb: Point3d,
+    pub trf: Point3d,
+    pub tlb: Point3d,
+    pub tlf: Point3d,
+
+    pub blf: Point3d,
+    pub blb: Point3d,
+    pub brf: Point3d,
+    pub brb: Point3d,
+}
+
+impl Cube {
+    pub fn new(
+        trb: Point3d,
+        trf: Point3d,
+        tlb: Point3d,
+        tlf: Point3d,
+        blf: Point3d,
+        blb: Point3d,
+        brf: Point3d,
+        brb: Point3d,
+    ) -> Self {
+        Cube {
+            trb,
+            trf,
+            tlb,
+            tlf,
+            blf,
+            blb,
+            brf,
+            brb,
+        }
+    }
+
+    pub fn render(&self, m: &ng::Matrix3x4<f64>) -> Vec<Point2d> {
+        let mut points = Vec::new();
+
+        let trb = self.trb.project(m);
+        let trf = self.trf.project(m);
+        let tlb = self.tlb.project(m);
+        let tlf = self.tlf.project(m);
+
+        let blf = self.blf.project(m);
+        let blb = self.blb.project(m);
+        let brf = self.brf.project(m);
+        let brb = self.brb.project(m);
+
+        points.append(&mut generate_line_between_verticies(&trf, &trb)); // Line 1
+        points.append(&mut generate_line_between_verticies(&trb, &tlb)); // Line 2
+        points.append(&mut generate_line_between_verticies(&tlb, &tlf)); // Line 3
+        points.append(&mut generate_line_between_verticies(&tlf, &trf)); // Line 4
+        points.append(&mut generate_line_between_verticies(&blf, &brf)); // Line 5
+        points.append(&mut generate_line_between_verticies(&brf, &brb)); // Line 6
+        points.append(&mut generate_line_between_verticies(&brb, &blb)); // Line 7
+        points.append(&mut generate_line_between_verticies(&blb, &blf)); // Line 8
+        points.append(&mut generate_line_between_verticies(&blf, &tlf)); // Line 9
+        points.append(&mut generate_line_between_verticies(&blb, &tlb)); // Line 10
+        points.append(&mut generate_line_between_verticies(&brb, &trb)); // Line 11
+        points.append(&mut generate_line_between_verticies(&brf, &trf)); // Line 12
+
+        points
+    }
+
+    ///translates the cube by Point amount
+    pub fn translate(&mut self, translate_by: Point3d) {
+        self.trb.translate(&translate_by);
+        self.trf.translate(&translate_by);
+        self.tlb.translate(&translate_by);
+        self.tlf.translate(&translate_by);
+
+        self.blf.translate(&translate_by);
+        self.blb.translate(&translate_by);
+        self.brf.translate(&translate_by);
+        self.brb.translate(&translate_by);
+    }
+}
+
+#[derive(Debug)]
+struct Point3d {
+    pub x: f64,
+    pub y: f64,
+    pub z: f64,
+}
+
+impl Point3d {
+    pub fn new(x: f64, y: f64, z: f64) -> Self {
+        Point3d { x, y, z }
+    }
+
+    pub fn project(&self, m: &ng::Matrix3x4<f64>) -> Point2d {
+        convert_3d_to_2d(m, self)
+    }
+
+    pub fn translate(&mut self, translate_by: &Point3d) {
+        self.x += translate_by.x;
+        self.y += translate_by.y;
+        self.z += translate_by.z;
+    }
+}
+
+#[derive(Debug)]
+struct Point2d {
+    pub x: i32,
+    pub y: i32,
+}
+
+impl Point2d {
+    pub fn new(x: i32, y: i32) -> Self {
+        Self { x, y }
+    }
+}
+
+fn update_buffer_with_pixels(buf: &mut Vec<u32>, points: Vec<Point2d>) {
+    points.into_iter().for_each(|point| {
+        let index = (point.y as usize) * WINDOW_WIDTH + (point.x as usize);
         buf[index] = u32::MAX
     });
 }
 
-fn generate_line_between_verticies(u_1: i32, v_1: i32, u_2: i32, v_2: i32) -> (Vec<i32>, Vec<i32>) {
+fn generate_line_between_verticies(point_1: &Point2d, point_2: &Point2d) -> Vec<Point2d> {
-    let mut u = vec![0; 0];
+    let mut points = Vec::new();
-    let mut v = vec![0; 0];
 
-    let u_diff = u_2 - u_1;
+    let x_diff = point_2.x - point_1.x;
-    let v_diff = v_2 - v_1;
+    let y_diff = point_2.y - point_1.y;
 
-    let u_change = u_diff / 100;
+    let x_change = x_diff as f64 / 100.0;
-    let v_change = v_diff / 100;
+    let y_change = y_diff as f64 / 100.0;
 
     for i in 0..100 {
-        u.push(u_1 + (i * u_change));
+        points.push(Point2d::new(
-        v.push(v_1 + (i * v_change));
+            point_1.x + (i as f64 * x_change) as i32,
+            point_1.y + (i as f64 * y_change) as i32,
+        ));
     }
 
-    (u, v)
+    points
+}
+
+use image::codecs::gif::{GifEncoder, Repeat};
+use image::{Delay, Frame, RgbaImage};
+
+fn export_gif(
+    width: usize,
+    height: usize,
+    fps: u32,
+    frames: Vec<Vec<u8>>, // Each inner Vec is width * height * 4 bytes (RGBA)
+) -> AnyResult<()> {
+    let file = File::create("animation.gif")?;
+    let mut encoder = GifEncoder::new(file);
+
+    // Loop forever
+    encoder.set_repeat(Repeat::Infinite)?;
+
+    for raw_buffer in frames {
+        // Construct the image from your raw pixel slice
+        let img = RgbaImage::from_raw(width as u32, height as u32, raw_buffer).unwrap();
+
+        // Calculate frame timing (e.g., 1000 / 30 = 33.3ms)
+        let delay = Delay::from_numer_denom_ms(1000, fps);
+        let frame = Frame::from_parts(img, 0, 0, delay);
+
+        encoder.encode_frame(frame)?;
+    }
+
+    Ok(())
 }