Интерполяция гладкой кривой на плоскости

[Sonic86]

У меня есть конечное множество точек $A_j(x_j;y_j)$ некоторой гладкой кривой $L$ на плоскости без самопересечений, возможно замкнутой. Мне надо в программе построить какую-нибудь интерполяционную кривую $M$, проходящую через все точки $A_j$ (изолинии). Как это делать?
Кривые Безье и В-сплайны мне, кажется не подойдут, т.к. там $A_j \not \in M$. Потом подумал строить кубические сплайны вида $y=P_3(x)$ или $x=P_3(y)$ для каждой пары точек $A_j, A_{j+1}$ с условием гладкости, но тогда мне не хватает условий: есть 2 условия $A_j \in M, A_{j+1} \in M$ и еще одно условие совпадения производной в первой точке с производной во второй точке с предыдущей кривой и все - больше нету.
Параметрические кривые $x=x(t),y=y(t)$, где $x(t),y(t)$ - многочлены 2-й или 3-й степени воде бы тоже не подходят - мало условий для определения коэффициентов, приходится брать не 2, а 3 или 4 точки и тогда начинаются проблемы с гладкостью кривой.
Как строить кривые?



Надеюсь, что мне не нужно читать книжки отсюда: topic667.html

[мат-ламер]

Есть такая разновидность сплайнов - NURBS. Может Вам подойдёт?

[ewert]

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

Правда, с двумя оговорками.

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

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

[Sonic86]

Вчера попробовал взять для каждой $j$-й кривой приближение

$$ \left\{ \begin{array}{cc} x_j(t)= a_jt^2+b_jt+c_j \\ y_j(t)=p_jt^2+q_jt+r_j \end{array} $$

(хотя это обычная парабола) с $0 \leq t \leq 1$ и условиями: 1. каждая кривая проходит через граничные точки + первые и вторые производные в граничных точках равны. Но тогда вылезли зависимости $j+1$-х коэффициентов от $j$-х, что меня совсем не радует, а равенство вторых производных дает страшную кубическую зависимость Пробовал выкинуть условие со второй производной и заменить его на условие "вершина параболы находится на серединном перпендикуляре к P_jP_{j+1}". Условие неестественное, я его мотивировал тем, чтобы вершина параболы вдруг не "выскочила" на дуге. Но формулы слишком страшные получились (вообще, для параболы $(x(t),y(t))$, заданной выше, есть простые формулы для ее вершины?)

мат-ламер писал(а):

Есть такая разновидность сплайнов - NURBS. Может Вам подойдёт?

Напрямую не получается - NURBS тоже не проходит через все узловые точки. Или я чего-то не понял?

ewert
я ж так не выживу - мне же это все в коде писать надо...
Ну узлов на каждой кривой вроде бы довольно много - по нескольку десятков узлов кривые точно будут, а м.б. и до сотни узлов дойдет.
Насчет равенства векторных производных надо подумать... Если других вариантов не найду, возьму этот.
Мне бы было просто замечательно, если бы уравнение кривой строилось локально - по паре или хотя бы по конечному числу точек. Это я хоть запрограммировать смогу. А рекуррентности типа $a_{j+1}=f(a_j)$ можно было бы решать методом итераций, пробежавшись несколько раз по кривой (хотя совсем не уверен в вопросах сходимости).

-- Чт июл 15, 2010 08:31:18 --

Можно наверное, знаете, как попробовать? Звенья ломаной $P_jP_{j+1}$ у меня уже есть. Берем, для каждого угла $P_{j-1}P_jP_{j+1}$ строим прямую $h_j$ так, чтобы она образовывала с отрезками$P_{j-1}P_j$ и $P_jP_{j+1}$ одинаковые углы. Тогда прямые $h_j, h_{j+1}$ задают для кривой $l_j$ направления (касательные, производную) в точках $P_j, P_{j+1}$ - уже 4 условия. И тогда наверное Безье можно построить через систему линейных уравнений $2 \times 2$.

-- Чт июл 15, 2010 08:33:59 --

Хотя нет, я так точку $\bar P_2$ переопределяю только.
Ну значит не Безье, а параболу, вроде бы она тогда подойдет...
Да, подходит!!!! Зашибись! Щас проверю на примерах еще...

[ewert]

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

Вчера попробовал взять для каждой $j$-й кривой приближение

$$ \left\{ \begin{array}{cc} x_j(t)= a_jt^2+b_jt+c_j \\ y_j(t)=p_jt^2+q_jt+r_j \end{array} $$

(хотя это обычная парабола) с $0 \leq t \leq 1$ и условиями: 1. каждая кривая проходит через граничные точки + первые и вторые производные в граничных точках равны. Но тогда вылезли зависимости $j+1$-х коэффициентов от $j$-х, что меня совсем не радует,

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

Впрочем, можно и локально обеспечить гладкость (хотя это уже и не будет сплайн по отдельным координатам). Например, задав наклон в каждом узле как наклон прямой, соединяющей соседние с ней точки. И выбрав третью опорную точку Безье для каждого отрезка на пересечении соотв. линий, проведённых через концы. Это тоже быстро (медленнее, чем в предыдущем случае, но порядок объёма вычислений тот же -- первый). Однако более-менее корректно это будет работать только для выпуклых кривых.

[Sonic86]

ewert
А если парабола - это Безье, то как тогда у меня с Безье не получилось, а с параболой получилось... ?

-- Чт июл 15, 2010 08:41:54 --

Ага, ну понятно все, буду кодить... Спасибо!

[Алексей К.]

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

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

"кубические сплайны вида $y=P_3(x)$ или $x=P_3(y)$..."
и
"кубические сплайны .... для каждой пары точек $A_j, A_{j+1}$".

Вам совсем не надо думать на тему $y=P_3(x)$ или $x=P_3(y)$. (А в чью пользу выбирать, когда и так и так годится?)
Сплайн строится сразу по всем точкам: все они вовлечены в построение системы линейных уравнений, которая потом решается методом прогонки.
Сплайн проходит ровно через ВСЕ заданные точки $A_j$. Вам следует только подумать о граничных условиях для крайних точек, и о параметризации: приписать ли значениям $t$ в узлах $t_j=j$ (что считается дурным тоном), или $t_j$ = (накопленная) длина ломаной $A_0A_1A_2\ldots A_j$.

Цитата:

Кривые Безье ... мне, кажется не подойдут, т.к. там $A_j \not \in M$.

Кубический сплайн будет состоять именно из (кубических) кривых Безье. На каждом участке та линейная система Вам даст, по сути, пару контрольных точек $P_j,Q_j$. Каждый кусочек сплайна типа $A_5A_6$ будет вписан в контрольный полигон $A_5P_5Q_5A_6$. Он не будет проходить через точки $P_5,Q_5$, но конечно пройдёт через $A_5,A_6$.

Мне кажется, я, несмотря на жару, не туфту написал, и не ненужное.

[ewert]

Алексей К. в сообщении #339370 писал(а):

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

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

[мат-ламер]

Однако, как предлагает Sonic86 - строить свой сплайн для каждой пары - абсурдно. Если 100 точек, то получается 10000 сплайнов - перебор явный. Можно строить не для всех точек, а по частям, если по смыслу в некотрых точках будет излом.

Цитата:

Кривые Безье и В-сплайны мне, кажется не подойдут

Sonic86. Вы может путаете точки интерполирования и управляющие точки? (Это я про сплайны - про Безье сказать ничего не могу). Интерполяционный сплайн конечно проходит через точки интерполяции.

[ewert]

мат-ламер в сообщении #339402 писал(а):

строить свой сплайн для каждой пары - абсурдно. Если 100 точек, то получается 10000 сплайнов - перебор явный.

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

[Sonic86]

Товарищи! Я не рублю в сплайнах! Не проходили мы их в школе! Поэтому давайте я просто поясню.
Алексей К.
Когда я писал $y=P_3(x), x=P_3(y)$ то имел ввиду интерполяционные многочлены 3-го порядка, не сплайны.
Думал-то я просто - я посмотрел, что кривые Безье не проходят через все точки, которые ее задают и подумал, что они мне не подойдут - у меня всего $n$ точек у ломаной - по паре точек на кривую, а уже для Безье 2-го порядка надо 3 точки. Откуда их взять? У меня их нету. А решать СЛУ для всех точек каждой изолинии - по-моему слишком сложно, а при обновлении точки изолинии - пересчитывать всю изолинию? Еще и метод решения СЛУ не дай бог придется программировать! У меня и так прога тормознутая! А я еще для каждой ломаной из $n$ точек буду решать СЛУ со скоростью $O(n^3)$. А если я еще и кубического Безье возьму - еще сложнее выйдет.

мат-ламер
Я когда говорил о том, что буду просто строить сплайн для каждой пары, я имел ввиду каждую пару точек, задающую отрезок ломаной. Иначе говоря, для ломаной $A_1A_2...A_nA_1$ из $n$ точек я беру $n$ отрезков $A_jA_{j+1}$ (а не $n^2$ отрезков A_iA_j) и для каждой пары строю изолинию.

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

[Алексей К.]

(Sonic86: 'А если я еще и кубического Безье возьму - еще сложнее выйдет' - ой-ой-ой!)

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

А решать СЛУ для всех точек каждой изолинии - по-моему слишком сложно, а при обновлении точки изолинии - пересчитывать всю изолинию? Еще и метод решения СЛУ не дай бог придется программировать! У меня и так прога тормознутая!

Похоже, я, наконец, всё же нарисую куб-сплайн программку на PostScripte (больше ни на чём не умею). Давно хотелось, да всё как-то лениво. Вот сяду в викэнд на велик, заеду на речку, возьму с собой шпаргалку про метод прогонки, два листочечка бумажки, карандашик, и нарисую...
И пусть вашей тормознутой проге будет стыдно!

Уж на ваших-то языках, на Javaх и Сях, эти штуки одним кликом гуглятся!

[Sonic86]

(Оффтоп)

Алексей К.
Ну если получится буду рад и со стыда прогу перепишу.
Так! Я погуглил. Оказывается метод прогонки проще метода Гаусса. Ненулевых коэффициентов всего 3 в каждой строке.
Тааак..., тогда надо подумать... тогда вполне возможно, что сделаю полноценные сплайны, может быть даже кубические...
Сейчас пока просто немного другие проблемы...

[ewert]

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

Оказывается метод прогонки проще метода Гаусса.

Это ровно метод Гаусса (с обратным ходом) и есть. Просто из-за трёхдиагональности объём вычислений оказывается очень мал.

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