Интерполяция замкнутых кривых

в 5:49, , рубрики: .net, C#, Алгоритмы, алгоритмы обработки данных, интерполяция, математика, метки: ,

Всем привет! Недавно возникла практическая необходимость использовать интерполяцию для замкнутых кривых. Проект разрабатывался под .Net на C#, а готовых реализаций алгоритма я не обнаружил, впрочем, как и для других языков и технологий. В результате пришлось самому изучить мат.часть существующих методов и написать свою библиотеку. Наработками и готовым решением готов поделиться с вами.

Интерполяция замкнутых кривых - 1


Поиск в гугле на тему «интерполяция замкнутых кривых» не дает особых результатов, кроме как на автореферат из диссертации от 2010 года господина Чашникова Н.В. Именно эту работу я взял за теоретическую основу. Должен сказать, что в процессе анализа математики возникали трудности, которые помог решить сам Николай Викторович, за что ему огромное спасибо! Собственно, дальше будет разбор этого автореферата с конкретными кусками кода.

Уточнение от Николая

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

Итак, задача стоит так: на входе есть набор двумерных точек, необходимо построить интерполяционную замкнутую кривую.

Для решения этой задачи можно использовать дискретные N-периодические сплайны с векторными коэффициентами.
Далее я буду использовать следующую терминологию и обозначения:

Полюса сплайна — набор исходных точек (векторов), m — количество полюсов сплайна, n — количество узлов сплайна между двумя соседними полюсами, N = n * m — период сплайна, r — порядок сплайна.

B-сплайн первого порядка вводится вполне явно:

Интерполяция замкнутых кривых - 2

Ничего сложного

        /// <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 в коде объясняется тем, что нам необходимо вычислять нормализованные коэффициенты, согласно формуле:

Интерполяция замкнутых кривых - 3

Вычисленные коэффициенты нужны для построения результирующего сплайна. Для этого воспользуемся формулой:

Интерполяция замкнутых кривых - 4

, здесь 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 с минимально необходимым набором методов и операций.

Вот он - 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);
        }

Ну, собственно, это вся реализация, которая дает нам вот такую картинку:

Интерполяция замкнутых кривых - 5

Стоп! Что это?! Сплайн не проходит через полюса! Незадача.

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

Интерполяция замкнутых кривых - 6

Интерполяция замкнутых кривых - 7

Интерполяция замкнутых кривых - 8

Здесь 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;
        }

Применяя формулы выше добиваемся решения поставленной задачи:

Интерполяция замкнутых кривых - 9

Интерполяция замкнутых кривых - 10

Стоит уточнить, что на вид сплайна непосредственное влияние оказывает порядок (нумерация) полюсов. Эксперименты с n и r оставляю на откуп читателю.

По-моему то, что надо!

P.s Еще раз огромное спасибо Николаю Чашникову за его работу и за помощь!

P.p.s Все исходники лежат здесь.

Автор: Denxc

Источник

Поделиться новостью

* - обязательные к заполнению поля