Pull to refresh

Точное вычисление геометрических предикатов

Reading time 7 min
Views 9.8K
Доброго вам дня, коллеги. Предлагаю вам прочитать статью о базовом аспекте вычислительной геометрии — точном вычислении предикатов.

Возможно, многие из вас сталкивались со «странным» поведением реализаций алгоритмов вычгеома: падениями и некорректными результатами. В общем, ничего удивительного в этих странностях нет — это издержки переноса непрерывной геометрии в суровые дискретные реалии арифметики с плавающей точкой. Рискну предположить, что кто-то из вас, отлаживая свои алгоритмы, решал подобные проблемы, подбирая магические константы. Скорее всего, если этот способ и приводил к результату, то, вероятно, к временному.

В этой статье я расскажу, как избавиться от недостатков вычислений с плавающей точкой практически без ущерба для эффективности. План статьи:
  • геометрический предикат «поворот»;
  • противоречивость «наивной» реализации поворота;
  • применение интервальной арифметики для упрощения расчета предикатов;
  • длинная арифметика и несколько альтернативных способов вычислений предикатов;
  • расчет погрешности вычисления поворота в числах с плавающей точкой.

Итак, если вы хотите знать, почему вычислительная геометрия не является «наукой о подборе эпсилонов» и как можно корректно и эффективно реализовать геометрический алгоритм, жмите сюда.

Поворот


Вообще, любой алгоритм, как дискретная структура, реализуется строго. Проблемы возникают в условных блоках, где проверяется значение какого-нибудь предиката (например, лежит ли точка внутри треугольника, пересекает ли отрезок четырехугольник и т.д.).

Рассмотрим, например, поворот, простейший геометрический предикат. Этот предикат на вход принимает три точки плоскости и возвращает 1, если c лежит слева от направленного отрезка ab; -1, если справа; и 0, если три точки лежат на одной прямой. Предикат реализуется с помощью векторного произведения:



Этот предикат, несмотря на свою простоту, очень важен: с его помощью можно проверить, например, пересекается ли пара отрезков, лежит ли точка внутри треугольника и т.п. Рассмотрим его реализацию:

enum turn_t {left = 1, right = -1, collinear = 0};

double cross(point_2 const & a, point_2 const & b)
{
   return a.x * b.y - b.x * a.y;
} 

turn_t turn(point_2 const & a, point_2 const & b, point_2 const & c)
{
   double det = cross(b - a, c - a);
   if (det > 0)
      return left;
   if (det < 0)
      return right;
      
   return collinear; 
}

Рассмотрим вызов функции turn со следующим входом:

point_2 a(3.0, 5.0);
point_2 b = -a;
point_2 c =  a * (1LL << 52);

turn_t res = turn(a, b, c); // ожидается collinear

Несмотря на то, что все точки лежат на одной прямой, результат, в силу ограниченной точности арифметики с плавающей точкой, не равен collinear. Часто, в качестве hot fix, полагают, что векторное произведение может считаться нулем, если меньше некоторой предопределенной константы e:

double det = cross(b - a, c - a);
if (det > e)
   return left;
if (det < -e)
   return right;
   
return collinear; 

Мысль неудачна, так как предикат будет несовместен уже на четырех точках: для любой такой константы e можно подобрать четыре точки так, что предикат будет возвращать противоречивые результаты. Например, рассмотрим два параллельных, не лежащих на одной прямой отрезка ab и cd (см. рисунок).



На этом рисунке треугольники и имеют одинаковую площадь, равную e/2. Поскольку удвоенная площадь треугольника равна произведению высоты на основание, все точки длинного отрезка ab образуют на основании cd треугольники одинаковой площади. Площадь эта будет меньше площадей аналогичным образом построенных треугольников на основании ab и вершинами на cd. Из этого следует, например, что при определенной константе e отрезок ab будет коллинеарен отрезку cd, но cd не будет коллинеарен отрезку ab.

Но заметим, что обычно collinear — достаточно редкий результат. Если бы мы могли большую часть левых и правых поворотов выдать, сравнив векторное произведение с некоторой константой, а оставшуюся часть ответов посчитать каким-нибудь менее эффективным, но точным способом (например, воспользовавшись длинной арифметикой), то мы могли бы считать нашу задачу выполненной, так как в среднем время расчета поворота увеличилось бы незначительно. Расчет константы e я приведу в приложении ниже, для тех, кто не испугается нудных, но тривиальных выкладок. После внесенных изменений предикат будет выглядеть так:

double l = (b.x - a.x) * (c.y - a.y);
double r = (c.x - a.x) * (b.y - a.y);

double e = (abs(l) + abs( r)) * std::numeric_limits<double>::epsilon() * 4;

double det = l - r;

if (det > e)
   return left;
if (det < -e)
   return right;
   
long_point_2 la(a), lb(b), lc( c);
long_t ldet = cross(lb - la, lc - la);
if (ldet > 0)
   return left;
if (ldet < 0)
   return right;
   
return collinear;

В общем, мы уже получили абсолютно корректный и достаточно быстрый предикат. Есть две тонкости: во-первых, его можно еще немного ускорить, а, во-вторых, вывод формулы для константы e в общем случае — занятие достаточно тяжелое, несмотря даже на существенную помощь со стороны символьной арифметики (sympy, sage и т.д.). Тут самое время вспомнить об интервальной арифметике.

Интервальная арифметика


Основная идея интервальной арифметики в том, чтобы зафиксировать вещественное число в некотором отрезке с floating-point границами. Любое точно представимое в плавающей арифметике число будет представлено вырожденным отрезком — точкой.

При сложении двух интервалов складываются их границы: верхние границы с округлением вверх (к плюс бесконечности), а нижние — с округлением вниз (к минус бесконечности). Несложно определить остальные арифметические операции из следующих соображений: точное значение вещественного числа неизвестно, оно может быть любым в отрезке. Значит, результирующий отрезок должен содержать результаты арифметического действия над всеми парами чисел из отрезков-операндов.

Существует несколько реализаций интервальной арифметики, например boost/interval.

Используя интервальную арифметику, можно намного точнее определить «доверительный интервал» расчета предиката в арифметике с плавающей точкой. Интервальная арифметика по понятным причинам в два-пять раз медленнее обычной арифметики, но порядка на три-четыре быстрее длинной арифметики. Так что в реализации предикатов имеет смысл уточнять доверительный интервал после обычной арифметики с плавающей точкой перед вызовом длинной арифметики:

double l = (b.x - a.x) * (c.y - a.y);
double r = (c.x - a.x) * (b.y - a.y);

double e = (abs(l) + abs( r)) * std::numeric_limits<double>::epsilon() * 4;

double det = l - r;

if (det > e)
   return left;
if (det < -e)
   return right;

interval_point_2 ia(a), ib(b), ic( c);
interval<double> idet = cross(ib - ia, ic - ia);

// если отрезок не содержит нуля, значит, 
// знак определен достоверно
if (!zero_in(idet))
{
   if (idet > 0)
      return left;
      
   return right;
}
   
long_point_2 la(a), lb(b), lc( c);
long_t ldet = cross(lb - la, lc - la);
if (ldet > 0)
   return left;
if (ldet < 0)
   return right;
   
return collinear;

Обратите внимание на то, что если убрать обычную арифметику перед интервальной, скорость расчетов упадет в два-пять раз, но не на три-четыре порядка, как в случае, если бы мы оставили только длинную арифметику. Иногда оценить доверительный интервал предиката существенно труднее, чем в данном случае, а производительность длинной арифметики не позволяет использовать ее в качестве основного средства вычислений. В этом случае интервальная арифметика будет являться разумным компромиссом.

О длинной арифметике


Вообще, длинная арифметика — не единственный способ точно определить знак арифметического выражения. Есть несколько алгоритмов эффективнее длинной арифметики на порядок-полтора. Речь идет об алгоритмах ESSA и Adaptive Precision Arithmetic. Приводить эти алгоритмы здесь я не буду, благо в интернете легко найти подробные описания. Сделаю лишь замечание, которое может спасти некоторое время при отладке: часто флаги сопроцессора выставляются таким образом, чтобы расчеты велись в десятибайтовых вещественных числах, которые при присваивании выталкиваются с округлением в восьми- или четырехбайтовые вещественные числа. За счет этого достигается большая точность вычислений, но это негативно сказывается на упомянутых алгоритмах ESSA и Adaptive Precision Arithmetic. В остальном же эти алгоритмы вполне переносимы и достаточно просто реализуемы.

Выводы


В данной статье я вас познакомил с методом фильтрованного вычисления предиката. На первом шаге фильтра (арифметикой с плавающей точкой) эффективно отсеивается большая часть входных данных. На втором шаге (интервальная арифметика) отсеивается значительная часть входных данных, прошедших первый фильтр. На третьем шаге (длинная арифметика, ESSA или Adaptive Precision) обрабатываются оставшиеся данные, прошедшие через предыдущие этапы. На наших тестах (равномерное распределение в квадрате) их ста миллионов входных данных, примерно двести тысяч прошло на интервальную арифметику. До длинной арифметики дошло всего несколько входов, что позволяет сделать оптимистичный вывод об эффективности и простоте данного подхода. Данный подход является общепринятым: его использует, например, библиотека вычислительной геометрии CGAL. В своих задачах вы вполне можете использовать и свои фильтры, сообразно с характером своих входных данных.

Ссылки



Примечание. Вычисление константы e для поворота


Здесь уже много написано про числа с плавающей точкой, поэтому я просто перечислю некоторые факты.

Двоичное число с плавающей точкой представляется в виде: . Операции с погрешностью над числами с плавающей точкой (в противовес обычным операциям над вещественными числами) обычно обозначают так: . «Машинный эпсилон» разные авторы определяют по-разному, я воспользуюсь таким определением: . Обратите внимание, в stl, например, эпсилон в два раза больше. Тогда погрешности операций при условии округления к ближайшему можно выразить так:





Посчитаем поворот в числах с плавающей точкой:

.

Перейдем от арифметики с плавающей точкой к вещественной арифметике:



Теперь сформулируем достаточное условие корректности вычисления знака . Представив себе шар с центром в точке и радиусом e, можно убедиться в корректности этого условия: шар не будет содержать нулевой точки, следовательно, все точки этого шара будут иметь один знак, в том числе и точка v. Осталось оценить модуль разности . Раскроем скобки с дельтами и вспомним, что модуль разности не превышает суммы модулей, а модуль произведения равен произведению модулей:



Мы получили нижнюю границу числа e в вещественных числах, но e является числом с плавающей точкой. Заметим, что



Отсюда немедленно следует:



Несложно получить ограничение дроби сверху:



Значит, нужная нам константа e может быть вычислена так:



Обратите внимание на то, что последнее умножение хоть и является вещественным, не выводит число из множества чисел с плавающей точкой, так как множитель — степень двойки.

UPD. Спасибо Portah и fsgs за указанную опечатку: функция cross возвращает не point_2, а double.
UPD2. Mrrl привел пример, когда имеет смысл перестроить фильтры для увеличения производительности.
Tags:
Hubs:
+49
Comments 14
Comments Comments 14

Articles