Я отлаживаю написанный мной алгоритм, который преобразует сложный самопересекающийся многоугольник в простой многоугольник путем идентификации всех точек пересечения и обхода внешнего периметра.
Я написал серию случайных данных, генерирующих стресс-тесты, и в этом я столкнулся с интересной ситуацией, которая приводила к сбою моего алгоритма после правильной работы тысячи и тысячи раз.
double fRand(double fMin, double fMax)
{
double f = (double)rand() / RAND_MAX; // On my machine RAND_MAX = 2147483647
return fMin + f * (fMax - fMin);
}
// ... testing code below:
srand(1);
for(int j=3;j<5000;j+=5) {
std::vector<E_Point> geometry;
for(int i=0;i<j;i++) {
double radius = fRand(0.6,1.0);
double angle = fRand(0,2*3.1415926535);
E_Point pt(cos(angle),sin(angle));
pt *= radius;
geometry.push_back(pt); // sending in a pile of shit
}
// run algorithm on this geometry
Это был обзор того, насколько это актуально и как я попал туда, где я сейчас.Есть намного больше деталей, которые я опускаю.
Что я смог сделать, так это сузить проблему до кода пересечения сегмента с сегментом, который я использую:
bool intersect(const E_Point& a0, const E_Point& a1,
const E_Point& b0, const E_Point& b1,
E_Point& intersectionPoint) {
if (a0 == b0 || a0 == b1 || a1 == b0 || a1 == b1) return false;
double x1 = a0.x; double y1 = a0.y;
double x2 = a1.x; double y2 = a1.y;
double x3 = b0.x; double y3 = b0.y;
double x4 = b1.x; double y4 = b1.y;
//AABB early exit
if (b2Max(x1,x2) < b2Min(x3,x4) || b2Max(x3,x4) < b2Min(x1,x2) ) return false;
if (b2Max(y1,y2) < b2Min(y3,y4) || b2Max(y3,y4) < b2Min(y1,y2) ) return false;
float ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3));
float ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3));
float denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1);
// check against epsilon (lowest normalized double value)
if (fabs(denom) < DBL_EPSILON) {
//Lines are too close to parallel to call
return false;
}
ua /= denom;
ub /= denom;
if ((0 < ua) && (ua < 1) && (0 < ub) && (ub < 1)) {
intersectionPoint.x = (x1 + ua * (x2 - x1));
intersectionPoint.y = (y1 + ua * (y2 - y1));
return true;
}
return false;
}
Что происходит, так как у меня есть два пересечения, которые возвращает эта функцияточно такое же значение точки пересечения для.Вот очень увеличенное изображение соответствующей геометрии:
Вертикальная линия определяется точками (0.3871953044519425, -0.91857980824611341), (0.36139704793723609, 0.91605957361605106)
зеленая почти горизонтальная линияпо точкам (0.8208980020500205, 0.52853407296583088), (0.36178501611208552, 0.88880385168617226)
и белой линии по (0.36178501611208552, 0.88880385168617226), (-0.43211245441046209, 0.68034202227710472)
Как видно, две последние линии имеют общую точку.
Моя функция дает мне решение (0.36178033094571277, 0.88880245640159794) для ОБА этих пересечений (одно из которых вы видите на рисунке красной точкой).
Причина, по которой это большая проблема, заключается в том, что мой алгоритм периметра зависит от сортировки точек пересечения на каждом ребре.Поскольку обе эти точки пересечения были рассчитаны на одинаковое точное значение, сортировка поместила их в неправильную ориентацию.Путь по периметру идет сверху и следовал за белой линией слева, а не зеленой линией слева, что означает, что я больше не следую за внешним периметром моей формы многоугольника.
Чтобы исправить эту проблему, я, вероятно, могу сделать много вещей, но я не хочу искать во всех списках пересечений, чтобы проверить по другим точкам, чтобы увидеть, совпадают ли позиции.Лучшим решением было бы попытаться повысить точность моей функции пересечения.
Итак, что я спрашиваю, почему точка решения такая неточная?Это потому, что одна из линий почти вертикальная?Должен ли я сначала выполнить какое-то преобразование?В обоих случаях почти вертикальная линия была передана как a0 и a1.
Обновление: Эй, посмотрите на это:
TEST(intersection_precision_test) {
E_Point problem[] = {
{0.3871953044519425, -0.91857980824611341}, // 1559
{0.36139704793723609, 0.91605957361605106}, // 1560
{-0.8208980020500205, 0.52853407296583088}, // 1798
{0.36178501611208552, 0.88880385168617226}, // 1799
{-0.43211245441046209, 0.6803420222771047} // 1800
};
std::cout.precision(16);
E_Point result;
intersect(problem[0],problem[1],problem[2],problem[3],result);
std::cout << "1: " << result << std::endl;
intersect(problem[0],problem[1],problem[3],problem[2],result);
std::cout << "2: " << result << std::endl;
intersect(problem[1],problem[0],problem[2],problem[3],result);
std::cout << "3: " << result << std::endl;
intersect(problem[1],problem[0],problem[3],problem[2],result);
std::cout << "4: " << result << std::endl;
intersect(problem[2],problem[3],problem[0],problem[1],result);
std::cout << "rev: " << result << std::endl;
intersect(problem[3],problem[2],problem[0],problem[1],result);
std::cout << "revf1: " << result << std::endl;
intersect(problem[2],problem[3],problem[1],problem[0],result);
std::cout << "revf2: " << result << std::endl;
intersect(problem[3],problem[2],problem[1],problem[0],result);
std::cout << "revfboth: " << result << std::endl;
}
Вывод:
Starting Test intersection_precision_test, at Polygon.cpp:1830
1: <0.3617803309457128,0.8888024564015979>
2: <0.3617803309457128,0.8888024564015979>
3: <0.3617803314022162,0.8888024239374175>
4: <0.3617803314022162,0.8888024239374175>
rev: <0.3617803635476076,0.8888024344185281>
revf1: <0.3617803313928456,0.8888024246235207>
revf2: <0.3617803635476076,0.8888024344185281>
revfboth: <0.3617803313928456,0.8888024246235207>
На самом деле у меня заканчиваются биты мантиссы или я могу значительно лучше работать с более умным алгоритмом?
Проблема в том, что у меня нет простого способа определить, когда вершина установлена действительно очень близко кдругая строка.Я не возражал бы переместить это или даже полностью обстреливать это, потому что ни один из них не испортит мой вид!