Поиск в гугле на тему «интерполяция замкнутых кривых» не дает особых результатов, кроме как на автореферат из диссертации от 2010 года господина Чашникова Н.В. Именно эту работу я взял за теоретическую основу. Должен сказать, что в процессе анализа математики возникали трудности, которые помог решить сам Николай Викторович, за что ему огромное спасибо! Собственно, дальше будет разбор этого автореферата с конкретными кусками кода.
По поводу интерполяции замкнутых кривых: ещё можно использовать периодические кривые Безье, они рассматривались на нашем семинаре. Там результат получается более традиционный, в виде тригонометрического полинома, то есть обычной функции вещественного аргумента. В некоторых случаях это может быть удобнее, чем дискретные сплайны (которые задаются функцией целочисленного аргумента).
Итак, задача стоит так: на входе есть набор двумерных точек, необходимо построить интерполяционную замкнутую кривую.
Для решения этой задачи можно использовать дискретные N-периодические сплайны с векторными коэффициентами.
Далее я буду использовать следующую терминологию и обозначения:
Полюса сплайна — набор исходных точек (векторов), m — количество полюсов сплайна, n — количество узлов сплайна между двумя соседними полюсами, N = n * m — период сплайна, r — порядок сплайна.
B-сплайн первого порядка вводится вполне явно:
/// <summary>
/// Вычисляет коэфициенты дискретного периодического Q-сплайна 1-ого порядка, согдасно
/// формулам http://www.math.spbu.ru/ru/mmeh/AspDok/pub/2010/Chashnikov.pdf (страница 6).
/// </summary>
/// <param name="n">Число узлов между полюсами.</param>
/// <param name="m">Число полюсов.</param>
/// <returns>Коэфициенты дискретного периодического Q-сплайна 1-ого порядка.</returns>
private static double[] CalculateQSpline(int n, int m) {
var N = n * m;
var qSpline = new double[N];
for (var j = 0; j < N; ++j) {
if (j >= 0 && j <= n - 1) {
qSpline[j] = (1.0 * n - j) / n;
}
if (j >= n && j <= N - n) {
qSpline[j] = 0;
}
if (j >= N - n + 1 && j <= N - 1) {
qSpline[j] = (1.0 * j - N + n) / n;
}
}
return qSpline;
}
Деление на n в коде объясняется тем, что нам необходимо вычислять нормализованные коэффициенты, согласно формуле:
Вычисленные коэффициенты нужны для построения результирующего сплайна. Для этого воспользуемся формулой:
/// <summary>
/// Вычисляет вектора дискретного периодического сплайна с векторными коэфициентами, согласно
/// формулам http://www.math.spbu.ru/ru/mmeh/AspDok/pub/2010/Chashnikov.pdf (страница 7).
/// </summary>
/// <param name="vectors">Полюса сплайна.</param>
/// <param name="qSpline">Коэффициенты B-сплайна 1-ого порядка.</param>
/// <param name="r">Порядок сплайна.</param>
/// <param name="n">Число узлов между полюсами сплайна.</param>
/// <param name="m">Число полюсов.</param>
/// <returns></returns>
private static Vector2[] CalculateSSpline(Vector2[] aVectors, double[] aQSpline, int r, int n, int m) {
var N = n * m;
var sSpline = new Vector2[r + 1][];
for (var i = 1; i <= r; ++i) {
sSpline[i] = new Vector2[N];
}
for (var j = 0; j < N; ++j) {
sSpline[1][j] = new Vector2(0, 0);
for (var p = 0; p < m; ++p) {
sSpline[1][j] += aVectors[p] * aQSpline[GetPositiveIndex(j - p * n, N)];
}
}
for (var v = 2; v <= r; ++v) {
for (var j = 0; j < N; ++j) {
sSpline[v][j] = new Vector2(0, 0);
for (var k = 0; k < N; ++k) {
sSpline[v][j] += aQSpline[k] * sSpline[v - 1][GetPositiveIndex(j - k, N)];
}
sSpline[v][j] /= n;
}
}
return sSpline[r];
}
Здесь, для удобства вычислений, был реализован класс Vector2 с минимально необходимым набором методов и операций.
/// <summary>
/// Реализует методы и операции для работы с 2-д вектором.
/// </summary>
public class Vector2 {
public double X;
public double Y;
/// <summary>
/// Конструктор по-умолчанию.
/// </summary>
public Vector2() {
this.X = 0;
this.Y = 0;
}
/// <summary>
/// Конструктор. Принимает координаты.
/// </summary>
/// <param name="x">Координата Х.</param>
/// <param name="y">Координата Y.</param>
public Vector2(double x, double y) {
this.X = x;
this.Y = y;
}
/// <summary>
/// Конструктор. Принимает Другое вектор.
/// </summary>
/// <param name="v">Исходный вектор.</param>
public Vector2(Vector2 v) {
X = v.X;
Y = v.Y;
}
/// <summary>
/// Реализует сложение векторов.
/// </summary>
/// <param name="a">Первый вектор.</param>
/// <param name="b">Второй вектор.</param>
/// <returns>Результат сложения.</returns>
public static Vector2 operator +(Vector2 a, Vector2 b) {
return new Vector2(a.X + b.X, a.Y + b.Y);
}
/// <summary>
/// Реализует вычитание векторов.
/// </summary>
/// <param name="a">Первый вектор.</param>
/// <param name="b">Второй вектор.</param>
/// <returns>Результат вычитания.</returns>
public static Vector2 operator -(Vector2 a, Vector2 b) {
return new Vector2(a.X - b.X, a.Y - b.Y);
}
/// <summary>
/// Реализует унарный минус.
/// </summary>
/// <param name="a">Исходный вектор.</param>
/// <returns>Результат применения унарного минуса.</returns>
public static Vector2 operator -(Vector2 a) {
return new Vector2(-a.X, -a.Y);
}
/// <summary>
/// Реализует умножение вектора на число.
/// </summary>
/// <param name="a">Исходный вектор.</param>
/// <param name="d">Число.</param>
/// <returns>Рузельтат умножения вектора на число.</returns>
public static Vector2 operator *(Vector2 a, double d) {
return new Vector2(a.X * d, a.Y * d);
}
/// <summary>
/// Реализует умножение числа на вектор.
/// </summary>
/// <param name="d">Число.</param>
/// <param name="a">Исходный вектор.</param>
/// <returns>Результат умножения.</returns>
public static Vector2 operator *(double d, Vector2 a) {
return a * d;
}
/// <summary>
/// Реализует умножение вектора на число.
/// </summary>
/// <param name="a">Исходный вектор.</param>
/// <param name="f">Число.</param>
/// <returns>Рузельтат умножения вектора на число.</returns>
public static Vector2 operator *(Vector2 a, float f) {
return a * (double)f;
}
/// <summary>
/// Реализует умножение числа на вектор.
/// </summary>
/// <param name="f">Число.</param>
/// <param name="a">Исходный вектор.</param>
/// <returns>Результат умножения.</returns>
public static Vector2 operator *(float f, Vector2 a) {
return a * (double)f;
}
/// <summary>
/// Реализует деление вектора на число.
/// </summary>
/// <param name="a">Исходный вектор.</param>
/// <param name="d">Число.</param>
/// <returns>Результат деления вектора на число.</returns>
public static Vector2 operator /(Vector2 a, double d) {
return new Vector2(a.X / d, a.Y / d);
}
/// <summary>
/// Реализует деление вектора на число.
/// </summary>
/// <param name="a">Исходный вектор.</param>
/// <param name="f">Число.</param>
/// <returns>Результат деления вектора на число.</returns>
public static Vector2 operator /(Vector2 a, float f) {
return a / (double)f;
}
}
А для обеспечения периодичности сплайна реализован небольшой метод GetPositiveIndex. Прошу сильно не ругаться из-за этого, просто не хотелось городить огород из-за мелочи.
/// <summary>
/// Обеспечивает периодичность для заданного множества.
/// </summary>
/// <param name="j">Индекс элемента.</param>
/// <param name="N">Количество элементов множества.</param>
/// <returns>Периодический индекс элемента.</returns>
private static int GetPositiveIndex(int j, int N) {
if (j >= 0) {
return j % N;
}
return N - 1 + ((j + 1) % N);
}
Ну, собственно, это вся реализация, которая дает нам вот такую картинку:
Стоп! Что это?! Сплайн не проходит через полюса! Незадача.
Для того, чтобы решить эту проблему и сделать интерполяционный сплайн, необходимо сделать предподготовку, а именно, воспользоваться следующей формулой для пересчета полюсов, но уже из другой работы.
Здесь apv — новые полюса, над которыми будем проводить все расчеты, изложенные выше, а последняя часть — это дискретное преобразование Фурье m-периодического сигнала.
/// <summary>
/// Пересчитывает коэффициенты сплайна для того, чтобы результирующий сплайн проходил через полюса.
/// http://dha.spb.ru/PDF/discreteSplines.pdf (страница 6 и 7).
/// </summary>
/// <param name="aPoints">Исходные точки.</param>
/// <param name="r">Порядок сплайна.</param>
/// <param name="n">Количество узлов между полюсами сплайна.</param>
/// <param name="m">Количество полюсов.</param>
/// <returns>Новые вектора.</returns>
private static Vector2[] RecalculateVectors(Vector2[] aPoints, int r, int n, int m) {
var N = n * m;
// Вычисляем знаменатель.
var tr = new double[m];
tr[0] = 1;
for (var k = 1; k < m; ++k) {
for (var q = 0; q < n; ++q) {
tr[k] += Math.Pow(2 * n * Math.Sin((Math.PI * (q * m + k)) / N), -2 * r);
}
tr[k] *= Math.Pow(2 * Math.Sin((Math.PI * k) / m), 2 * r);
}
// Вычисляем числитель.
var zre = new Vector2[m];
var zim = new Vector2[m];
for (var j = 0; j < m; ++j) {
zre[j] = new Vector2(0, 0);
zim[j] = new Vector2(0, 0);
for (var k = 0; k < m; ++k) {
zre[j] += aPoints[k] * Math.Cos((-2 * Math.PI * j * k) / m);
zim[j] += aPoints[k] * Math.Sin((-2 * Math.PI * j * k) / m);
}
}
// Считаем результат.
var result = new Vector2[m];
for (var p = 0; p < m; ++p) {
result[p] = new Vector2(0, 0);
for (var k = 0; k < m; ++k) {
var d = (zre[k] * Math.Cos((2 * Math.PI * k * p) / m)) - (zim[k] * Math.Sin((2 * Math.PI * k * p) / m));
d *= 1.0 / tr[k];
result[p] += d;
}
result[p] /= m;
}
return result;
}
/// <summary>
/// Вычисляет узловые точки дискретного N-периодического сплайна с векторными коэффициентами.
/// </summary>
/// <param name="aPoints">Полюса сплайна (исходные точки). Должно быть не менее 2-х полюсов.</param>
/// <param name="r">Порядок сплайна.</param>
/// <param name="n">Число узлов между полюсами сплайна.</param>
/// <param name="aIsIncludeOriginalPoints">True - сплайн будет проходить через полюса, false - сплайн не будет проходить через полюса.</param>
/// <returns></returns>
public static Vector2[] Calculate(Vector2[] aPoints, int r, int n = 5, bool aIsIncludeOriginalPoints = true) {
if (aPoints == null) {
throw new ArgumentNullException("aPoints");
}
if (aPoints.Length <= 2) {
throw new ArgumentException("Число полюсов должно быть > 2.");
}
if (r <= 0) {
throw new ArgumentException("Порядок сплайна должен быть > 0.");
}
if (n < 1) {
throw new ArgumentException("Число узлов между полюсами сплайна должно быть >= 1.");
}
var m = aPoints.Length;
var N = n * m;
Vector2[] vectors;
if (aIsIncludeOriginalPoints) {
vectors = RecalculateVectors(aPoints, r, n, m);
} else {
vectors = new Vector2[m];
aPoints.CopyTo(vectors, 0);
}
var qSpline = CalculateQSpline(n, m);
var resultPoints = CalculateSSpline(vectors, qSpline, r, n, m);
return resultPoints;
}
Применяя формулы выше добиваемся решения поставленной задачи:
Стоит уточнить, что на вид сплайна непосредственное влияние оказывает порядок (нумерация) полюсов. Эксперименты с n и r оставляю на откуп читателю.
По-моему то, что надо!
P.s Еще раз огромное спасибо Николаю Чашникову за его работу и за помощь!
P.p.s Все исходники лежат здесь.
Комментарии (18)
Boomburum
06.09.2016 09:36+3Смотрел на картинку до ката и только потом понял, что на ней написано Habr :)
prostofilya
06.09.2016 09:45а готовых реализаций алгоритма я не обнаружил, впрочем, как и для других языков и технологий.
Python + scipy
Правда, за реализацией придётся лезть внутрь, либо в доки scipy.
hombit
06.09.2016 10:28+1В этом примере ничего сложнее классических кубических сплайнов не используется, в документации scipy даны ссылки на классические книги
appsforlife
06.09.2016 10:48+2А catmull rom сплайны не пробовал использовать? По одному на сегмент кривой.
Denxc
06.09.2016 11:18+1Не пробовал. После вашего комментария за 10 минут накидал реализацию, так как много информации по ним. Удивлен.
Спасибо!
Denxc
06.09.2016 13:28Поигрался с этим алгоритмом. Вполне рабочий, легкий в реализации. Но, к моей практической задачи не подошел из-за некоторых особенностей построения сплайна. Здесь исходники.
vertu77
06.09.2016 11:46Спасибо.
Может кто подскажет ответ на вопрос по близкой тематике?
Есть набор точек (отсчеты АЦП), расположенных на кривой, несколько напоминающей вершинку синусоиды (но крутизна переднего и заднего фронтов разная). Хотелось бы интерполировать эту кривую и найти ее экстремум.
В сторону каких алгоритмов смотреть с учетом ограниченных вычислительных ресурсов (предполагается использование в микроконтроллере)?
fireSparrow
06.09.2016 13:29А вы не подскажете, что вообще стоит почитать с точки зрения отрисовки кривых?
Интересуют не сколько конкретные алгоритмы, сколько база.
Хочется разобраться, чтобы понимать, что вообще происходит в этих алгоритмах и потом дорабатывать и писать новые под свои нужды.
Сейчас, например, передо мной стоит задача отрисовки кольца (то есть двух концентрических окружностей с заполнением пространства между ними). Сложность в том, что координаты центра не целочисленные, ширина кольца произвольная и тоже не целочисленная, и контуры должны быть сглаженными. Ничего из нагугливающихся стандартных алгоритмов не могу прикрутить к этой задаче.
В какую сторону читать и по каким словам гуглить?
И ещё не могу понять, как отрисовывать кривые Безье произвольной ширины со сглаживанием контуров.
Aleksanderis
06.09.2016 15:10Как раз на днях искал, что-то подобное про «рисование плавных траекторий». Увидел нечаянно это и осенило — похоже это что-то близкое к тому, что мне надо. Сам, хоть и будучи программистом, никогда не был силён в «высшей математике» (во всяких этих кривых, полях и прочих интеграллах), так что немного сомневаюсь… Вроде как и нужный результат у вашего алгоритма, но вроде как параметры другие.
Как вам кажется — годится ли подобный алгоритм для нахождения _плавной_ траектории между заданными координатами, и что бы траектория проходила по этим точкам? Или тут всё же другая история?
Смотрел про кривые Безье (как более часто встречаемые), так они неособо подошли — там суть в самих кривых, а не в траектории, проходящей через точки.
IIvana
07.09.2016 03:36Я бы по рабоче-крестьянски написал так: вводим систему координат (она по умолчанию у нас и так присутствует), параметризуем точки по каждой координате — то есть получаем массив иксов и массив игреков, к каждому из этих массивов добавляем массив параметра — можно равномерный шаг, можно пропорциональный расстоянию между точками на плоскости — и интерполируем эти 2 набора точек от параметра ЛЮБЫМ понравившимся нам способом с нужной степенью гладкости — хоть глобальным сплайном с непрерывными 0,1,2 производными, хоть локальным с непрерывными 0,1 — тот же Катмулл-Ром и прочие, имя им легион. Все, пробегаемся с любым шагом по параметру и получаем точки нашей гладкой кривой.
Leo5700
13.09.2016 18:13+1Спасибо за статью. Тема и правда в сети освещена мало, инструментарий пригодится мне для интерполяции номограмм.
Mendel
Спасибо, а не подскажете примеры где это применимо на практике?
Сразу захотелось подумать в сторону того как сделать перестановку точек чтобы выбирать минимальную длину, но следующей мыслью было — если применимо для реальной задачи. А реальную задачу как-то и не придумать.
Denxc
Автоматическое проектирование, геометрическое моделирование и тд.
Конкретно я буду использовать этот алгоритм для контроля геометрических характеристик цилиндрических алюминиевых слитков на одном из российский заводов. Есть лазерные датчики, которые дают множество точек поверхности, затем по этим точкам генерируется интерполяционная поверхность.