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

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

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

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

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

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

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

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

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

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

        /// 
        /// Вычисляет коэфициенты дискретного периодического Q-сплайна 1-ого порядка, согдасно 
        /// формулам http://www.math.spbu.ru/ru/mmeh/AspDok/pub/2010/Chashnikov.pdf (страница 6).
        /// 
        /// Число узлов между полюсами.
        /// Число полюсов.
        /// Коэфициенты дискретного периодического Q-сплайна 1-ого порядка.
        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 — полюса сплайна.

Реализуем вычисление

        /// 
        /// Вычисляет вектора дискретного периодического сплайна с векторными коэфициентами, согласно
        /// формулам http://www.math.spbu.ru/ru/mmeh/AspDok/pub/2010/Chashnikov.pdf (страница 7).
        /// 
        /// Полюса сплайна.
        /// Коэффициенты B-сплайна 1-ого порядка.
        /// Порядок сплайна.
        /// Число узлов между полюсами сплайна.
        /// Число полюсов.
        /// 
        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

    /// 
    /// Реализует методы и операции для работы с 2-д вектором.
    /// 
    public class Vector2 {

        public double X;
        public double Y;

        /// 
        /// Конструктор по-умолчанию.
        /// 
        public Vector2() {
            this.X = 0;
            this.Y = 0;
        }

        /// 
        /// Конструктор. Принимает координаты.
        /// 
        /// Координата Х.
        /// Координата Y.
        public Vector2(double x, double y) {
            this.X = x;
            this.Y = y;
        }        

        /// 
        ///  Конструктор. Принимает Другое вектор.
        /// 
        /// Исходный вектор.
        public Vector2(Vector2 v) {
            X = v.X;
            Y = v.Y;
        }

        /// 
        /// Реализует сложение векторов.
        /// 
        /// Первый вектор.
        /// Второй вектор.
        /// Результат сложения.
        public static Vector2 operator +(Vector2 a, Vector2 b) {
            return new Vector2(a.X + b.X, a.Y + b.Y);
        }

        /// 
        /// Реализует вычитание векторов.
        /// 
        /// Первый вектор.
        /// Второй вектор.
        /// Результат вычитания.
        public static Vector2 operator -(Vector2 a, Vector2 b) {
            return new Vector2(a.X - b.X, a.Y - b.Y);
        }

        /// 
        /// Реализует унарный минус.
        /// 
        /// Исходный вектор.
        /// Результат применения унарного минуса.
        public static Vector2 operator -(Vector2 a) {
            return new Vector2(-a.X, -a.Y);
        }

        /// 
        /// Реализует умножение вектора на число.
        /// 
        /// Исходный вектор.
        /// Число.
        /// Рузельтат умножения вектора на число.
        public static Vector2 operator *(Vector2 a, double d) {
            return new Vector2(a.X * d, a.Y * d);
        }

        /// 
        /// Реализует умножение числа на вектор.
        /// 
        /// Число.
        /// Исходный вектор.
        /// Результат умножения.
        public static Vector2 operator *(double d, Vector2 a) {
            return a * d;
        }

        /// 
        /// Реализует умножение вектора на число.
        /// 
        /// Исходный вектор.
        /// Число.
        /// Рузельтат умножения вектора на число.
        public static Vector2 operator *(Vector2 a, float f) {
            return a * (double)f;
        }

        /// 
        /// Реализует умножение числа на вектор.
        /// 
        /// Число.
        /// Исходный вектор.
        /// Результат умножения.
        public static Vector2 operator *(float f, Vector2 a) {
            return a * (double)f;
        }

        /// 
        /// Реализует деление вектора на число.
        /// 
        /// Исходный вектор.
        /// Число.
        /// Результат деления вектора на число.
        public static Vector2 operator /(Vector2 a, double d) {
            return new Vector2(a.X / d, a.Y / d);
        }

        /// 
        /// Реализует деление вектора на число.
        /// 
        /// Исходный вектор.
        /// Число.
        /// Результат деления вектора на число.
        public static Vector2 operator /(Vector2 a, float f) {
            return a / (double)f;
        }
    }

А для обеспечения периодичности сплайна реализован небольшой метод GetPositiveIndex. Прошу сильно не ругаться из-за этого, просто не хотелось городить огород из-за мелочи.

GetPositiveIndex

        /// 
        /// Обеспечивает периодичность для заданного множества.
        /// 
        /// Индекс элемента.
        /// Количество элементов множества.
        /// Периодический индекс элемента.
        private static int GetPositiveIndex(int j, int N) {
            if (j >= 0) {
                return j % N;
            }

            return N - 1 + ((j + 1) % N);
        }

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

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

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

Здесь apv — новые полюса, над которыми будем проводить все расчеты, изложенные выше, а последняя часть — это дискретное преобразование Фурье m-периодического сигнала.

Пересчет полюсов сплайна

        /// 
        /// Пересчитывает коэффициенты сплайна для того, чтобы результирующий сплайн проходил через полюса.
        /// http://dha.spb.ru/PDF/discreteSplines.pdf (страница 6 и 7). 
        /// 
        /// Исходные точки.
        /// Порядок сплайна.
        /// Количество узлов между полюсами сплайна.
        /// Количество полюсов.
        /// Новые вектора.
        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;
        }

Итоговый метод вычисления сплайна

        /// 
        /// Вычисляет узловые точки дискретного N-периодического сплайна с векторными коэффициентами.
        /// 
        /// Полюса сплайна (исходные точки). Должно быть не менее 2-х полюсов.
        /// Порядок сплайна.
        /// Число узлов между полюсами сплайна.
        /// True - сплайн будет проходить через полюса, false - сплайн не будет проходить через полюса.
        /// 
        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 Все исходники лежат здесь.

Комментарии (4)

Mendel

6 сентября 2016 в 09:13

0

↑

↓

Спасибо, а не подскажете примеры где это применимо на практике?
Сразу захотелось подумать в сторону того как сделать перестановку точек чтобы выбирать минимальную длину, но следующей мыслью было — если применимо для реальной задачи. А реальную задачу как-то и не придумать.
- Denxc
  
  6 сентября 2016 в 09:20
  
  0
  
  ↑
  
  ↓
  
  Автоматическое проектирование, геометрическое моделирование и тд.
  Конкретно я буду использовать этот алгоритм для контроля геометрических характеристик цилиндрических алюминиевых слитков на одном из российский заводов. Есть лазерные датчики, которые дают множество точек поверхности, затем по этим точкам генерируется интерполяционная поверхность.
Boomburum

6 сентября 2016 в 09:36

0

↑

↓

Смотрел на картинку до ката и только потом понял, что на ней написано Habr:)
prostofilya

6 сентября 2016 в 09:45 (комментарий был изменён)

0

↑

↓

а готовых реализаций алгоритма я не обнаружил, впрочем, как и для других языков и технологий.

Python + scipy
Правда, за реализацией придётся лезть внутрь, либо в доки scipy.