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.