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;
};