find-ellipse-axes.js

function getDescription() {
  return `
    <p>
      Use the sx/sy sliders to adjust the x/y scale of our unit circle,
      and use the shx/shy sliders to adjust the x/y shear components. The
      code will generate three points on the circle (aqually spaced
      angularly), transform those using the slider parameters, then 
      construct a covariance matrix, calculate the eigenvalues and
      eigenvectors, and use those to determine what the ellipse looks
      like that we get when we transform the entire circle according
      to the parameters from our sliders.
    </p>
    <p>
      The base "unit circle" is drawn in black, the transformed
      circle (minus the shear) is drawn in blue, and the fully
      transformed result, which is an ellipse, is drawn in red,
      with the ellipse axes drawn in purple for the major axis,
      and green for the minor axis.
    </p>
    `;
}

const { default: numeric } = await import(
  "https://cdn.jsdelivr.net/npm/[email protected]/+esm"
);

const scaleMatrix = new Matrix([[1, 0], [0, 1]]);
const shearMatrix = new Matrix([[1, 0], [0, 1]]);
const rotationMatrix = new Matrix([[1, 0], [0, 1]]);

function setup() {
  setSize(600, 400);
  setBorder(1, `black`);
  setGrid(20, `grey`);
  addSlider(`sx`, { min: 0, max: 200, value: 100, step: 1, transform: v => scaleMatrix.set(0, 0, v) });
  addSlider(`sy`, { min: 0, max: 200, value: 100, step: 1, transform: v => scaleMatrix.set(1, 1, v) });
  addSlider(`shx`, { min: -3, max: 3, value: 1, step: 0.01, transform: v => shearMatrix.set(0, 1, v) });
  addSlider(`shy`, { min: -3, max: 3, value: 0, step: 0.01, transform: v => shearMatrix.set(1, 0, v) });
  addSlider(`angle`, { min: -PI, max: PI, value: 0, step: 0.01, transform: (v, c = cos(v), s = sin(v)) => rotationMatrix.setData([[c, -s], [s, c]]) });
}

function draw() {
  makeEmpty();
  drawBaseShapes();
  const { a, b, v1, v2 } = fitEllipse();
  drawAxis(a, v1, `purple`);
  drawAxis(b, v2, `green`);
}

// We only need three coordinates for this to work,
// so we can pretty much hard-code everything:
function fitEllipse() {
  const θ = TAU / 3, θ2 = 2 * θ;
  const [[x1, y1], [x2, y2], [x3, y3]] = [
    transformCoords(1, 0),
    transformCoords(cos(θ), sin(θ)),
    transformCoords(cos(θ2), sin(θ2)),
  ];

  // No need for loops:
  const X = (x1 * x1 + x2 * x2 + x3 * x3) / 3;
  const Y = (y1 * y1 + y2 * y2 + y3 * y3) / 3;
  const XY = (x1 * y1 + x2 * y2 + x3 * y3) / 3;

  // Get the eigenvalues:
  const trace = X + Y;
  const gap = sqrt((X - Y) ** 2 + (2 * XY) ** 2);
  let λ1 = (trace - gap) / 2;
  let λ2 = λ1 + gap;

  // Ensure 1=major, 2=minor:
  if (λ1 < λ2) [λ1, λ2] = [λ2, λ1];

  // Then get the associated eigenvectors:
  const v1 = normalize(X - λ2, XY);
  const v2 = normalize(X - λ1, XY);

  // And then we're done, we have our semi-axes:
  const a = sqrt(2 * λ1);
  const b = sqrt(2 * λ2);
  return { a, b, v1, v2 };
}

// ===================================================

function makeEmpty() {
  clear(`white`);
  center();
  setStroke(`black`);
  noFill();
}

function drawBaseShapes(base = [], ellipse = []) {
  // base circle
  setStroke(`black`);
  circle(0, 0, 100);
  showVector([100, 0]);
  showVector([0, 100]);

  // generate transformed data
  for (let t = 0; t <= TAU; t += 0.1) {
    base.push(baseTransform(cos(t), sin(t)));
    ellipse.push(transformCoords(cos(t), sin(t)));
  }

  // transformed circle
  setStroke(`blue`);
  poly(base);
  showVector(baseTransform(1, 0));
  showVector(baseTransform(0, 1));

  // resultant ellipse
  setStroke(`red`);
  poly(ellipse);
  showVector(transformCoords(1, 0));
  showVector(transformCoords(0, 1));
}

function drawAxis(m, [x, y], color) {
  const [ax, ay] = [m * x, m * y];
  setStroke(color);
  line(ax, ay, -ax, -ay);
}

function showVector(p) {
  line(0, 0, ...p);
  circle(...p, 3);
}

function baseTransform(x, y) {
  return rotationMatrix.transform(scaleMatrix.transform([x,y]));
}

function transformCoords(x, y) {
  return baseTransform(...(shearMatrix.transform([x,y]));
}

function normalize(x, y) {
  const m = (x * x + y * y) ** 0.5;
  return [x / m, y / m];
}