Сглаживание и интерполяция сплайнами

[Unling]

Всем привет. У меня имеется набор точек (траектория) в трёхмерном пространстве и я хочу её умно сгладить. Точки эти в пространстве распределены сильно неравномерно, на прямых участках между соседними точками может быть 500 м, а на поворотах между соседними точками 20 см. Расскажу, что уже попробовал:
1) kernel_smoothing. Не подошёл, т.к. либо сжирает углы на поворотах, либо обеспеичвает недостаточное сглаживание на длинных прямых.
2) Наивный следующий подход: прохожу по траектории скользящим окном размера +-10 точек от центральной точки окна. Внутри окна строю по точкам полином второго порядка. И по этому полиному рассчитываю набор точек около центра окна. Затем сдвигаю окно и повторяю процедуру. Это мне подходит, но это работает слишком медленно. Очевидно, можно как-то экономить при вычислениях на соседних окнах. Вопрос, как именно.
3) Фильтр Савицкого-Голая из scipy. Он в целом делает что-то похожее на наивный подход п.2, но не делает интерполяцию, а она мне нужна.
4) интерполяция сплайнами, типа между i и i+1 точками кривая задаётся как полином 2 порядка и гарантированно проходит через все входные точки. Но мне наоборот это не нужно. Порядок моего полинома очевидно должен быть меньше размера окна.
На картинке нарисовал пример входных данных и что хочется получить:

Посоветуйте, пожалуйста, куда копать.

[StaticZero]

Сплайн степени 3 дефекта 1 пробовали? По-моему очень хорошо бы заехал сюда.

[Unling]

StaticZero, мм, а можно поодробнее, что вы имеете в виду? На всякий случай поясню, на картинке только маленький кусочек траектории, на самом деле в ней множество прямых участков и поворотов. Если есть какие-то ссылки на код python или С++, это было бы хорошо

[StaticZero]

Нет, извините, код вы напишете сами, однако соображения следующие. (Я туманно написал, это правда)

Будем задавать траекторию в виде $\mathbf r = \mathbf r(t)$, где $t$ --- это такой параметр. Будем считать, что точки у вас параметрические $t_0, t_1, \ldots, t_N$, а вы знаете в них значения $\mathbf r(t_i) = \mathbf r_i$. Попробуем построить пространственный сплайн покоординатным методом. Для этого мы рассмотрим отдельно, скажем, координату $z$.

Обозначим кубический полином, центрированный на $t_i$ и проводимый через соседние две точки $t_{i - 1}$ и $t_{i + 1}$, как $P_i$; он задается, как

$$ P_i(t) = p_{i0} + p_{i1} (t - t_i) + p_{i2} (t - t_i)^2 + p_{i3} (t - t_i)^3. $$

Обозначим $t_{i + 1} = c_i$, $t_{i - 1} = a_i$ (концы отрезка, соответствующие $i$-му полиному).

Наложим условия непрерывности на концах:

$$ P_i(a_i) = P_{i - 1}(c_{i - 1}), \qquad P_i(c_i) = P_{i + 1}(a_{i + 1}). $$

Такие же условия на производные -- ещё одно уравнение. Такие же условия на вторые производные (сиречь кривизны) -- третье уравнение.
В точке $t_i$ $P_i(t_i) = P_{i0} = z(t_i) = z_i$, это четвёртое уравнение. Если сплайн, скажем, состоит из двух полиномов, то
в их общей точке три уравнения непрерывности это общие для них уравнения, условие интерполяции в ней это четвёртое, плюс ещё по два для внутренней точки каждого полинома, всего шесть уравнений (а коэффициентов восемь), так что придётся дописать ещё два уравнения для концов,
чтобы замкнуть систему уравнений для коэффициентов сплайна.

Это каноническое построение сплайна $S^2_3$ вы используете для движения в поворотах. Когда вы выходите на длинный прямой участок, вам нужна аппроксимация, которая:
1) обеспечивает непрерывность между конечными точками сплайна;
2) обеспечивает непрерывность первых производных (скоростей по сути) (если вы аппроксимируете траекторию механического тела, то это естественное требование)
3) собственно, аппроксимирует данные.

Условие непрерывности на каждом конце два уравнения, непрерывности производных ещё по два, и остальные степени свободы в вашем распоряжении. Но непрерывность самой функции со сплайном вы, допустим, запишете легко (там понятно, чему должны быть равны сами функции).

Как быть с непрерывностью производных, ведь сплайн ещё не замкнут? С одной стороны, уравнение на непрерывность производных с обоих концов сплайна вы можете записать, тогда сплайн замкнётся. С другой стороны, сами значения производных вам станут известны только тогда, когда вы сделаете аппроксимацию. Один из возможных путей --- вы её сначала делаете, потом записываете на листочек значения производной в концах (допустим, $5{,}48$ и $-9{,}76$, потом накладываете это в качестве условия $\dot S(\text{нужный конец}) = 5{,}48$, $\dot S(\text{другой конец}) = -9{,}76$) ну и т. д. Но это не единственный метод.

Подумайте об этом сначала; код это меньшая из проблем.

[Unling]

StaticZero в сообщении #1473329 писал(а):

Когда вы выходите на длинный прямой участок, вам нужна аппроксимация, которая:

Спасибо! Правильно ли я понял, что в первой части вашего ответа описана процедура, которая построит сплайн, проходящий через все заданные точки? А для прямых участков нужно модифицировать её (ослабить, как бы), чтобы не требовать чёткого прохождения всех точек? Если так, нужно как-то детектировать поворот/не поворот. Сомнительно, что это можно сделать надёжно. Поправьте, пожалуйста, если так. Для прямых нужно заменить условия прохождения через точки на какие-то выражения типа суммы квадратов отклонений?

[StaticZero]

Unling в сообщении #1473332 писал(а):

Правильно ли я понял, что в первой части вашего ответа описана процедура, которая построит сплайн, проходящий через все заданные точки? А для прямых участков нужно модифицировать её (ослабить, как бы), чтобы не требовать чёткого прохождения всех точек?

Понимаете на первый взгляд правильно. Но на второй взгляд вам нужно очень чётко осознавать, что аппроксимация и интерполяция это вещи весьма разные.

Unling в сообщении #1473332 писал(а):

Если так, нужно как-то детектировать поворот/не поворот.

Ну возьмите в качестве условия густоту точек на единицу длины шкалы параметра $t$.

Unling в сообщении #1473332 писал(а):

Для прямых нужно заменить условия прохождения через точки на какие-то выражения типа суммы квадратов отклонений?

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

[Unling]

StaticZero в сообщении #1473334 писал(а):

Ну возьмите в качестве условия густоту точек на единицу длины шкалы параметра $t$.
Вам искусственным образом надо стыковать результаты.

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

[StaticZero]

Unling в сообщении #1473353 писал(а):

По поводу густоты точек - тенденция, что на поворотах точек больше, она есть, но это не всегда так.

Смысл в том, чтобы не получить жёсткую задачу. Вы говорите, что на прямых участках порядка 500 м может быть, а на поворотах порядка 0,2 м. Разница на три порядка это не очень страшно, но строить сплайны по точкам с одним порядком расстояний, а аппроксимацию -- по точкам с другим порядком это правильно.

Unling в сообщении #1473353 писал(а):

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

Посмотрите описание в википедии; повторите его своими руками, тогда точно разберётесь. Инструмент полезный.