2D Gaussian Tight Bounding Quad

by Matt

HTML

<base href="https://rawcdn.githack.com/mrdoob/three.js/r156/examples/" />
<script async src="https://cdn.jsdelivr.net/npm/[email protected]/dist/es-module-shims.js"></script>
<script type="importmap">
	{
		"imports": {
			"three": "https://cdn.jsdelivr.net/npm/[email protected]/build/three.module.js",
			"three/addons/": "https://cdn.jsdelivr.net/npm/[email protected]/examples/jsm/",
			"lil-gui": "https://cdn.jsdelivr.net/npm/[email protected]/dist/lil-gui.esm.min.js"
		}
	}
</script>

CSS

canvas {
	position: fixed;
	inset: 0;
}

JavaScript

import * as THREE from 'three'
import { OrbitControls } from 'three/addons/controls/OrbitControls.js'
import { GUI } from 'lil-gui'

/**
 * Supporting example for https://twitter.com/the_ross_man/status/1704331583786778895
 * by Matt Rossman
 *
 * This renders a Gaussian from the covariance matrix in a (roughly) fullscreen shader,
 * and draws the bounds as a transformed quad to show that they line up.
 * 
 * In practice you could draw the Gaussian on the rotated quad and simplify the shader
 * to a unit Gaussian.
 */

let scene, camera, renderer
let params, uniforms
let meshGaussian, meshBounds

init().then(animate)

/**
 * Given a covariance matrix:
 * 
 * | a  b |
 * | c  d |
 *
 * Decompose this into eigenvectors and eigenvalues.
 * Eigenvectors point along the axes of the ellipse. Their magnitude doesn't matter.
 * Eigenvalues represent the square of the scaling factors along those axes.
 * 
 * From these we can extract rotation and scaling values to transform our quad.
 *
 * References:
 * - https://www.youtube.com/watch?v=e50Bj7jn9IQ
 * - https://en.wikipedia.org/wiki/Eigenvalue_algorithm#2%C3%972_matrices
 * - https://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html
 */
function decomposeCovariance(a, b, d) {
  
  const det = a * d - b * b;  // matrix is symmetric, so "c" is same as "b"
	const trace = a + d;

	const mean = 0.5 * trace;
	const dist = Math.sqrt(mean * mean - det);

	const lambda1 = mean + dist;  // 1st eigenvalue
	const lambda2 = mean - dist;  // 2nd eigenvalue
  
  let v1; // 1st eigenvector
  let v2; // 2nd eigenvector
  
  if (b === 0) {
  	// https://twitter.com/the_ross_man/status/1706342719776551360
  	if (a > d) v1 = [1, 0];
    else v1 = [0, 1];
    
  } else v1 = [b, d - lambda2];
  
  // We don't actually need this, but for posterity the 2nd eigenvector is just a
  // 90 degree rotation of the first since Gaussian axes are orthogonal
  
  v2 = [v1[1], -v1[0]];
  
  // angle
  
  const theta =...