Сравнение вещественных чисел

[Blackwanderer]

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

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

Edit: не хотите пока заводить новые структуры данных - ладно.

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

[Igors]

теоретически эти точки совпадают, на практике их координаты не совпадают абсолютно точно. Если взять слишком маленькую 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 это не решает  

[Blackwanderer]

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

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

[Igors]

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

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

[popper]

Причем линия обязательно должна быть непрерывной.

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

Код:

if (fabs(Z1 - Z2) < epsilon) {
  Z = fabs(Z1); // Z = fabs( min (Z1, Z2) );
  if (fabs(X1 - X2) < epsilon * Z1 && fabs(Y1 - Y2) < epsilon * Z1)
     //OK!!!
}

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

[Blackwanderer]

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

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

Скорее всего, формула основана на линейной интерполяции, в которой производится деление на величину (z2-z1) (это значения функции в вершинах ребра). Хорошо бы избежать операции деления, т.е. в результаты расчета заносить только числитель, и хранить дополнительно третьим параметром Z этот знаменатель.

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

[Igors]

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

Да. Детали:

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

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

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

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

[Blackwanderer]

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

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

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

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

[Tonal]

Есть ещё такая тема как Интервальная арифметика.
Она как раз позволяет более точно учесть проблемы с точностью представлений и позгешностью.
Реализация есть в том же Boost-е.