Russian Qt Forum

Широко известно, что сравнивать напрямую числа с плавающей точкой на равенство нельзя:

C++ (Qt)
double a;
double b;
...
if(a == b) // Wrong!!!
...
 

А вот как их сравнивать надо - освещено достаточно плохо. Наиболее распространенный (и дешевый) подход - сравнить разность двух чисел с неким малым значением:

C++ (Qt)
double a;
double b;
...
if(fabs(a - b) < epsilon) // True!!!
...
 

Вот только как выбрать эту самую epsilon? Преимущественное большинство берет значение из головы, в крайнем случае злосчастную DBL_EPSILON. Мне же хотелось бы строго математического обоснования. Кое-что близкое нашел тут (http://www.rsdn.ru/forum/cpp/2640596.1.aspx) и тут (http://stackoverflow.com/questions/3109941/when-to-use-dbl-epsilon-epsilon), но приведенное обоснование так и не понял (с двоичной арифметикой туговато, надо поднимать). Хотелось бы с вашей помощью разобраться правы или нет авторы материалов по приведенным ссылкам и если правы - то почему. А может вы знаете какой-нибудь еще (более правильный) способ сравнения вещественных чисел?

У qt-ов есть метод, но на чем они основываются хз

C++ (Qt)
Q_DECL_CONSTEXPR static inline bool qFuzzyCompare(double p1, double p2)
{
    return (qAbs(p1 - p2) <= 0.000000000001 * qMin(qAbs(p1), qAbs(p2)));
}

Вообще непонятное магическое число.

На эти грабли я наступал много раз :) Напр нужно вычислить нормаль к треугольнику. На первый взгляд дело нехитрое (треугольник задан 3-мя точками p0, p1, p2)

C++ (Qt)
Point3D nrml = CrossProduct(p1 - p0, p2 - p1);
double len = nrml.length();
if (len < 1.0e-6) ..  
 

Это нормально во многих случаях, напр для рендера: треугольник достаточно мал, значит он НЕ планарен, выкинуть его нафиг - и все дела. Но не всегда. Напр для preview выкидывать его нельзя. Приходится изощряться - попытаться увеличить исходные вектора. Если все равно не помогает - принять нормаль "от фонаря" (напр 0, 0, -1) лучшего нет. Др. словами нет какого-то общего решения, надо смотреть конкретно по задаче.

Конечно есть ф-ции nextafter (чтобы найти следующее значение) но на практике от них мало толку, epsilon почти всегда есть. Существует еще более мерзкая проблема - бой/неустойчивость флота. Ну то если интересно обсудим

Рассказывай!

Ну вроде как претендует на универсальность следующее:

C++ (Qt)
if (fabs(a - b) <= 16 * DBL_EPSILON * fmax(fabs(a), fabs(b)))
{
    // ...
}
 

Ни на что оно не претендует :) Это довольно прямолинейная попытка вычислить релятивную погрешность. Все зависит от того что Вы собираетесь делать с полученным результатом в задаче. Напр если Вам (хотя бы) потребуется возвести разницу в квадрат - получите проблемы. Более практично установить общий предел точности (напр 1.0e-6) для всех вычислений в задаче - и от него уже плясать.

Проблема в том, что мне не известен удовлетворительный предел точности. Собственно задача стоит так:
Я пробегаюсь по всем треугольникам и нахожу сегменты линии. Неприятности начинаются тогда, когда я их соединяю. А соединяю я их от одного конечного сегмента до другого, присоединяя по одному за раз. Если epsilon слишком большая - то иногда происходит "перескок" через сегмент, т.е. к текущему сегменту подсоединяется не следующий сегмент, а сегмент после следующего. В итоге линия разрывна. Если epsilon слишком маленькая, то иногда два смежных сегмента не считаются смежными. В итоге линия не достраивается.

Ну те "узлы треугольников" называются "вертексы" :) И нормально надо строить edges (ребра). Т.е. каждый треугольник знает от 0 до 3 смежных треугольников. Тогда ясно какой треугольник брать следующим. Это только на первый взгляд сложновато, но все равно на то и выйдет как только задача хоть немного усложнится. Напр Вы полагаете что у Вас всего одна поверхность, а если их 2 или больше? А если ф-ция пересекает ребро в 2 и более точках? Или просто ф-ция достаточно затратна и лупить по всем треугольникам медленновато. Типичная структура edge

C++ (Qt)
struct CEdge {
 int mVertex[2]  // индексы в контейнере вертексов 
 int mFacet[2];   // индексы треугольников которые шерят это ребро (второй -1 если это край)
};
 
struct CTriangle {
 ...
 int mEdge[3];  // индексы в контейнере CEdge
};
 

+1

Igors, поделитесь, пожалуйста, опытом.

аттапырил ушко.
а че там с qreal? такие же проблемы?
у олдтролей есть там еще один тип внутренний не помню как зовут... толи фиксед то ли еще чето. не помню.

C++ (Qt)
typedef double qreal;

Ну вот допустим мы удаляем непланарные треугольники (пост #3). И на каждом шаге (или запуске) нам приходят одни и те же исходные треугольники, но повернутые по-разному. Матрица поворота известна, мы можем их развернуть в начальное положение, но при этом данные все же отличаются от исходных (пусть в каком-то дальнем знаке). Этого достаточно чтобы уже другие треугольники удовлетворяли условию планарности. Если есть данные привязанные к ним по индексам - тогда "ой" :'( И можно сколько угодно выбирать epsilon, но всегда найдется такой который "перескочит" через него (и, как правило находится он удивительно скоро :))

Поэтому всегда лучше опираться на "порядок построения" чем "склеивать куски" - это ненадежно

За идею спасибо! В принципе, уже морально подготовился к тому, что придется переделывать все алгоритмы/структуры данных. Вот только процесс это длительный, т.к. озвученная задача - лишь один из модулей программного комплекса. А пока приходится работать с тем что есть. Так что если у кого-то есть мысли/советы/предложения по поводу первоначально поставленного вопроса - милости прошу.

В терминологии моей предметной области они называются "узлы расчетной области" :) Я численник.

epsilon обычно выбирается на основании минимальной длины ребра. Однако читая Ваш пост #7 я не пойму откуда взялась эта проблема? Ну считайте 1 ребро 1 раз - и ничего сравнивать не нужно. Да еще и расчетов в 3 раза меньше :)  

Edit: не хотите пока заводить новые структуры данных - ладно. Но никто не мешает напр сбрасывать обсчитанные ребра в локальный хэш - который живет только во время пробега по треугольникам

См. вложение. Черным нарисованы треугольники. Синим, красным и зеленым - независимо найденные сегменты линии уровня. Проблема с красным сегментом. Если epsilon слишком большая, то с концом синего сегмента "совпадают" и начало красного сегмента и начало зеленого сегмента. Сетка увеличена во много раз. В действительности координаты начал красного и зеленого сегмента могут отличаться и в восьмой значащей цифре, и в десятой... С другой стороны, координаты конца синего сегмента и координаты начала красного сегмента рассчитываются по разному: координаты синего сегмента - через координаты вершин верхнего треугольника, координаты красного - через координаты вершин среднего треугольника. Т.е., хотя теоретически эти точки совпадают, на практике их координаты не совпадают абсолютно точно. Если взять слишком маленькую epsilon, получается, что линия прерывается на синем сегменте т.к. нет ни одного сегмента начало которого совпадало бы с концом синего сегмента.

Я хочу. Но на это уйдет не меньше месяца. А в текущую программу кардинальные изменения вносить нельзя: вызовет цепную реакцию. Сейчас треугольники хранятся набором вершин и ничего о соседях не знают. Переделывать структуру данных для треугольника - переписывать, как минимум, триангулятор и конечноэлементную сетку. Что может привести к еще более обширным изменениям. Программа уже давно держится только на костылях :(.

Все понятно, но из 2 смежных треугольников один всегда будет считаться  после другого. Значит можно взять уже имеющуюся точку на общем ребре, и это совсем не сложно

Пусть "ни кола ни двора", даже индексов нет, но координаты-то всегда есть. Простой "местный" класс, напр

C++ (Qt)
struct CTempEdge {
 CTempEdge( const Point * p0, const Point * p1 ) : 
 mCalcFlag(false),
 mInters(0)
 {
  if (*p0 < *p1) {   // "меньший" вертекс первый
   mPt[0] = p0;
   mPt[1] = p1;
  }
  else {
   mPt[1] = p0;
   mPt[0] = p1;
  }
 }
 
 bool operator < ( const CTempEdge & sec ) const
 {
  if (*mPt[0] < *sec.mPt[0]) return true;
  if (*mPt[0] > *sec.mPt[0]) return false;
  return *mPt[1] < *sec.mPt[1];
 }
 
 bool mCalcFlag;    // флаг "рассчитано" 
 const Point * mInters;  // указатель на точку пересечения (или NULL)
 const Point * mPt[2];  // указатели на 2 вертекса
};
 

Оператор < для вертекса - просто последовательное сравнение для x, y, (z). Ну и любой ассоциативный контейнер (хоть std::set), и там уже блистать техникой ООП :). Не вижу как это вызовет глобальные переделки.

Также по Вашему чертежу хорошо видно что отсекаемый сегмент может быть сколь угодно мал, и, значит, меньше погрешности вычисления. Др. словами "правильная следующая точка" - необязательно "ближайшая". Любой выбор epsilon это не решает  

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

Та пусть считаются в каком угодно порядке - Вы всегда можете посмотреть было ли ребро уже обсчитано из std::set, ключ-то есть

На мой взгляд правильнее говорить, что линия должна быть неразрывной, а реально с учетом применяемого алгоритма, разрывность линии не должна превышать некоторого эпсилон. Если предметная область, в которой решается задача, не задает этого эпсилона, алгоритм нужно менять.
Другая неприятная сторона этого алгоритма, как я понял, заключается в том, что каждое ребро просчитывается 2 раза. Причем вычисленные два раза координаты точки ребра не совпадают при одинаковых параметрах, входящих в расчет и разница в результатах значительна для данной задачи. Понятно, что это не мистика, но признак того, что результат сильно зависит от ошибок "машинного вычисления" (не знаю, как точно называется это понятие).
Поэтому, если изменить алгоритм на данном этапе невозможно, предлагаю изучить формулу, по которой проводится расчет координат точки ребра. Скорее всего, формула основана на линейной интерполяции, в которой производится деление на величину (z2-z1) (это значения функции в вершинах ребра). Хорошо бы избежать операции деления, т.е. в результаты расчета заносить только числитель, и хранить дополнительно третьим параметром Z этот знаменатель. Тогда алгоритм сравнения двух точек (X1, Y1) и (X2, Y2), дополненных третьим параметром Z1 и Z2, может быть такой:

PS. Наверняка в алгоритме расчета координат стоит предусловие, в котором сравниваются значения функции в вершинах ребра между собой и с заданным значением уровня. Там должно использоваться такое же значение epsilon.

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

К сожалению, не получится. Делится только часть выражения, определяющего координату точки.

Да. Детали:

- ребра надо создавать даже если они не пересечены ф-цией (с указателем на сегмент NULL). Иначе опять точность подгадит ( для одного треугольника пересечение есть, для др нет)

- я предлагал "указатель на точку пересечения", но "указатель/ссылка на сегмент" не хуже, а может и лучше (напр пересечение проходит по 2 точкам одного ребра)

Класс такого рода задач довольно широк, и ф-ция необязательно аналитическая - напр просто "нечто" возвращающее +/-, ну и поехали по ребру делением пополам. А если и есть формулы, то они могут быть слишком сложны. В общем в "подробности ф-ции" лучше не влезать

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

В текущей реализации функция на треугольнике линейна. Т.е. дважды пересекать одно и то же ребро изолиния не может по определению.

Есть ещё такая тема как Интервальная арифметика (http://ru.wikipedia.org/wiki/%D0%98%D0%BD%D1%82%D0%B5%D1%80%D0%B2%D0%B0%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D0%B0%D1%80%D0%B8%D1%84%D0%BC%D0%B5%D1%82%D0%B8%D0%BA%D0%B0).
Она как раз позволяет более точно учесть проблемы с точностью представлений и позгешностью.
Реализация есть в том же Boost-е. (http://www.boost.org/doc/libs/1_48_0/libs/numeric/interval/doc/interval.htm)