* Simple Line Fitting

 * Keith J. O'Hara <kohara2@swarthmore.edu>

 * March 2024

 */

​

// uses ML.js https://github.com/mljs/ml?tab=readme-ov-file

// because jsfeat svd is broken :(

​

// the points we are fitting

let p = [];

​

const USE_SVD = false;

​

function setup() {

  createCanvas(640, 480);

  noStroke();

  background(255);

  print(ML.SVD);

}

​

function mousePressed() {

  p.push(new p5.Vector(mouseX, mouseY, 1));

}

​

function estimate_least_squares() {

  let x = new p5.Vector(0, 0, 0);

​

  if (p.length >= 2) {

    as = [];

    bs = [];

    // Creates the estimation matrix

    for (let i = 0; i < p.length; i++) {

      as.push([p[i].x, 1]);

      bs.push([p[i].y]);

    }

    let mls_a = new ML.Matrix(as);

    let mls_b = new ML.Matrix(bs);

    let mls_x = ML.MatrixLib.solve(mls_a, mls_b);

    let a = mls_x.data[0][0];

    let b = -1;

    let c = mls_x.data[1][0];

    x = new p5.Vector(a, b, c);

    //print(x);

  }

  return x;

}

​

function estimate_total_least_squares() {

  let x = new p5.Vector(0, 0, 0);

​

  if (p.length >= 2) {

    as = [];

    bs = [];

    // Creates the estimation matrix

    for (let i = 0; i < p.length; i++) {

      as.push([p[i].x, p[i].y, 1]);

    }

    if (USE_SVD) {

      let svd = new ML.SVD(as);

      let a = svd.V.data[0][2];

      let b = svd.V.data[1][2];

      let c = svd.V.data[2][2];

      x = new p5.Vector(a, b, c);

      //print("SVD " + svd.V + "\n" + x);

    } else {

      let A = new ML.Matrix(as);

      let E = new ML.EVD(A.transpose().mmul(A));

      a = E.V.data[0][0];

      b = E.V.data[1][0];

      c = E.V.data[2][0];

      x = new p5.Vector(a, b, c);

      //print("E:" + E.V + "\n" + x);

    }

  }

  return x;

}

​

function draw() {

  background(255);

  stroke(0);

  fill(0);

  for (let i = 0; i < p.length; i++) {

    ellipse(p[i].x, p[i].y, 5, 5);

  }

  if (p.length >= 2) {

    stroke(255, 0, 0);

    plotLine(estimate_least_squares());

    stroke(0, 0, 255);

    plotLine(estimate_total_least_squares());

  }

}

​

function plotLine(v) {

  if (abs(v.x) < abs(v.y)) {

    // mostly horizontal

    line(0, v.z / -v.y, width, (v.x * width + v.z) / -v.y);

  } else {

    //mostly vertical

    line(v.z / -v.x, 0, (v.y * height + v.z) / -v.x, height);

  }

}

​