I've been working towards a Gaussian splat renderer. But first I wanted a basic 2D version to learn on. It's actually remarkably simple. A Gaussian, by definition, is $e^{-x^2}$. Moving into 2D just gives you one term per coordinate. You control it by shifting and scaling with constants, e.g. $2e^{-(x + 5)^2}$, or the classic normal distribution. In code this translates to something very simple:
struct Gaussian2D {
glm::vec2 mu;
glm::mat2 sigma;
glm::vec4 color;
};
The fancy statistical names map onto plainer ones: $\boldsymbol{\mu}$ (mu) is the position,
$\Sigma$ (sigma) is the rotation and scale. The full formula is
$$ G(\mathbf{p}) = \exp\!\left( -\tfrac{1}{2}\,(\mathbf{p} - \boldsymbol{\mu})^{\top}\, \Sigma^{-1}\, (\mathbf{p} - \boldsymbol{\mu}) \right) $$
with an optional $\frac{1}{2\pi\sqrt{|\Sigma|}}$ out front if you want it to integrate to 1. For rendering
we don't care about that, so it's dropped. In code that's:
glm::vec2 p = glm::vec2( x, y );
glm::vec2 diff = p - g.mu;
float res = glm::exp( -0.5f * glm::dot( diff, sigmaInv * diff ) );
AddPixel( image, x, y, g.color * res );
The -0.5f comes straight from the normal distribution's $e^{-x^2 / (2\sigma^2)}$. Pulling
the $\tfrac{1}{2}$ out is what keeps sigma behaving like a variance instead of half of one.
That's it, surprisingly simple. The hard part is making it fast and training the Gaussians in the first place. A topic for another day.