Всем привет! Недавно возникла практическая необходимость использовать интерполяцию для замкнутых кривых. Проект разрабатывался под .Net на C#, а готовых реализаций алгоритма я не обнаружил, впрочем, как и для других языков и технологий. В результате пришлось самому изучить мат.часть существующих методов и написать свою библиотеку. Наработками и готовым решением готов поделиться с вами.
Поиск в гугле на тему «интерполяция замкнутых кривых» не дает особых результатов, кроме как на автореферат из диссертации от 2010 года господина Чашникова Н.В. Именно эту работу я взял за теоретическую основу. Должен сказать, что в процессе анализа математики возникали трудности, которые помог решить сам Николай Викторович, за что ему огромное спасибо! Собственно, дальше будет разбор этого автореферата с конкретными кусками кода.
Итак, задача стоит так: на входе есть набор двумерных точек, необходимо построить интерполяционную замкнутую кривую.
Для решения этой задачи можно использовать дискретные N-периодические сплайны с векторными коэффициентами.
Далее я буду использовать следующую терминологию и обозначения:
Полюса сплайна — набор исходных точек (векторов), m — количество полюсов сплайна, n — количество узлов сплайна между двумя соседними полюсами, N = n * m — период сплайна, r — порядок сплайна.
B-сплайн первого порядка вводится вполне явно:
Деление на n в коде объясняется тем, что нам необходимо вычислять нормализованные коэффициенты, согласно формуле:
Вычисленные коэффициенты нужны для построения результирующего сплайна. Для этого воспользуемся формулой:
, здесь v — порядок сплайна, ap — полюса сплайна.
Здесь, для удобства вычислений, был реализован класс Vector2 с минимально необходимым набором методов и операций.
А для обеспечения периодичности сплайна реализован небольшой метод GetPositiveIndex. Прошу сильно не ругаться из-за этого, просто не хотелось городить огород из-за мелочи.
Ну, собственно, это вся реализация, которая дает нам вот такую картинку:
Стоп! Что это?! Сплайн не проходит через полюса! Незадача.
Для того, чтобы решить эту проблему и сделать интерполяционный сплайн, необходимо сделать предподготовку, а именно, воспользоваться следующей формулой для пересчета полюсов, но уже из другой работы.
Здесь apv — новые полюса, над которыми будем проводить все расчеты, изложенные выше, а последняя часть — это дискретное преобразование Фурье m-периодического сигнала.
Применяя формулы выше добиваемся решения поставленной задачи:
Стоит уточнить, что на вид сплайна непосредственное влияние оказывает порядок (нумерация) полюсов. Эксперименты с n и r оставляю на откуп читателю.
По-моему то, что надо!
P.s Еще раз огромное спасибо Николаю Чашникову за его работу и за помощь!
P.p.s Все исходники лежат здесь.
Поиск в гугле на тему «интерполяция замкнутых кривых» не дает особых результатов, кроме как на автореферат из диссертации от 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 с минимально необходимым набором методов и операций.
Вот он - 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. Прошу сильно не ругаться из-за этого, просто не хотелось городить огород из-за мелочи.
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 Все исходники лежат здесь.
