Всем привет! Недавно возникла практическая необходимость использовать интерполяцию для замкнутых кривых. Проект разрабатывался под .Net на C#, а готовых реализаций алгоритма я не обнаружил, впрочем, как и для других языков и технологий. В результате пришлось самому изучить мат.часть существующих методов и написать свою библиотеку. Наработками и готовым решением готов поделиться с вами.
Поиск в гугле на тему «интерполяция замкнутых кривых» не дает особых результатов, кроме как на автореферат из диссертации от 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 в коде объясняется тем, что нам необходимо вычислять нормализованные коэффициенты, согласно формуле:
Вычисленные коэффициенты нужны для построения результирующего сплайна. Для этого воспользуемся формулой:
, здесь v — порядок сплайна, ap — полюса сплайна.
/// <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 Все исходники лежат здесь.
Автор: Denxc
