Как не нужно вычислять центр масс(центроид).

[sdikiy]

Привет,

недавно на просторах интернет обнаружил "замечательный" алгоритм вычисления площади и центроида полигона.
http://paulbourke.net/geometry/polyarea/
С кучей примеров на различных языках программирования.

Эта штука смогла "просочиться" в несколько проектов ГИС:
https://github.com/openlayers/openlayer ... ng.js#L193
https://github.com/qgis/Quantum-GIS/blo ... r.cpp#L224

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

[Игорь Белов]

Формула вычисления площади безупречна (почти). Формулы координат центроида безупречны. На плоскости.

P. S. Автор статьи предупреждает, что формула площади не работает для контура с самопересечениями. Действительно, посчитанное "в лоб" значение окажется равным разности площадей подконтуров, пройденных против часовой стрелки, и подконтуров, пройденных по часовой (в традиционной системе координат).

P. P. S. Однако именно это неправильное значение A и следует подставлять в формулы координат центроида.

P. P. P. S. Конечно, если оно не равно нулю

[sdikiy]

Подсказка номер раз:

С математической точки зрения, с учетом ограничений указанных автором, формулы безупречны.

Но попробуйте посчитать координаты центроида, в проекции 900913, для туалета типа "сортир" расположенного на Аляске. т.е. думать нужно в направлении погрешностей, вносимых алгоритмом, при его реализации для реальных расчетов.

PS. Как то так http://en.wikipedia.org/wiki/User:Sdiki ... onCentroid

[Игорь Белов]

Спасибо, sdikiy, за подсказку.
С учётом того, что далеко не все посетители форума имеют блестящее образование по теме "Как нужно вычислять центр масс (центроид)", думаю, Вы могли бы принести пользу, если бы написали статью с обзором современных подходов.

[Максим Дубинин]

как обычно в таких случаях говорят "неистово плюсую".

[Игорь Белов]

Простите мою тупость. Всё не могу понять, в чём корень проблемы. Если в недопустимом росте погрешностей, так методы борьбы описаны в горах учебников полувековой давности. В случае "подсказки раз" (большие значения координат и малый экстент полигона) применИм элементарный подход:

• найти центроид-0 как середину экстента;
• перенести начало координат в центроид-0;
• вычислить A без потери точности;
• если A = 0, координаты центроида нули, иначе вычислить их по формулам;
• вернуть начало координат на место.

Математика остаётся математикой.
Может, есть другие подсказки?

[sdikiy]

@ErnieBoyd, огромное спасибо за Ваш пост. Да, проблема алгоритма в огромных погрешностях, в случае (большие значения координат и малый экстент полигона).

Метод борьбы - крайне прост и очевиден... вот пример реализации алгоритма описанного Вами:
https://github.com/sdikiy/openlayers/co ... fab0810bee

Подсказок больше нет, все действительно крайне банально.

Вопрос в другом, почему люди не применяют эти методы, описанные в учебниках полувековой давности?
Почему сопротивляются, когда им явно указывают на ошибки?
http://hub.qgis.org/issues/5256

Возможно их плохо учат в институтах?

[Игорь Белов]

Каждый живёт в своём мирке comfortably numb. До нас непросто достучаться, сочувствую.

[Максим Дубинин]

Почему сопротивляются, когда им явно указывают на ошибки?
http://hub.qgis.org/issues/5256

Хаб у меня висит, но судя по тому, что я вижу в кэше, никто особо не "сопротивлялся", просят того же самого, объяснить, вроде довольно вежливо.

Изложите суть проблемы и решение в виде последовательной четкой аргументации в одном месте (разбросанные ссылки на оракл, википедию, сайты людей, патчи в OL/QGIS это хорошо, но это хаотично), снабдите иллюстрациями, короче говоря - убедите. Если для вас это важно, сделайте это так, чтобы для того, чтобы понять не нужно было сидеть и проверять вашу математику.

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

Ну а мы можем сделать так, что это будет исправлено в QGIS.

[sdikiy]

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

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

[Максим Дубинин]

ну так скооперовались бы с ErnieBoyd и написали вместе, вы суть, он обоснование

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

Нечто подобное вы сейчас и делаете, пишете в форум, поднимаете проблему. Разработчики видят, я вас уверяю. Мое пожелание было только в том чтобы сделать процесс более эффективным. Но если кугис и проблема вас на самом деле не волнуют (тут я до конца не понял, так что ваши действия вроде показывают обратное , то я конечно пас.

[sdikiy]

"Может, что-то в консерватории подправить?" (с) сами знаете кто

https://github.com/mapserver/mapserver/ ... ve.c#L1232

Похоже все тащат этот код в свои проекты...

PS. Кто знает, как в русскоязычной литературе называется этот алгоритм Shoelace formula

[Voltron]

Мда... даже не знаю, что и сказать

[Игорь Белов]

Между тем можно найти и более простые подходы, устраняющие проблему.
В формуле «ботиночного шнурка» площадь полигона A представляется суммой площадей треугольников ΔAi, как в формуле (1).
Предложение 1
Объявим тип промежуточных переменных в кодах на хабе long double вместо double. Даже в Intel'овской ущербной арифметике с 80-битным представлением это даст 63-битную мантиссу. В результате для полигона с характерным размером стороны 1 м неприятности начнутся только при значениях координат порядка 2³², т.е. около четырёх миллиардов, что с лихвой покрывает практические потребности.
Предложение 2
К выражению в скобках прибавим и вычтем Xi Yi. После перегруппировки членов получим выражение (2), в котором считается разность произведений XΔY и YΔX, порядок величины которых не вызывает проблем. Имеем то же векторное произведение, только не двух длинных сторон треугольника, а длинной и короткой.