Комментарии 42
Какая боль...
А ведь можно было свести задачу к пересечению плоскостей – и всё стало бы просто, линейно (пардон, почти: из-за сферы квадратное уравнение вылезает) и без всей этой жуткой тригонометриии.
Набросок:
Вычисляем декартовы координаты x, y, z концов "отрезка" (можно принять радиус сферы за 1, можно R – неважно)
"Отрезок" заменяем на плоскость, проходящую через его концы и центр сферы – она задаётся линейным уравнением вида a*x+b*y+c*z = 0
Два таких уравнения (для двух отрезков) задают прямую. Дополняем их уравнением сферы x²+y²+z²=R² и находим координаты двух точек пересечения.
Проверяем, попадает ли одна из двух точек пересечения на оба отрезка (она в плоскости отрезка, так что всё сведётся к проверке знаков векторных произведений, если надо – могу расписать подробнее)
Переводим декартовы координаты в эйлеровы или какие нам нужны.
Как-то у Вас всё однозначно. Вопрос же всегда в соотношении точности и скорости. Алгоритмы надо по этому принципу сравнивать, а не "просто"/"сложно".
Любые преобразования крадут точность, да и на скорости положительно не скажутся.
Если исходить из этих характеристик, то точность того что я описал (поска точки пересечения) будет близка к точности алгоритма Винсети, а скорость... ну ассимптотическая сложность логарифмическая
Да, у меня всё однозначно: тригонометрические уравнения – всегда боль, если можно их избежать, всего лишь увеличив число переменных на 1 или 2 – следует это сделать. Посмотрите, к примеру, на 3d-графику: все представляют поворот матрицами или кватернионами, никто не использует углы в вычислениях.
Правда, это всё про решение для шара. Для эллипсоида возникнут сложности – на нём, кажется, ортодромия может не лежать в одной плоскости. И, кстати, возможно, это сломает и ваш алгоритм (требует проверки) – не факт, что можно использовать прямую геодезическую задачу для вычисления координат промежуточной точки: не факт, что угол для промежуточной точки будет таким же, как угол для конечной.
Вы какие-то странные вещи пишете про ортодромию и геодезические задачи.
Мне кажется факт что ортодромии в одной плоскости НЕ лежат. Из этого и исходим... но ничего не сломалось
Прямую геодезическую задачу очень даже можно использовать для вычисления координат это факт - в этом смысл задачи.
То что угол для промежуточной и конечной НЕ совпадёт, так это тоже факт, но что это меняет?!
Используя обратную геодезическую задачу, определим расстояние между двумя известными точками и угол.
Полученное расстояние разделим пополам, и, применив прямую геодезическую задачу, по найденному ранее углу, и полученному расстоянию, найдём координаты средней точки.
Так вот, если угол для промежуточной точки не тот же, что для конечной – то нельзя в пункте 2 использовать угол из пункта 1. В этом сложность.
(надо проверять, может, для эллипсоида угол и тот же – я не знаю).
Так здесь в обоих случаях же считаем от самой первой точки, только один раз мы берём полную дистанцию (т.е. вычисляем её), а во второй раз 1/2 прикладываем.
Если формулировка вас запутала, предложите вариант - я исправлю
П.С. Если на следующей итерации мы возьмём среднюю точку и конечную - то углы считаем заново
Да нет, формулировка как раз не запутала. Но, кажется, тут действительно нет проблемы, "азимут" будет одним и тем же. Теперь надо проверять алгоритм решения прямой геодезической задачи :-)
Ещё хотел написать: для эллипсоида можно использовать ваш же подход (проверив его) с алгоритмом половинного деления (или лучше заменить его на метод Ньютона), но всё, что только возможно, вычислять в декартовых координатах.
Upd: И да, любопытно, какова будет погрешность вычисления координат точки пересечения, если её вычислять, пренебрегая "сплющенностью" Земли. Возможно, если точности не хватит – выгодней будет вначале найти точку пересечения "на сфере", а потом скорректировать.
Но ведь и так применяется метод половинного деления: он же, насколько я помню, заключается в том, чтобы искать тот подинтревал где функция меняет знак.
Теперь попарно сравним
с
и
с
ЕСЛИ у
и
разные знаки, ТО
ЕСЛИ у
и
разные знаки, ТО
Различие знаков говорит о том, что пересечение произошло на этом участке, между соответствующими значениями долготы. Алгоритм рекурсивно продолжается, переходим на ШАГ 2.
Вспомнился мультик: "Лучше день потерять… Запомни: лучше день потерять, потом за пять минут долететь."
Так и тут, лучше преобразования сначала и в конце сделать, зато все основные вычисления без тригонометрии в векторах делать, точнее будет.
Что вы несёте? На чём вообще основано это утверждение что тригонометрия снижает как-то точность?
На опыте и оценках погрешностей вычислений.
Тогда от Филдсовской премии вас отделяет один шаг. Доказательство того, что тригонометрические функции врут - перевернёт мир вычислений
Врут не тригонометрические функции, а численные вычисления с ограниченной точностью.
Как при вычислении корней квадратного уравнения, есть абсолютно правильная школьная формула. Но знакомые с численными вычислениями знают что для точных вычислений надо, вдобавок, использовать симметричную с дискриминантом в знаменателе.
А какая область применения?
Для какой задачи важна эллиптичность Земли? И какая ошибка будет если ей пренебречь?
Мы использовали для триангуляции, применяли жадный алгоритм триангуляции, а он основан на поиске пересечений.
Если грубо, то я хотел чтобы алгоритм разбивал на треугольники как в первом примере (слева), а он разбивал как во втором (справа).
Это ещё красивый пример, на самом деле всё было гораздо хуже, получались треугольники с углами 170° , 7°, 3°, и высота у них была несколько метров при длине основания в несколько километров.
Тогда решили всеми правдами и неправдами повысить точность.
П.С. от проблемы полностью не избавились, но триангуляция стала строиться иначе
П.П.С. мы в итоге остановились на упрощении до сферы, но наработки остались
Точность от всей этой тригонометрии сильно страдает. В статье есть оптимистическая фраза
можно использовать координаты любой из двух пересекающихся линий - значение совпадёт
Они не совпадут конечно, погрешность может быть существенная, особенно для отрезков около полюсов, почти-параллельных отрезков и так далее.
Для вычисления на сфере есть гугловская библиотека S2 geometry, рекомендую глянуть как они это считают. Там оценивается точность, выбирается порядок вычислений дающий меньшую погрешность, и иногда все равно вынуждено использующую exact float для получения достаточно точного результата. Правда, она только для правильной сферы.
Совпадут на сфере, решение же аналитическое, если не совпали - значит ошибка в расчётах. С чего вообще в тригонометрии точность страдает?
Вы конечно можете сказать что если 10 разделить на 3, а потом умножить на 3, то 10 мы уже не получим, - да, из-за подобных выкрутасов с плавающей точкой не совпадут, но я бы не называл эту погрешность существенной.
П.С. я не понял что значит "правильная сфера".
да, из-за подобных выкрутасов с плавающей точкой не совпадут, но я бы не называл эту погрешность существенной.
У меня требования когда эта погрешность бывает очень существенной. Например, из-за такой погрешности в вычислениях можно неправильно определить с какой стороны от отрезка лежит точка.
Тогда у меня для Вас плохие новости, вы только посмотрите сколько тригонометрии в преобразовании координат: ГОСТ 32453—2017
Можно очень аккуратно использовать тригонометрию, рассчитывая и указывая как в ГОСТ получаемую точность вычислений. Главное не считать это "аналитической" формулой, и понимать что не любая правильная с аналитической точки зрения формула даст такую же точность в любой точке глобуса.
"Аналитическая" формула для пересечения на сфере у вас падает с делением на ноль если один из отрезков идёт строго вдоль меридиана, или если скажем результат попал на 90 градусов долготы. Уже это должно наводить на мысли что и рядом с этими точками точность вычислений сильно страдает.
Вектора обычно лучше хотя бы тем, что избегают таких выделенных значений.
Не падают, для этого есть проверки. Получается что пишем if и вектора теряют преимущества?
Так нет проверок, ни здесь, ни на GitHub. Но попробуйте добавить проверки, попробуйте, очень интересно будет посмотреть.
Выше же был намек - просто добавить в SpheroidLongitude проверку и специальный случай для coord12.LonR == coord11.LonR и подобных случаев не поможет, численная неустойчивость очевидно останется.
Заколебал уже. Мне посыл: "давайте избавляться от тригонометрии ибо неточно, а для этого завезём ещё больше тригонометрии на преобразование координат, но аккуратно" - непонятен.
Проверки то дело житейское, всегда можно где-то улучшить, было бы что, а вот сама логика использовать тригонометрию для избавления от тригонометрии немного хромает.
Ну отличный метод, вместо того чтобы учиться чему-то - минусовать карму пытающимся объяснить основы численных методов и показать очевидные ошибки. Это лучший способ избавиться от численной неустойчивости формул :)
Вы не пытаетесь объяснять, вы делаете набросы, так что метод норм
Особенно приятно услышать про "набросы" от человека начинающего общение с "Что вы несёте?", когда он сам не понял простейших вещей.
Ну да ладно, оставайтесь с "аналитической" формулой, для которой за минуту легко придумать очевидный пример, где она выдаст ошибку на десятки градусов широты :).
я так понял "очевидного примера" не последует?!
UPD (немного проясню свою позицию, видимо вы меня не понимаете):
Вы пишете про "очевидные примеры", пишете как точнее, а как нет, но это строго ИМХО. Экстраординарные утверждения требуют экстраординарных доказательств, иначе это наброс, поэтому не удивляйтесь тому что услышьте в ответ "Что вы несёте?".
Я Вам пишу про то, что есть явное противоречие, потому что трансформация сама требует тригонометрии, а Вы отвечаете про "аккуратно можно", а потом говорите что я вам несправедливо ставлю минус.
Прошу прощения если кажусь резким, но полагаю что минус был заслужен.
И не исключаю что моя реализация может содержать ошибки, буду значит стараться "аккуратно" использовать тригонометрия, а не избавляться от неё.
Это для вас это были экстраординарные утверждения, для людей сколько-нибудь работавших с численными методами на практике - очевидные. Поэтому первый совет и был - учиться у умных людей (было предложено учиться у Гугла, там умные люди, - была упомянута S2 geometry, в ней подробные комментарии когда какой метод вычислений, когда и почему выбран, и оценки ошибок). Писать лекцию в комментариях я не собирался и после откровенного хамства и минусования кармы точно не буду.
Пример очевиден, выше дважды описано как его сконструировать, и ошибка в десятки градусов проверена на коде из GitHub.
Но нам то нужен эллипсоид, ибо Земля сплюснута.
Если точнее, то Земля - геоид, а не эллипсоид.
Земля то конечно геоид, но математической формулой геоид не описать, поэтому эллипсоид принят как наиболее близкий к поверхности геоида
Почему это не описать? Если вы не знаете такой формулы - это не значит невозможно, нет никаких причин для этого.
А кто-то смог? Можно как-то ознакомиться хоть.
А причины есть, - неравномерное распределение масс Земли. Измерить геоид в каждой конкретной точке возможно, и возможно построить 3D модель. Можно даже сделать математическую формулу которая будет лишь апроксимацией для конкретного региона, но не более.
П.С. да и зачем, мы же координаты считаем а не гравитационный потенцииал?
Не знаю, геодезией не интересуюсь. Но знаю, что кватернионы в том числе и для подобного рода задач были придуманы, и для множества достаточно простых моделей готовых формул в сети не найти. В 2D, например, я самостоятельно выводил формулы в комплексных числах для механизма Чебышева и промежуточных фигур между окружностью и правильными многоугольниками. Тоже вроде всем давно хорошо известные модели, а готовых формул не нашлось.
Ну геоид - "подобный Земле", это по сути ни о чем не говорит.
Если хочется просто и элегантно, нужно сводить к эллипсоиду или уж тогда просто фигячить модель на реальных данных с заданной точностью.
И да, я его придумал потому что не смог нагуглить, может он где‑то и без меня описан был, может за эти 7 лет кто‑то написал что‑то похожее, а может я придумал какую‑то фигню, которую умные люди изобретать не станут...
Вот есть статья 2017 года про алгоритм для поиска пересечения двух ортодромий: https://www.researchgate.net/publication/321358300_Intersection_and_point-to-line_solutions_for_geodesics_on_the_ellipsoid
Спасибо.
В своё оправдание скажу что я выложил своё решение в июне 2017 года (и это уже готовый результат, понятно что потом ещё приводился в какой-то более менее приемлемый вид):
А этот документ от 29 марта 2018 (если это не перезалив):
Я не к тому что они скопипастили, подготовка публикации тоже требует времени большого, а к тому что загуглить именно эту я не мог.
Но в целом да, они описывают тот же самый итеративный подход, с той разницей, что я для решения геодезических задач применяю алгоритмы Винсенти, а они сводят всё к сфере и применяют Cotangent four-part formulae чтобы это ни значило.
(where we will use a suitable value for the radius of the sphere R as e.g. the semimajor axis of the ellipsoid)
...
Now we resort to a rarely used formula of spherical trigonometry known as cotangent four-part formula
На мой взгляд это странно, если сводить всё к сфере - то итерационный алгоритм не нужен.
Простите, не смог внимательно прочитать вашу статью, но, если не ошибаюсь, вы пропустили одно ключевое свойство:
Rsin(A)=cte, где А - азимут, а R - радиус параллели.
В точке пересечения параллель будет одна, и, следовательно, её радиус тоже. Азимут вычисляется из прямой геодезической проблемы многими приложениями, отсюда же значение константы (радиус начальной параллели нам известен).
А начиная с каких расстояний имеет смысл начинать использовать ортодромию вместо плоских прямых? На расстонии пары километров ошибка незначительна, а когда она накапливается достаточно?
Ну как вариант можем опираться на вычисление длин параллелей.
Например для сферы (без учёта эллиптического сжатия, для упрощения), если я ни где не ошибся в расчётах, то:
45 параллель имеет длину 28337 км
46 параллель - 27838 км
т.е. разница 499 км, или примерно 1.4 км на каждый градус отложенный вдоль параллели.
В прямоугольных координатах эта разница теряется, тут уже сами решайте готовы ли вы к такой погрешности в 1,4 км на каждый градус.
При том чем выше - тем больше погрешность, если между 45 и 46 параллелями различие в длине 499 км, то между 46 и 47 уже - 507, соответственно погрешность будет нарастать.
Ну и компенсировать погрешность как-то сложно, потому что, если длинны широт различаются, то длины меридианов - нет, и прямоугольные координаты не делают таких различий.
Изобретаю свой сложный способ поиска координат точки пересечения двух линий