Ошибка в программе из-за вычислений с погрешностью

[Sonic86]

Я написал программу. Программа принимает по одной некоторое множество точек и строит по ней триангуляцию (добавляет новую точку - достраивает триангуляцию и т.д.). По некоторым причинам я обрабатываю координаты точек с точностью $\epsilon$. В частности, если уже в триангуляции есть отрезок из точек $P_0,P_1$ и мы проверяем для текущей точки $P$ верно ли, что $P \in [P_0P_1]$, то мы проверяем это не точно, а с погрешностью $\epsilon$. При этом в работе программы обнаружена такая ошибка. Случается так, что уже есть отрезок $[P_0P_1]$, мы берем новую точку $P_2$, она лежит $\epsilon$-близко к $[P_0P_1]$ и тогда мы считаем, что $P_2 \in [P_0P_1]$. При этом на самом деле прямые $P_2P_0$ и $P_2P_1$ уже непараллельны и это отражается Далее, берем еще одну точку $P_3$ $\epsilon$-близкую к $[P_2P_0]$. И тоже после проверки считаем, что $P_3$ лежит на $[P_2P_0]$. Но при этом может оказаться, что $P_3$ не является $\epsilon$-близкой к $[P_2P_1]$. И тогда триангуляция перестраивается и разрушается, сначала в ней появляются одинаковые треугольники, потом треугольники с совпадающими вершинами, а потом начинаются страшные глюки.
Вопрос в том, что делать. Должны же быть какие-то стандартные общие рецепты связанные с машинной точностью вычислений. Я их просто не знаю. Посоветуйте книжку или сайт или еще что-то.

[creative]

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

[Circiter]

2Sonic86
Слишком странные симптомы для численных нестабильностей. Возможно, проблема кроется в самом алгоритме. Попробуйте вбить в поисковик что-нибудь вроде incremental triangualtion.

-- Ср авг 04, 2010 03:14:08 --

Вот даже в википедии есть немного информации об on-line разбиении Делоне + куча ссылок.

[Yuri Gendelman]

Sonic86,

Обычно проблемами, связанными с машинной точностью вычислений, называют проблемы неустойчивости в численных методах решения уравнений, нахождения СЗ и СВ, пр. Ваша проблема лишь косвенно связана с точностью вычислений. Дело в том, что при любом сколь угодно малом $\epsilon$ > 0 используемое Вами определение близости нетранзитивно: если точка $Q_0$ близка к $Q_1$, а $Q_1$ близка к $Q_2$, то это не значит, что $Q_0$ близка к $Q_2$.

Совет: если обнаруживается, что точка $P_2$ $\epsilon$-близка к отрезку $[P_0P_1]$, то отбросьте отрезок $[P_0P_1]$ и добавьте отрезки $[P_0P_2]$ и $[P_2P_1]$.
Вариант: оставьте отрезок $[P_0P_1]$ и просто отбросьте точку $P_2$.

[Sonic86]

creative писал(а):

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

Если будет сильно плохо, выложу. Сейчас я все еще нахожу ошибки.

circiter писал(а):

2Sonic86
Слишком странные симптомы для численных нестабильностей. Возможно, проблема кроется в самом алгоритме. Попробуйте вбить в поисковик что-нибудь вроде incremental triangualtion.
Вот даже в википедии есть немного информации об on-line разбиении Делоне + куча ссылок.

Алгоритм я брал с хорошей книжки. Там было и про численные ошибки, но не про это. Впрочем, еще посмотрю.
Блин, англоязычные сайты-то я и не приметил. Посмотрю.

Yuri Gendelman писал(а):

Sonic86,

Обычно проблемами, связанными с машинной точностью вычислений, называют проблемы неустойчивости в численных методах решения уравнений, нахождения СЗ и СВ, пр. Ваша проблема лишь косвенно связана с точностью вычислений. Дело в том, что при любом сколь угодно малом $\epsilon$ > 0 используемое Вами определение близости нетранзитивно: если точка $Q_0$ близка к $Q_1$, а $Q_1$ близка к $Q_2$, то это не значит, что $Q_0$ близка к $Q_2$.

Совет: если обнаруживается, что точка $P_2$ $\epsilon$-близка к отрезку $[P_0P_1]$, то отбросьте отрезок $[P_0P_1]$ и добавьте отрезки $[P_0P_2]$ и $[P_2P_1]$.
Вариант: оставьте отрезок $[P_0P_1]$ и просто отбросьте точку $P_2$.

Думаете, что это не из-за точности?... Я понимаю, что нетранзитивно и в данном конкретном случае способ есть - взять $\epsilon$ меньше (и тогда отношение с большей вероятностью будет транзитивно) или вообще взять $\epsilon=0$...
Просто не хочу уходить от проблемы за счет всяких обстоятельств. Думаю, что потом когда-нибудь еще столкнусь с этим.

Т.е. проблема в алгоритме все-таки. Попробую посмотреть по источникам, если не получится - выложу код.

[Yuri Gendelman]

Sonic86 в сообщении #342487 писал(а):

Думаете, что это не из-за точности?... Я понимаю, что нетранзитивно и в данном конкретном случае способ есть - взять $\epsilon$ меньше (и тогда отношение с большей вероятностью будет транзитивно) или вообще взять $\epsilon=0$...

Я не говорил, что не из-за точности. Но при любом сколь угодно малом $\epsilon>0$ ситуацию, о которой Вы писали, нужно в алгоритме предусматривать. В численных методах такая ситуация - скорее исключение.

И даже если Вы будете считать, что $\epsilon=0$, на самом деле $\epsilon$ будет больше 0. В программах Уилкинсона (для СЗ и СВ), например, используется машинно-зависимая оценка точности представления чисел: такое $EPS>0$, что $1+EPS > 1$.

Разумеется, чем меньше $\epsilon$, тем с большей вероятностью Ваш алгоритм будет правильно работать. Вас устроит, если Ваша программа будет работать почти всегда? Кроме того, чем меньше $\epsilon$, тем выше вероятность, что при триангуляции у Вас появятся артефакты - треугольники с почти нулевой площадью. Я бы все-таки советовал скорректировать алгоритм.

[Sonic86]

Yuri Gendelman писал(а):

Разумеется, чем меньше $\epsilon$, тем с большей вероятностью Ваш алгоритм будет правильно работать. Вас устроит, если Ваша программа будет работать почти всегда? Кроме того, чем меньше $\epsilon$, тем выше вероятность, что при триангуляции у Вас появятся артефакты - треугольники с почти нулевой площадью. Я бы все-таки советовал скорректировать алгоритм.

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

[Circiter]

И все-таки попробуйте триангуляцию Делоне. Она гарантирует, что треугольники будут достаточно "жирными". Там в вики и на русском статья есть...

[sergey83]

И при такой триангуляции (Делоне) к вам не будут приходить точки, почти лежащие на отрезке. Ваш набор точек никуда не денется, но будут строится другие треугольники.

[Sonic86]

circiter писал(а):

И все-таки попробуйте триангуляцию Делоне. Она гарантирует, что треугольники будут достаточно "жирными". Там в вики и на русском статья есть...

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

sergey83 писал(а):

И при такой триангуляции (Делоне) к вам не будут приходить точки, почти лежащие на отрезке. Ваш набор точек никуда не денется, но будут строится другие треугольники.

Не понимаю. Берем точки

$$P_0(0;0), P_1(0;100), P_2(100;0), P_3(100;100), P_4(50;50), P_5(52;50), P_6(51;50+\delta)$$

и все - вот и длинный тонкий треугольник. К тому же дело-то в том, что у меня триангуляция не появляется сразу готовая, а постепенно строится, перестраивается при добавлении каждой новой точки. М.б. после каждой перестройки триангуляции получаются только хорошие треугольники, но до перестройки - не обязательно.

Sonic86 писал(а):

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

Код выложить не смогу после того, как я оттуда выкинул все лишнее, все равно осталось 500 строк. К тому же его еще протестить надо - вдруг я его неправильно исправил, а я уже 4 часа убил
М.б. все-таки попробуем отвлеченно проблему решить. Например, была такая идея - записывать точки $P_0,P_1,...$ для которых $\rho(P_{n-1}P_n,P_{n+1})<\epsilon$ в специальные массивы и предполагать там, что они все лежат на одной прямой (т.е. считать отношение транзитивным). Кода станет еще больше, но сильнее тормозить не будет - случай все равно редкий, массивы не будут много занимать. И вроде бы таким способом сильно далекие точки в этот массив не будут включены

[Circiter]

2Sonic86

Цитата:

У меня и так Делоне.

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

Проблема в том, что on-line вариант этого алгоритма требует небольшой модификации уже построенной триангуляции при поступлении новой точки и эта перестройка обходится примерно в $\mathcal{O}(\log n)$, а при реализации "в лоб" так и вообще в $\mathcal{O}(n)$ (определяется и фрагментируется треугольник, в который попала точка, после чего применяется т.н. flipping technique). Но если этого не делать, то это будет уже не Делоне-триангуляция.

Использование этой техники подразумевает задействования хитрых структур данных. Вопросы численной стабильности возникают только в сложных случаях. Впрочем, работ на эту тему немало, например S.Fortune, Numerical Stability of Algorithms for 2D Delaunay Triangulations.

[Sonic86]

circiter писал(а):

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

Ув. circiter, у меня действительно используется триангуляция Делоне, я сам не стал ничего придумывать, поскольку вовремя заметил, что в книге Скворцова "Триангуляция Делоне и ее применение" есть гораздо больше материала, чем я сам выдумал. Я взял все из книги. Выбрал самый лучший по мнению автора алгоритм и реализовал. И на случайных данных алгоритм прекрасно работает. И я вполне понимаю, почему он так хорошо работает (асимптотику).
Единственное, в чем я сомневаюсь, так это в проверке условия Делоне для треугольников $P_1P_2P_3$ и $P_1P_3P_0$, вместо $\angle P_1P_2P_3 + \angle P_1P_0P_3 \leq \pi$ я взял $P_1P_3 \leq P_2P_0$.
А вот насчет "в любой момент времени" - нет. Сразу при вставке новой точки триангуляция не обязана быть оптимальной. А если только вставить точку, без ребер, так и самой триангуляции пока еще тоже нету.
Книжку попробую найти.

[Circiter]

2Sonic86

Цитата:

И я вполне понимаю, почему он так хорошо работает (асимптотику).

Ну я бы не сказал, что это такой уж и монструозный и неповоротливый алгоритм. По-моему он вполне шустрый, даже для real-time сгодится...

Цитата:

Единственное, в чем я сомневаюсь, так это в проверке условия Делоне для треугольников $P_1P_2P_3$ и $P_1P_3P_0$, вместо $\angle P_1P_2P_3 + \angle P_1P_0P_3 \leq \pi$ я взял $P_1P_3 \leq P_2P_0$.

Ну тогда неудивительно. :) Сделайте как в книжке (через скалярные произведения, к примеру), наверняка ситуация сразу-же улучшится.

Цитата:

А вот насчет "в любой момент времени" - нет. Сразу при вставке новой точки триангуляция не обязана быть оптимальной. А если только вставить точку, без ребер, так и самой триангуляции пока еще тоже нету.

Это понятно. :) Я имел ввиду не столь тривиальные случаи.

В общем, мне самому стало интересно. Попробую реализовать свой простейший Делоне-триангулятор и погонять его на плохих данных. О результатах сообщу.

[Sonic86]

circiter писал(а):

Ну я бы не сказал, что это такой уж и монструозный и неповоротливый алгоритм. По-моему он вполне шустрый, даже для real-time сгодится...

Я не иронизировал вроде. С++-ный вариант алгоритма 100 точек менее чем за секунду обрабатывал, а моя прога тормозит из-за того, что работает с БД, а не с массивами.

[Sonic86]

circiter, Вы меня натолкнули на проверку, я нашел несколько глупых ошибок и исправил. Сейчас проверяю вычисление с погрешностью, возможно, что и сработает.