Как эффективно нарисовать изолинии 2D нормального распределения при определенных значениях функции плотности?
У меня есть двумерное нормальное распределение, представленное его средним значением и ковариационной матрицей. Теперь я хочу нарисовать линию вокруг каждой точки, где функция плотности
превышает определенное пороговое значение. (Термин нормализации опущен, поэтому пороговое значение может применяться ко всему распределению независимо от его размера.)
Конечно, это может быть сделано путем перебора всего изображения и проверки формулы для каждого пикселя (см. drawVariant2()). Тем не менее, это невероятно медленно, и я хочу, чтобы механизм рисования был более или менее в режиме реального времени.
Другой метод, который я нашел, вычисляет собственные значения и векторы ковариационной матрицы и использует их для преобразования изображения круга (см. drawVariant1()). Проблема в том, что у меня нет представления о том, что означает расстояние нарисованной линии от среднего в терминах функции плотности.
Можно ли как-нибудь нарисовать изолинию моего распределения при определенных значениях функции плотности?
Это код двух методов рисования, которые у меня есть до сих пор. Он использует собственную библиотеку для матричных и векторных операций.
#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;
};
Мне жаль, что код довольно длинный, но я надеюсь, что приведение полного примера поможет ответить на этот вопрос.
Решение
Что ж, после инвестирования в этот день я думаю, что наконец нашел подходящее решение.
Идея состоит в том, чтобы вычислить точку на оси X, которая имеет желаемую вероятность w.r.t. стандартному нормальному распределению. Это можно сделать с помощью алгоритма деления пополам. После этого я могу использовать расстояние этой точки в качестве радиуса круга, который трансформируется с помощью среднего и ковариационной матрицы.
Весь исходный код здесь. Я надеюсь, что это поможет, если кто-то когда-либо столкнется с подобной проблемой. Если кто-нибудь найдет какие-то ошибки, я буду рад их исправить.
#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;
};
Другие решения