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

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

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

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

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

Поворот


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

Рассмотрим, например, поворот, простейший геометрический предикат. Этот предикат на вход принимает три точки плоскости и возвращает 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 привел пример, когда имеет смысл перестроить фильтры для увеличения производительности.
Share post

Similar posts

Comments 14

    0
    Круто. Внятно и толково изложили материал, знать и разбираться в котором должен в идеале уметь любой программист.
      +3
      Зачем это любому программисту?
        0
        Затем, чтобы понимать, как устроена арифметика с плавающей точкой. Важно не столько умение реализовать вышеописанный предикат, сколько понимание того, как работать с дробными числами, чтобы равные числа были равными.
          +2
          Программист за всю жизнь может написать много неплохих программ и ни разу не столкнуться с дробными или вещественными числами. Всё зависит от того, по какому профилю он работает.
            0
            чтобы равные числа были равными


            Это относится не только к числам с плавающей точкой.
        0
        Хорошо написано!
          +2
          Еще можно указать, что если на входе у алгоритма рациональные числа, то сравнить результат с нулем можно с помощью многомодульной арифметики — посчитать значение в конечных полях по нескольким небольшим простым модулям (5-10 чисел в окрестности 32749 подойдут), и если каждый результат оказался либо равным нулю, либо неопределенным (из-за деления на 0), то наши точки наверняка лежат на прямой, и интервальная арифметика ничего не даст. Если же результат отличен от нуля, можно запускать интервальную, а потом и длинную арифметику. Многомодульная арифметика работает гораздо быстрее, чем длинная рациональная. К сожалению, определять с ее помощью знак числа напрямую не очень просто, обычно это делают, когда известно, что результат — целое число (а промежуточные значения рациональны). Например, при вычислении определителя целочисленных матриц.
            0
            В общем случае проверка на равенство с нулем ничего не дает: вероятность того, что три произвольные точки лежат на одной прямой, ничтожна. А как посчитать знак выражения в модульной арифметике, я пока представить не могу.

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

            Может, я неверно понял ваш комментарий?
              0
              Нет, все правильно. Но вероятность того, что точки окажутся на одной прямой, или в одной плоскости, в некоторых случаях огромна. Например, если решается задача, связанная с разбиениями пространства на многогранники. Там постоянно возникают вопросы: образуют ли такие-то точки одну грань, проходит ли очередная плоскость через уже построенную вершину, являются ли два многогранника соседними (хранить эту информацию в виде какого-нибудь мультиграфа намного сложнее, чем держать просто набор координат вершин и плоскостей). Чаще всего хватает правильного выбора эпсилона, но в ответственных случаях могут потребоваться точные ответы.
                0
                Я понял, спасибо. Дам ссылку в статье на ваш комментарий.
            0
            Часть картинок умерло. Обновите статью, пожалуйста.
              0
              Не вижу. Может, дело в том, что к статье давно не обращались и latex-renderer не успел сгенерировать формулы? Сейчас, вроде бы, всё хорошо.
                0
                Хм. Похоже, что проблема в мобильном приложении. В браузере нормально. Сейчас напишу баг-репорт.
                  0
                  А всё равно, может быть сделать картинки?

              Only users with full accounts can post comments. Log in, please.