Поиск готового решения

[Igors]

Добрый день

Есть N точек (x, y, z) с регулярным шагом по времени. Порядок N - сотни, тысячи и больше. Требуется заменить их кубическим сплайном Безье с минимальным числом контрольных точек. Ну т.е. кривой проходящей через все N точек, разумеется с какой-то заданной погрешностью.

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

Нахожу это http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.96.5193&rep=rep1&type=pdf
По названию вроде "точно в масть", но выкладки довольно обширны и как реализовать - хз. Это, кстати, типичный рез-т для источников с "академическим уклоном".

Ладно, нашел это  http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gems/FitCurves.c
Тут все конкретно, причем по своему опыту знаю: джемсовский код работает железно. Но.. как Вам нравится такой код?   К тому же он 30-летней давности. (а как любят кричать "велик!", "костыль!").
Попробовал найти тех кто ссылается на этот алгоритм - ну неск работ в стиле первой ссылки.

И что делать? Понятно что просторы инета безбрежны и всегда можно сказать "плохо искал". Но практически.. ну погуглю еще пару дней, останусь там же. Дальше что? Решение принимать все равно придется. Поэтому спрошу так:

Стоит ли связываться с найденным готовым решением и пытаться тащить его в рабочий код ?

Спасибо

[ssoft]

Решал такую задачу лет 15 назад). Попробую найти материалы.
Пока что в общем виде так:
- Выбираем вид сплайна. Кубический сплайн или Безье 3-го порядка - это всего 4 управляющие точки.
- Далее с помощью метода наименьших квадратов выражаем искомые величины - коэффициенты или контрольные точки.

Метод наименьших квадратов не сложный. Возможно, пока ищу материалы уже решите задачу)).

[Igors]

ssoft, спасибо за желание помочь

- Выбираем вид сплайна. Кубический сплайн или Безье 3-го порядка - это всего 4 управляющие точки.

Однозначно кубик (2 управляющие точки на сегмент), ну или каждая контрольная точка сплайна имеет 2 управляющие.

Метод наименьших квадратов не сложный. Возможно, пока ищу материалы уже решите задачу)).

МНК навскидку не помню, но разберусь, тоже когда-то делал, НО он ведь ищет прямую проходящую "наилучшим образом" через N точек. Ну может и не только прямую, но "заданную" ф-цию.  А для Безье что задавать? Там же с управляющими точками кривые могут быть любые.

А мож все-таки велик? Ну вот "душонка" (как говорил складской работник) у меня такая - все к велику тянется    Да и ковыряться в том кале мамонта ой несладко.

[Igors]

Следующие 2 дня гугления ничего принципиально нового не внесли, нашел с десяток перепечаток того же алгоритма 90 года. Ну думаю, значит это классика, надо его брать. Перевел на плюсы без проблем, все скомпилил и, наконец, попытался разобраться как он работает. Во-первых, он вычисляет вектора по паре первых/последних точек, а значит если в исходных точках есть "шумы" -  будет работать плохо. Ну то ладно, как-то загладить можно. Но "во-вторых" - просто облом. Он ведь только для "одномерной" кривой типа y = f(t), и на мои 3D точки (x, y, z) не распространяется   А сделать 3 кривые - да, можно, но их никак не удается "склеить". В общем, пока сижу "у разбитого корыта"

[ssoft]

Пока не нашел материалов, еще по архивам поищу ...

Он ведь только для "одномерной" кривой типа y = f(t), и на мои 3D точки (x, y, z) не распространяется ...

А почему для одномерной? Ведь в общем виде кривая представляет собой P = f(t), для кубического сплайна P = P0 + P1 * t + P2 * t^2 + P3 * t^3, где t - вещественное число, P - любой вектор состояния, например, координаты {x, y, z}.
Нужно просто определить для типа точек P алгебраические операции "+", "-" между точками и "*", "/" на число. Формулы те же, результат в виде точек.

[m_ax]

Нужно просто определить для типа точек P алгебраические операции "+", "-" между точками и "*", "/" на число.

Да можно проще.. Всё в конечном счёте сводится к СЛАУ на P0, P1, и т.д..

[ssoft]

Нашел давние материалы (см. во вложении). С тех пор навряд ли что-то изменилось)).

[Igors]

Нашел давние материалы (см. во вложении). С тех пор навряд ли что-то изменилось)).

Довольно общие выкладки, до конкретики еще далеко, но все равно спасибо.

Нашел я одну неказистую ссылочку (файл в аттаче), где товарищ Муртаза Хан предлагает реализацию МНК для кубического Безье. См описание в маленьклм pdf. Правда исходники на языке какого-то мат пакета, но понять можно. И очень мне его метод понравился. Во-первых, четко заявляется что исходные точки могут быть любой размерности. Во-вторых, нет никакой завязки на касательные в крайних точках как в алгоритме Шнейдера (гарантированный геморрой). Ну и вообще гораздо проще.

Ладно, решил я провериться. Строю свою кривую (пока всего 2 точки сплайна + 2 управляющие), потом перевожу ее в N точек и скармливаю алгоритму Хана. РАБОТАЕТ! Получаю те же  управляющие точки что и в исходной кривой. Хотелось похвастаться рисунком, но там просто 2 одинаковые кривые.

Это что же получается? А как же Academic Press, Ньютон, Рафсон и десятки статей которые, в принципе, всему этому вторят? Если вышел Муртаза и довольно простой алгеброй все решил - зачем все это нужно? Может его решение как-то ограничено, чего-то не учитывает? Не вижу чего

Правда с последним пунктом ("Fitting Strategy") я пока не разобрался. Ну побил все рекурсивно, и что дальше? Как работает "break-and-fit" не понял

К сожалению, это еще не все, есть еще проблема, напишу позже

[ssoft]

Довольно общие выкладки, до конкретики еще далеко, но все равно спасибо.

Написал пример получения СЛАУ (во вложении).

Обычно вопросы возникают при определении t_i, которые должны обеспечить наилучшую параметризацию.
http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=zvmmf&paperid=773&option_lang=rus
http://www.dissercat.com/content/nailuchshaya-parametrizatsiya-v-zadachakh-priblizheniya-krivykh-i-poverkhnostei

Вывод этих работ по поиску наилучшей параметризации такой:

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

[Racheengel]

А мне одному интересно стало, как же юзеру то с этими 100500 точек работать?
А главное-зачем?

[ssoft]

А мне одному интересно стало, как же юзеру то с этими 100500 точек работать?
А главное-зачем?

Обычно 100500 точек - это какие-то дискретные измерения, которые хотелось бы заменить простой непрерывной аналитической зависимостью.
Задач, где это необходимо большое количество - наука, графические редакторы, CAD, имитация плавного движения по данным измерения и т.п.
100500 точек аппроксимируются и заменяются, например, 4-мя, с которыми пользователь впоследствии и работает.

[Igors]

Написал пример получения СЛАУ (во вложении).

Да, спасибо, понял. Однако Ваш стартовый "знак удивления" не вполне уместен , въехать о чем речь не так уж просто. 

1) Кто такие "фи" которыми прямо-таки засыпан весь pdf?
Это коэффициенты сплайна (для каждого сегмента) которые зависят только параметра t. Напр все данные есть, рисуем кривую

b0(t) = (1-t)^3
b1(t) = 3 (1-t)^2 t,
b2(t) = 3 (1-t) t^2
b3(t) = t^3

bezier(t) = b0(t) * A0  +  b1(t) * A1  +  b2(t) * A2  +  b3(t) * A3
где
A0 и A3 - точки сплайна, кривая проходит через них
A1 и A2 - управляющие точки, перемещая их юзер настраивает кривизну

Меняя t от 0 до 1 с каким-то шагом получаем кривую

2)

Для аппроксимации требуется решить систему линейных алгебраических
уравнений относительно A вида

Тут я довольно долго "лупал глазками" - откуда эта система взялась, в чем ее смысл и кто такие phi(k)? Наконец сообразил что это рез-т решения МНК который принимается без доказательств. А "k" - индексы коэффициентов (b0 (1, 2, 3) выше)

3)

Если требуется, чтобы сплайн начинался в точке 0 p , а заканчивался в точке
Np p , то СЛАУ можно привести к виду

Ну так получится "масло масляное", т.е. первое и второе ур-я обратятся в тождества. Просто A0 и A3 известны, имеем СЛАУ 2x2, т.е. оставляем только 2-е и 3-е ур-е, и в них суммы при A0 и A3 переносим вправо. Получаем именно то что у Хана - оказывается он ничего не изобретал а воспользовался известным в математике решением (о котором я не знал).

Но то "легко критиковать", в целом - отличное объяснение, еще раз спасибо. Обычно ж вообще ни хрена не получиnь, махнет хвостиком (типа "можно проще")  - и все   

Обычно вопросы возникают при определении t_i, которые должны обеспечить наилучшую параметризацию.

Да, возникают, похоже я опять все неправильно понял. Напишу большой пост (с картинками  )

[Igors]

Еще раз как нарисовать кривую

b0(t) = (1-t)^3
b1(t) = 3 (1-t)^2 t,
b2(t) = 3 (1-t) t^2
b3(t) = t^3

bezier(t) = b0(t) * P0  +  b1(t) * P1  +  b2(t) * P2  +  b3(t) * P3
где
P0 и P3 - точки сплайна, кривая проходит через них
P1 и P2 - управляющие точки, перемещая их юзер настраивает кривизну

Меняя t от 0 до 1 с каким-то шагом получаем кривую

Это соответствует первому аттачу. Слева - путь, т.е. кривая на плоскости, справа - зависимости X и Y от времени t (x(t), y(t)). В верхней точке кривой скорость по Y = 0, а скорость по X везде постоянна. Ну такие вектора мы и даем драгать юзеру в левом окне. Хорошо, он взял и увеличил вектор по X (второй аттач). Путь как-то искривился, при этом у(t) остается без изменений, а вот x(t) - уже не движение с постоянной скоростью.

И вот это изменение X юзера может и не устроить, он хочет чтобы скорость по X оставалась постоянной. Очевидно это возможно (третий аттач), но как сделать это (интерактивной) настройкой? Ведь скорость по Y в верней точке уже 0.

Решается так: втихаря двигаем контрольные точки (в левом окне юзеру это не показываем), НО компенсируем это модуляцией t. Т.е. на вход рисования поступает уже не текущее время, а измененное

t = c3 * t^3 + c2 * t^2 + c1 * t + c0

Коэффициенты "с" вычисляются на основании доп параметра "W" (вес) который есть у каждой точки упр-я. Заметим что на 3-м аттаче вектора в правом  окне стали длиннее - это как раз и есть интерактивное упр-е весами. Чем длиннее вектор - тем больше W, это как бы "замедляет время" в окрестности точки (хотя и звучит бредово)

Ну и, естественно, кривые записанные с W != 1 под МНК не ложатся (не восстанавливаются). Так может Academic Press об этом и толковал? (вряд ли они не знали об МНК). В общем, пока не понял, думаю

[ssoft]

Ну и, естественно, кривые записанные с W != 1 под МНК не ложатся (не восстанавливаются).

Ложатся ложатся).

Для МНК абсолютно без разницы какой набор базисных функций будет использован - ряд Тейлора, ряд Эйлера, полиномы Бернштейна, NURBS и т.п.
Не важно с весами кривые или нет.

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

Ну так получится "масло масляное", т.е. первое и второе ур-я обратятся в тождества. Просто A0 и A3 известны ...

Значения A0 и A3 известны - это только для полиномов Бернштейна (Безье, NURBS и иже с ними).
A0=p0 и A3=p3

Для ряда Тейлора, например, будет
A0 = p0
A0 + A1 + A2 + A3 = p3.

Все зависит от вида сплайна, который определяется видом базисной функции "фи".

[Igors]

Если в процессе аппроксимации необходимо найти еще  эти самые веса, то необходимо сформулировать аналитические зависимости, через которые можно эти самые веса вычислить.

Вот они

P0, P1, P2. P3 - точки сплайна. P0, P3 известны, нужно найти управляющие P1, P2
W1, W2 - веса точек P1 и P2, их тоже нужно найти

Модуляция параметра t
c3 * t^3 + c2 * t^2 + c1 * t = user_T
где
 user_T - "исходное" время на сегменте (от 0 до 1)
 t - модулированное время для вычисления коэффициентов Безье

Т.е. сначала решается кубическое ур-е и находится t

Коэффициенты "c" зависят от W1, W2

c1 = 3.0 * W1;
c2 = 3.0 * (-2.0 * W1 + W2);
c3 = 1 + 3.0 * (W1 - W2);

При W1 = 1/3 и W2 = 2/3 получим user_T = t

И вот еще: наше обсуждение игнорирует (или никак не учитывает) эту ссылку
http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gems/FitCurves.c
Которая по сути перепечатывалась много раз. Наивно полагать что авторы не знали общей математики (МНК, метод последовательных приближений и.т.п). Однако мои попытки "разобраться" в этой работе пока безуспешны. Первое на чем я затыкаюсь - это параметризация по длине пути (хорде chord). Ну ладно, пусть t - расстояние по пути (тоже от 0 до 1), как-то посчитали P1 и P2, и.. что мне с этим делать? У меня входные точки (x, y, z)(t) и выходные должны быть "в том же формате". Можете прояснить этот момент?

Спасибо

Update: подправил вычисление t.