How to efficiently draw the isolines of a 2D Normal Distribution at specific values of the density function?

Question

I have a two dimensional normal distribution represented by its mean and covariance matrix. Now i want to draw a line around every point where the density function

exceeds a certain threshold value. (The normalization term is ommited so the threshold can be applied to all distribution regardless of their size.)

Of course, this can be done by iterating over the entire image and checking the formula for every pixel (see drawVariant2()). However, this is incredibly slow and I want the drawing mechanism to be more or less real-time capable.

Another method I have found computes the eigenvalues and -vectors of the covariance matrix and uses them to transform an image of a circle (see drawVariant1()). The problem with this is that I don't have a representation of what the distance of the drawn line from the mean means in terms of the density function.

Is there any way I can draw the isoline of my distribution at certain values of the density function?

This is the code of the two drawing methods I have up to now. It uses the Eigen Library for the matrix and vector operations.

#include "Eigen/Core"
#include "Eigen/Eigenvalues"

unsigned int width = 500, height = 500;
Eigen::Matrix<double, 2, 1> mean;
Eigen::Matrix<double, 2, 2> cov;

const double C_PI = 3.14159265358979323846;

void drawVariant1(unsigned char * img)
{
    const int isolineRadius = 3;
    const double baseCircleSteps = 100;
    const double circleLength = 2 * C_PI * isolineRadius * baseCircleSteps;

    Eigen::EigenSolver<Eigen::Matrix<double, 2, 2>> eigenSolver(cov);

    eigenSolver.compute(cov);

    const double eVal1 = eigenSolver.eigenvalues().real()(0);
    const double eVal2 = eigenSolver.eigenvalues().real()(1);

    const Vector2 eVec1 = eigenSolver.eigenvectors().real().col(0);
    const Vector2 eVec2 = eigenSolver.eigenvectors().real().col(1);

    for (double phi = 0; phi < 2 * C_PI; phi += 2 * C_PI / circleLength)
    {
        const double x = isolineRadius * std::cos(phi);
        const double y = isolineRadius * std::sin(phi);

        const Vector2 posProjected = x * std::sqrt(eVal1)  * eVec1  + y * std::sqrt(eVal2) * eVec2 + mean;

        const int xP = (int)posProjected(0);
        const int yP = (int)posProjected(1);

        if (xP >= 0 && xP < width && yP >= 0 && yP < height)
        {
            img[yP * width + xP] = 255;
        }
    }
}


void drawVariant2(unsigned char * img)
{
    const double threshold = 0.5;
    for(int x = 0; x < width; ++x)
    {
        for(int y = 0; y < height; ++y)
        {
            const Eigen::Matrix<double, 2, 1> point((double)x, (double)y);
            const double likelihood = std::exp( -0.5 * (point - mean).transpose() * cov.inverse() * (point - mean) );

            if(likelihood >= threshold)
                img[y * width + x] = 128;
        }
    }
}


void main()
{
    mean << 100, 100;
    cov << 50, 0, 0, 20;

    unsigned char * img = new unsigned char[width * height];
    memset(img, 0, width*height);

    drawVariant1(img);
    drawVariant2(img);

    writePGM("test.pgm", img, width, height);

    delete [] img;
};

I am sorry that the code is rather long, but I hope providing a complete example might help answering this question.

Answer 1

Well, after investing this day I think I finally found a suitable solution.

The idea is to compute the point on the x-axis which has the desired likelihood w.r.t. to the standard normal distribution. This can be done with a bisection-algorithm. Afterwards I can use the distance of this point as the radius of the circle which is transformed by the means of the mean and covariance matrix.

The entire source code is here. I hope this helps if someone might ever encounter a similar problem. If anyone finds some errors I would be happy to correct them.

#include "Eigen/Core"
#include "Eigen/Eigenvalues"

#include "PgmIO.h"

typedef Eigen::Matrix<double, 2, 1> Vector;
typedef Eigen::Matrix<double, 2, 2> Matrix;

const double C_PI = 3.14159265358979323846;

Vector mean;
Matrix cov;

int width = 500, height = 500;

void draw(unsigned char * img, double threshold)
{
    const double epsilon = 0.001;

    Matrix unity;
    Vector base;

    base << 1, 0;
    unity << 1, 0, 0, 1;

    double drawingDistance = 0;

    // Search for the Point on the X-Axis that has the desired likelihood with a bisection-method
    for(double a = 0, b = 5, likelihood = 0; std::abs(likelihood - threshold) > epsilon; )
    {
        const double currentDistance = (a + b) / 2;
        const Vector evalPoint = currentDistance * base;

        likelihood = std::exp(-0.5 * evalPoint.transpose() * unity * evalPoint); //unitiy.inverse() == unity

        // Suitable point found
        if(std::abs(likelihood - threshold) < epsilon)
        {
            drawingDistance = currentDistance;
            break;
        }

        if(likelihood > threshold)
        {
            a = currentDistance; // If the likelihood is too large search further away from the origin
        }
        else
        {
            b = currentDistance; // If the likelihood is too small search closer to the origin
        }
    }

    Eigen::EigenSolver<Matrix> eigenSolver(unity);

    eigenSolver.compute(cov);

    const double eVal1 = eigenSolver.eigenvalues().real()(0);
    const double eVal2 = eigenSolver.eigenvalues().real()(1);

    const Vector eVec1 = eigenSolver.eigenvectors().real().col(0);
    const Vector eVec2 = eigenSolver.eigenvectors().real().col(1);

    const double baseCircleSteps = 100;
    const double circleLength = 2 * C_PI * drawingDistance * baseCircleSteps;

    for (double phi = 0; phi < 2 * C_PI; phi += 2 * C_PI / circleLength)
    {
        // Compute the points on a circle within drawingDistance range
        const double x = drawingDistance * std::cos(phi);
        const double y = drawingDistance * std::sin(phi);

        // Project point to the eqivalent point on the isoline
        const Vector posProjected = x * std::sqrt(eVal1)  * eVec1  + y * std::sqrt(eVal2) * eVec2 + mean;

        const int xP = (int)posProjected(0);
        const int yP = (int)posProjected(1);

        // Set point in the image
        if (xP >= 0 && xP < width && yP >= 0 && yP < height)
        {
            img[yP * width + xP] = 255;
        }
    }
}


void main()
{
    mean << 100, 100;
    cov << 800, 100, 100, 500;

    unsigned char * img = new unsigned char[width * height];
    memset(img, 0, width*height);

    draw(img, 0.01);

    writePGM("img.pgm", img, width, height);

    delete [] img;
};