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

Reviewed-on: #1

.gitattributes

@@ -1,2 +1,5 @@
 *.doc filter=lfs diff=lfs merge=lfs -text
 *.bmp filter=lfs diff=lfs merge=lfs -text
+*.odt filter=lfs diff=lfs merge=lfs -text
+presentation/main.pdf filter=lfs diff=lfs merge=lfs -text
+presentation/main.pptx filter=lfs diff=lfs merge=lfs -text

presentation/main.typ

@@ -0,0 +1,239 @@
+#import "@preview/xwysyy:0.4.0": *
+
+#show: xwysyy-pre.with(
+  theme: "midnight",
+  config-info(
+    title: [Project 1: Linear Camera Calibration],
+    subtitle: [ENGR 4350 - Computer Vision],
+    author: "Zander Johnson",
+    date: datetime.today(),
+    institution: "University of Central Arkansas",
+  ),
+)
+
+#title-slide()
+
+#outline-slide()
+
+= Project Objective
+
+== Project Objective & Test Dataset
+
+#textbox(
+  [*Project Objectives*
+
+  - 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.],
+
+  [*Test Data & Development Tools*
+
+  - *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: [
+    #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

@@ -0,0 +1,10 @@
+[package]
+name = "project"
+version = "0.1.0"
+edition = "2024"
+
+[dependencies]
+anyhow = "1.0.104"
+image = "0.25.10"
+minifb = "0.28.0"
+nalgebra = "0.35.0"

project/src/main.rs

@@ -0,0 +1,482 @@
+use anyhow::Result as AnyResult;
+use minifb as fb;
+use std::ops::Mul;
+use std::{fs::File, io::Read};
+
+use nalgebra as ng;
+
+fn get_data(file_path: &str) -> AnyResult<Vec<i32>> {
+    let mut file = File::open(file_path)?;
+
+    let mut buf = "".to_string();
+    file.read_to_string(&mut buf)?;
+
+    Ok(buf
+        .split_whitespace()
+        .filter_map(|f| f.parse::<i32>().ok())
+        .collect())
+}
+
+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>> {
+    //Verify row counts are equal
+    if camera_points.nrows() != image_points.nrows() {
+        panic!("Must have an equal number of rows between the image and camera points");
+    }
+
+    //Cache for later use
+    let rows = camera_points.nrows();
+
+    let iter = camera_points
+        .row_iter() //Iterate over each row
+        .zip(image_points.row_iter()) //Zip up the rows of image points
+        .flat_map(|(camera, image)| {
+            //Extract x,y,z from cords
+            let x = camera[0];
+            let y = camera[1];
+            let z = camera[2];
+
+            //Extract uv from image
+            let u = image[0];
+            let v = image[1];
+
+            //Create a 12 long vector which is the first row
+            let mut vec = vec![x, y, z, 1, 0, 0, 0, 0, -u * x, -u * y, -u * z, -u];
+            //Append another 12 elements to create another row
+            vec.append(&mut vec![
+                0,
+                0,
+                0,
+                0,
+                x,
+                y,
+                z,
+                1,
+                -v * x,
+                -v * y,
+                -v * z,
+                -v,
+            ]);
+            //return the vec which can be converted into an iterator
+            vec
+        });
+    ng::Matrix::from_row_iterator_generic(ng::Dyn(2 * rows), ng::U12, iter)
+}
+
+/// Returns (K,R,t) in that order
+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);
+
+    // Reversal matrix J
+    let j = ng::Matrix3::new(0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
+
+    // RQ via QR on (J * A)^T
+    let ja_t = (j * a).transpose();
+    let qr = ja_t.qr();
+    let q = qr.q();
+    let r_qr = qr.r();
+
+    let mut k = j * r_qr.transpose() * j;
+    let mut r = j * q.transpose();
+
+    // Fix negative diagonal elements in K
+    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;
+            }
+        }
+    }
+
+    // Ensure proper rotation determinant
+    if r.determinant() < 0.0 {
+        r *= -1.0;
+        k *= -1.0;
+    }
+
+    // Extract translation vector t = K_raw^-1 * b
+    let t = k.lu().solve(&b).unwrap();
+
+    // Normalize K so K[2,2] == 1.0
+    let scale = k[(2, 2)];
+    k /= scale;
+
+    (k, r, t)
+}
+
+fn init() -> AnyResult<ng::Matrix3x4<f64>> {
+    let observe_data = get_data("../observe.dat")?;
+    let observe_mat = ng::MatrixXx2::from_row_iterator_generic(
+        ng::Dyn(observe_data.len() / 2),
+        ng::Const::<2>,
+        observe_data,
+    );
+
+    println!("Image:\n{:}", observe_mat);
+
+    let model_data = get_data("../model.dat")?;
+    let model_mat = ng::MatrixXx3::from_row_iterator_generic(
+        ng::Dyn(model_data.len() / 3),
+        ng::Const::<3>,
+        model_data,
+    );
+
+    println!("Camera:\n{:}", model_mat);
+
+    let q = create_matrix(model_mat.as_view(), observe_mat.as_view());
+
+    let q = ng::OMatrix::from_iterator_generic(
+        ng::Dyn(q.nrows()),
+        ng::U12,
+        q.iter().map(|num| *num as f64),
+    );
+
+    println!("Q:\n{:}", q);
+
+    let v_t = ng::SVD::new(q, false, true).v_t.expect("Expected");
+
+    println!("V:\n{:}", v_t);
+
+    let m_flat = v_t.row(v_t.nrows() - 1);
+
+    let m = ng::OMatrix::from_row_iterator_generic(ng::U3, ng::U4, m_flat.iter().copied());
+
+    println!("M:\n{:}", m);
+
+    let model_mat = ng::OMatrix::from_iterator_generic(
+        ng::U4,
+        ng::Dyn(model_mat.nrows()),
+        model_mat
+            .row_iter()
+            .flat_map(|row| vec![row[0], row[1], row[2], 1]),
+    );
+
+    println!("P:\n{:}", model_mat);
+
+    let predicted_mat = ng::OMatrix::from_iterator_generic(
+        ng::U2,
+        ng::Dyn(model_mat.ncols()),
+        model_mat.column_iter().flat_map(|column| {
+            let new_column = &m.mul(&column.map(|x| x as f64));
+            let x = new_column[0];
+            let y = new_column[1];
+            let z = new_column[2];
+
+            vec![(x / z) as i32, (y / z) as i32]
+        }),
+    );
+
+    println!("Predict:\n{:}", predicted_mat);
+
+    let (k, r, t) = decompose_projection_matrix(&m);
+
+    println!("K:\n{:}\nR:\n{:}\nt:\n{:}", k, r, t);
+
+    let alpha = k[0];
+    let theta = 1.0_f64.atan2(k[3] / -alpha);
+    let beta = k[4] * theta.sin();
+    let u_0 = k[6];
+    let v_0 = k[7];
+
+    println!(
+        "alpha: {:}, theta: {:}, beta: {:}, u_0: {:}, v_0: {:}",
+        alpha, theta, beta, u_0, v_0
+    );
+
+    Ok(m)
+}
+
+//Returns (u,v)
+fn convert_3d_to_2d(m: &ng::Matrix3x4<f64>, point: &Point3d) -> Point2d {
+    let dim_3 = ng::Matrix4x1::new(point.x, point.y, point.z, 1.0);
+
+    let dim_2 = m * dim_3;
+
+    Point2d::new((dim_2[0] / dim_2[2]) as i32, (dim_2[1] / dim_2[2]) as i32)
+}
+
+const WINDOW_HEIGHT: usize = 480;
+const WINDOW_WIDTH: usize = 640;
+
+fn main() -> AnyResult<()> {
+    //Handles all of the startup init code! Instead of keeping it in the way of the actual image drawing.
+    let m = init()?;
+    let mut window = fb::Window::new(
+        "Display Buffer",
+        WINDOW_WIDTH,
+        WINDOW_HEIGHT,
+        fb::WindowOptions::default(),
+    )?;
+
+    let mut points_3d = Vec::new();
+
+    for i in 0..(10 * 10) {
+        points_3d.push(Point3d::new((i / 10) as f64, (i % 10) as f64, 0.0));
+    }
+
+    for i in 0..(10 * 10) {
+        points_3d.push(Point3d::new((i % 10) as f64, 10.0, (i / 10) as f64));
+    }
+
+    for i in 0..(10 * 10) {
+        points_3d.push(Point3d::new(10.0, (i % 10) as f64, (i / 10) as f64));
+    }
+
+    let points = points_3d
+        .into_iter()
+        .map(|point| convert_3d_to_2d(&m, &point))
+        .collect::<Vec<Point2d>>();
+
+    let mut buffer = vec![0 as u32; WINDOW_WIDTH * WINDOW_HEIGHT];
+
+    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();
+
+    window.update_with_buffer(&buffer, WINDOW_WIDTH, WINDOW_HEIGHT)?;
+    window.set_target_fps(30);
+
+    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
+        let mut buffer = background_buffer.clone();
+
+        cube.translate(Point3d::new(X_TRANS_RATE, Y_TRANS_RATE, Z_TRANS_RATE));
+
+        update_buffer_with_pixels(&mut buffer, cube.render(&m));
+
+        x_trans += X_TRANS_RATE;
+        y_trans += Y_TRANS_RATE;
+        z_trans += Z_TRANS_RATE;
+
+        if x_trans > MAX_TRANS {
+            cube.translate(Point3d::new(-x_trans, 0.0, 0.0));
+            x_trans = 0.0;
+        }
+        if y_trans > MAX_TRANS {
+            cube.translate(Point3d::new(0.0, -y_trans, 0.0));
+            y_trans = 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;
+        if window.is_key_down(fb::Key::Q) {
+            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(())
+}
+
+///A Cube defined by 8 3d points.
+#[derive(Debug)]
+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(point_1: &Point2d, point_2: &Point2d) -> Vec<Point2d> {
+    let mut points = Vec::new();
+
+    let x_diff = point_2.x - point_1.x;
+    let y_diff = point_2.y - point_1.y;
+
+    let x_change = x_diff as f64 / 100.0;
+    let y_change = y_diff as f64 / 100.0;
+
+    for i in 0..100 {
+        points.push(Point2d::new(
+            point_1.x + (i as f64 * x_change) as i32,
+            point_1.y + (i as f64 * y_change) as i32,
+        ));
+    }
+
+    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(())
+}

project1.md

@@ -0,0 +1,100 @@
+**ENGR4350 Computer Vision**
+
+**Project 1. Linear Approach to Camera Calibration**
+
+**Due Data: September 20, 2026**
+
+**Objective:** To use a linear approach to calibrate a camera and to use
+the camera parameters to predict the 2D image-domain location of a set
+of 3D points or a 3D moving object.
+
+Procedure:
+
+1.  Study the lecture note by showing the use of projection matrix that
+    maps objects from 3-D to 2-D, the definition of intrinsic and
+    extrinsic parameters of a camera, and the linear approach to
+    geometric camera calibration.
+2.  Given an image *test_image.bmp* captured by a camera which is posed
+    toward a 3D chess board, we manually selected 27 points whose 2-D
+    image coordinates are saved in the file *observe.dat*, and whose 3-D
+    coordinates in a world coordinate are stored in the file
+    *model.dat*. Please use the linear approach to calibrate the camera
+    system by computing the projection matrix (M) from which you are
+    required to compute the intrinsic and extrinsic parameters,
+    including , u0, v0, α, β, and the rotation matrix (R), and the
+    shift vector (t).
+3.  In order to verify the correctness of the camera calibration, please
+    use the projection matrix to map the three set of 3D points, whose
+    3-D coordinates are {x,y,0\|x,y=0,...,10}, {x,10,z\|x,z=0,...,10},
+    {10,y,z\|y,z=0,...,10}, respectively, into the 2-D image space.
+    Discuss your results.
+4.  Create a video file that shows a 3D object moving in the 3D scene
+    along a specific pre-defined path. For example, a 3D cube (1x1x1)
+    can be displayed by nine lines connecting seven vertexes (see the
+    example below). Each line can be drawn with around 100 samples. (You
+    can see a video example from Slide 13 of the Lecture 8 handout in
+    the slideshow mode.)
+
+**Report Requirements:**
+
+1.  Briefly discuss the basics of geometric camera modeling.
+2.  Briefly discuss the linear approach to camera calibration.
+3.  Show the simulation results by figures (Part 3).
+4.  Include some sample frames of the video (Part 4) in the report.
+5.  The Python source code should be included as the appendix of the
+    report with detailed comments. A separate Python file should also be
+    provided for testing.
+6.  Zip all files (DOC, AVI, M-file) into one package and upload it to
+    Blackboard by the due date.
+
+**Useful Python functions:**
+
+You need to import the required libraries:
+
+> import torch\
+> import numpy as np\
+> import matplotlib.pyplot as plt\
+> from PIL import Image\
+> import time
+
+> import cv2\
+> import os
+
+- Inner product: torch.dot(t1,t2), Cross product : torch.cross(t1,t2).
+
+> v0 = (s \*\* 2) \* torch.dot(a2, a3)\
+> \
+> cross_a1_a3 = torch.cross(a1, a3)
+
+> 
+
+- Singular value decomposition (SVD) :
+
+> \# \-\-- Solve for the projection matrix M \-\--\
+> \# The solution to Qm=0 is the last column of V from the SVD of Q\
+> \# In PyTorch, torch.linalg.svd returns U, S, Vh (V transpose)\
+> U, S, Vh = torch.linalg.svd(Q)\
+> \
+> \# The solution is the last row of Vh, which corresponds to the
+> smallest singular value\
+> m = Vh\[-1, :\]\
+> \
+> \# Reshape the 12x1 vector m into the 3x4 projection matrix M\
+> M = m.reshape(3, 4)
+
+- Norm of a vector: torch.linalg.norm(vector).
+
+> torch.linalg.norm(cross_a1_a3)
+
+> 
+
+- Inverse of a matrix: torch.linalg.inv(M)
+
+> torch.linalg.inv(M)
+
+> 
+
+- Reshape of a vector to a matrix:
+
+> \# Reshape the 12x1 vector m into the 3x4 projection matrix M\
+> M = m.reshape(3, 4)