RU: Curves. Interpolation and Approximation

Геометрические примитивы

Кривые: интерполяция и аппроксимация

В этом разделе мы кратко напомним постановку задач интерполяции и аппроксимации набора точек гладкой кривой.

Интерполяционный многочлен

Пусть в пространстве заданы $n+1$ точек $P_0, P_1, P_2,\dots, P_n$, которые имеют радиус-векторы $\mathbf{r}_0, \mathbf{r}_1, \mathbf{r}_2,\dots, \mathbf{r}_n$. Задача интерполяции состоит в построении кривой, проходящей через указанные точки в указанном порядке.

Конечно, через фиксированный набор точек можно провести бесконечное число кривых, так что задача интерполяции не имеет единственного решения. Далее мы построим несколько типов кривых и обсудим их достоинства и недостатки. Будем строить кривые в виде $\mathbf{r}(t)$, где параметр $t$ изменяется на некотором отрезке $[a, b]$: $a\leq t \leq b$. Введем на отрезке $[a, b]$ сетку $\{t_i\}$ из $n+1$ точек: $a=t_0 < t_1 < t_2<\dots < t_n=b$ и потребуем, чтобы при значении параметра $t=t_i$ кривая проходила через точку $P_i$, так что $\mathbf{r}(t_i)=\mathbf{r}_i$. Введение параметризации и сетки может быть выполнено различными способами. Обычно выбирают либо равномерную сетку, полагая $a=0$, $b=n$, $t_i=i$, либо, что более предпочтительно, соединяют точки отрезками и в качестве разности значений параметра $t_{i+1}-t_i$ берут длину отрезка $\mathbf{r}_{i+1}-\mathbf{r}_{i}$.

Одним из распространенных способов интерполяции является использование кривой в виде многочлена от $t$ степени $n$, т.е. в виде функции

(1)
\begin{align} \mathbf{r}(t) = \mathbf{p}^{(n)}(t) =\sum_{k=0}^n \mathbf{a}_k t^k \end{align}

Многочлен имеет $n+1$ коэффициентов $\mathbf{a}_k$, которые можно найти из условий $\mathbf{r}(t_i)=\mathbf{r}_i$. Эти условия приводят к системе линейных уравнений для коэффициентов $\mathbf{a}_k$:

(2)
\begin{align} \sum_{k=0}^n \mathbf{a}_k t_i^k = \mathbf{r}_i, \qquad 0\leq i \leq n \end{align}

или в матричном виде

(3)
\begin{align} \begin{pmatrix} 1 & t_0 & t_0^2 & \ldots & t_0^n \\ 1 & t_1 & t_1^2 & \ldots & t_1^n \\ 1 & t_2 & t_2^2 & \ldots & t_2^n \\ \vdots & \vdots& \vdots& & \vdots \\ 1 & t_n & t_n^2 & \ldots & t_n^n \end{pmatrix} \begin{pmatrix} \mathbf{a}_0 \\ \mathbf{a}_1 \\ \mathbf{a}_2 \\ \vdots \\ \mathbf{a}_n \end{pmatrix}= \begin{pmatrix} \mathbf{r}_0 \\ \mathbf{r}_1 \\ \mathbf{r}_2 \\ \vdots \\ \mathbf{r}_n \end{pmatrix} \end{align}

Подчеркнем, что нужно решить три системы уравнений: для $x$, $y$ и $z$ координат. Все они имеют одну матрицу коэффициентов, обращая которую, по значениям радиус-векторов $\mathbf{r}_i$ точек вычисляются векторы $\mathbf{a}_k$ коэффициентов многочлена. Определитель матрицы

(4)
\begin{align} W(t_0, t_1, \ldots , t_n) = \begin{vmatrix} 1 & t_0 & t_0^2 & \ldots & t_0^n \\ 1 & t_1 & t_1^2 & \ldots & t_1^n \\ 1 & t_2 & t_2^2 & \ldots & t_2^n \\ \vdots & \vdots& \vdots& & \vdots \\ 1 & t_n & t_n^2 & \ldots & t_n^n \end{vmatrix} = \prod_{\substack{j,j\\i>j}} (t_i -t_j) \end{align}

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

Обсудим достоинства и недостатки интерполяционного многочлена. К достоинствам можно отнести его следующие свойства:

  • Для заданного набора точек и сетки параметра кривая строится однозначно.
  • Кривая имеет непрерывные производные любого порядка.

Основные недостатки состоят в том, что:

  • С ростом числа точек порядок многочлена возрастает, а вместе с ним возрастает число операций, которое нужно выполнить для вычисления точки на кривой.
  • С ростом числа точек у интерполяционной кривой могут возникнуть осцилляции (см. пример ниже).

Классическим примером, показывающим возникновение осцилляций у интерполяционного многочлена, служит нтерполяция на равномерной сетке значений функции

(5)
\begin{align} f(x) = \frac{1}{1+x^2} \end{align}

Введем на отрезке $[-5, 5]$ равномерную сетку $x_i=-5+ih$, $h=10/n$, $0 \leq i \leq n$ и рассмотрим поведение многочлена $y(x) = \sum_{i=0}^n a_i x^i$ который в точках $x_i$ принимает значения $y_i=1/(1+x_i^2)$. На рисунке справа представлено поведение самой функции (штрих-пунктирная линия) и трех интерполяционных кривых при $n=2, 4, 8$:

  • интерполяция на сетке $x_0=-5$, $x_1=0$, $x_2=5$ — квадратичная парабола;
  • интерполяция на сетке $x_0=-5$, $x_1=-2.5$, $x_2=0$, $x_3=2.5$, $x_4=5$ — многочлен четвертой степени;
  • интерполяция на сетке $x_0=-5$, $x_1=-3.75$, $x_2=-2.5$, $x_3=-1.25$, $x_4=0$, $x_5=1.25$, $x_6=2.5$, $x_7=-3.75$, $x_8=5$ — многочлен восьмой степени.
Runge.jpg

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

Интерполяционный многочлен в форме Лагранжа

Коэффициенты интерполяционного многочлена линейно зависят от положения точек, так что, группируя члены, можно записать

(6)
\begin{align} \mathbf{p}^{(n)}(t) = \sum_{i=0}^n L_i^{(n)}(t)\mathbf{r}_i \end{align}

где $L_i^{(n)}(t)$ — некоторые многочлены степени $n$. Очевидно, они обладают следующими свойством: $L_i^{(n)}(t_j) = \delta_{ij}$, т.е. многочлен $L_i^{(n)}(t)$ равен нулю во всех узлах сетки, кроме узла $t_i$, в котором он имеет значение единица. Многочлены $L_i^{(n)}(t)$ называются многочленами Лагранжа.

В OpenCascade вычисление значения интерполяционного многочлена выполняется процедурой

    Standard_Integer  PLib::EvalLagrange(const Standard_Real Parameter,
           const Standard_Integer                DerivativeRequest,
           const Standard_Integer                Degree,
           const Standard_Integer                Dimension, 
           Standard_Real&                        Values,
           Standard_Real&                        Parameters,
           Standard_Real&                        Results)

Поясним назначение параметров. Процедура возвращает код ответа: 0 в случае успешного выполнения.
  • Parameter — входная величина, значение параметра, при котором нужно вычислить интерполяционную функцию;
  • DerivativeRequest — входная величина, число производных, которое нужно вычислить. Процедура может вычислить не только функцию, но и значения производных. Если производные вычислять не нужно, значение данного параметра — 0;
  • Degree — входная величина, степень интерполяционного многочлена $n$. Она должна быть на единицу меньше, чем число точек (см. далее);
  • Dimension — входная величина, размерность (число компонент у) точек. Если интрполируется функция $f(x)$, то Dimension=1, Если интерполируются точки на плоскости $\mathbf{r}_i=(x_i, y_i)$, то Dimension=2, если интерполируются точки в пространстве $\mathbf{r}_i=(x_i, y_i, z_i)$, то Dimension=3, и т.д.;
  • Values — входная величина, значения координат точек. Значения передаются в виде одномерного массива, где координаты хранятся последовательно. Для функции это последовательность значений $f_0, f_1, \ldots, f_n$, для точек на плоскости это последовательность $x_0, y_0, x_1, y_1 \ldots, x_n, y_n$, всего $2(n+1)$ значений, для точек в пространстве это $x_0, y_0, z_0, x_1, y_1, z_1 \ldots, x_n, y_n, z_n$, всего $3(n+1)$ значений и так далее. Обратим внимание на то, что если точки хранятся в программе в виде одномерного массива, то данный параметр есть ссылка на массив его значений:
    TColgp_Array1OfPnt PointsArray;
    Standard_Real *point_array;
    . . . 
    point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ; 
    . . . 
    PLib::EvalLagrange(<>, <>, <>, <>, point_array[0], ...) ;
  • Parameters — входная величина, значения параметров $t_i$ точек;
  • Results — выходная величина, значение интерполяционной функции и производных. Они хранятся последовательно, аналогично описанному выше.

В качестве примера вычислим значение интерполяционного многочлена для функции предыдущего примера

// Функция f(x) = 1 / (1+x^2)
    int i;
    int degree = 4;    // Степень многочлена
    Standard_Real values[5];    // Значения функции
    Standard_Real parameters[5];    // Значения узлов сетки параметра 
    Standard_Real results;        // Результат интерполяции
    Standard_Real parameter;    // Значение параметра, при котором нужно вычислить функцию
    for( i=0; i<5; i++ ) {
        parameters[i] = -5. + i*2.5;
        values[i] = 1. / (1. + pow(parameters[i],2) );
    }
    parameter = 3.8;
    PLib::EvalLagrange(parameter, 0, degree, 1, values[0], parameters[0], results);
//     results    -0.36433103448276

Интерполяция Эрмита

Интерполяция Эрмита является обобщением интерполяции Лагранжа, учитывающим не только значения функции, но и ее производных. Напомним, что первая производная векторной функции связана с направлением касательной, а вторая — с кривизной.

Рассмотрим случай заданных первых производных. Итак, пусть для $n+1$ значений параметра $t$, т.е. для $t_0 < t_1 < \ldots <t_n$ заданы положения точек в пространстве $\mathbf{r}(t_i)=\mathbf{r}_i$ и значения первых производных по параметру $\mathbf{r}'(t_i) =\mathbf{r}'_i$. Требуется построить векторную функцию $\mathbf{r}(t)$, удовлетворяющую этим условиям. Запишем ее в виде

(7)
\begin{align} \mathbf{r}(t)=\sum_{i=0}^n\left(H_{i0}(t)\mathbf{r}(t_i)+H_{i1}(t)\mathbf{r}'(t_i)\right) \end{align}

где $H_{i0}(t)$ и $H_{i1}(t)$ многочлены cтепени $2n+1$, которые должны быть найдены из сформулированных выше условий.

Рассмотрим подробнее наиболее употребительный на практике частный случай интерполирования по двум точкам, в которых заданы значения функции и производной. Полагая единичным диапазон изменения параметра: $0 \leq t \leq 1$, а узлы сетки равными $t_0=0$, $t_1=1$, получим вид многочленов Эрмита

(8)
\begin{align} H_{00}(t)&=(2t+1)(t-1)^2=2t^3-3t^2+1 & H_{10}(t)&=(3-2t)t^2=-2t^3+3t^2 \nonumber\\ H_{01}(t)&=t(t-1)^2=t^3-2t^2+t & H_{11}(t)&=(t-1)t^2=t^3-t^2 \end{align}

Интерполяционная формула в этом случае примет вид:

(9)
\begin{align} \mathbf{r}(t)=H_{00}(t)\mathbf{r}(0)+H_{10}(t)\mathbf{r}(1)+ H_{01}(t)\mathbf{r}'(0)+H_{11}(t)\mathbf{r}'(1) \end{align}

В OpenCascade вычисление значения интерполяционного многочлена Эрмита третьего порядка выполняется процедурой

    Standard_Integer PLib::EvalCubicHermite( const Standard_Real  Parameter,
        const Standard_Integer       DerivativeRequest,
        const Standard_Integer       Dimension, 
        Standard_Real&               Values,
        Standard_Real&               Derivatives,
        Standard_Real&               theParameters,
        Standard_Real&               Results)

Поясним назначение параметров. Процедура возвращает код ответа: 0 в случае успешного выполнения.
  • Parameter — входная величина, значение параметра, при котором нужно вычислить интерполяционную функцию;
  • DerivativeRequest — входная величина, число производных, которое нужно вычислить. Процедура может вычислить не только функцию, но и значения производных. Если производные вычислять не нужно, значение данного параметра — 0;
  • Dimension — входная величина, размерность (число компонент у) точек. Если интрполируется функция $f(x)$, то Dimension=1, Если интерполируются точки на плоскости $\mathbf{r}_i=(x_i, y_i)$, то Dimension=2, если интерполируются точки в пространстве $\mathbf{r}_i=(x_i, y_i, z_i)$, то Dimension=3, и т.д.;
  • Values — входная величина, значения координат двух точек. Значения передаются в виде одномерного массива, где координаты хранятся последовательно. Для функции это последовательность значений $f_0, f_1$, для точек на плоскости это последовательность $x_0, y_0, x_1, y_1$, для точек в пространстве это $x_0, y_0, z_0, x_1, y_1, z_1$ и так далее.
  • Derivatives— входная величина, значения производных в двух точках. Они хранятся аналогично координатам.
  • Parameters — входная величина, значения параметров $t_0, t_1$ (два значения);
  • Results — выходная величина, значение интерполяционной функции и производных. Они хранятся последовательно, аналогично описанному выше.

Сплайн

Как мы видели выше, задача интерполяции состоит в нахождении $n+1$ функций $c_i(t)$ параметра $t$, таких, что векторная функция $\mathbf{r}(t) = \sum_{i=0}^n c_i(t)\mathbf{r}_i$ принимает при заданных значениях параметра $t_i$ заданные значения $\mathbf{r}_i$. Предполагается, что коэффициенты зависят от параметра и сетки, но не зависят от положения точек $\mathbf{r}_i$. Перечислим некоторые свойства коэффициентов $c_i(t)$.

  • поскольку коэффициенты не зависят от положения точек, должны выполняться условия $c_i(t_j) = \delta_{ij}$, т.е. функция $c_i(t)$ должна быть равна нулю во всех узлах $t_j$ сетки, кроме узла $t_i$, в котором она должна иметь значение единица. Если это свойство не выполнено, то не будет выполнено условие $\mathbf{r}(t_i)=\mathbf{r}_i$. Такая ситуация возможна для сглаживающих кривых, которые будут рассматриваться ниже;
  • естественно требовать, чтобы в случае набора совпадающих точек $\mathbf{r}_0= \mathbf{r}_1 =\mathbf{r}_2=\dots= \mathbf{r}_n$ интерполяционная кривая представляла собой эту точку. Отсюда следует условие $\sum_{i=0}^n c_i(t) = 1$;
  • в качестве коэффициентов $c_i(t)$ мы будем расматривать многочлены (быть может различные на разных интервалах), т.е., например, мы не будем здесь расматривать интерполяцию тригонометрическими функциями (суммы Фурье) или экспонентами.
  • желательно, чтобы коэффициенты $c_i(t)$ были непрерывны и имели непрерывные производные.

Наглядный смысл функции $c_i(t)$ состоит в том, что она определяет влияние положения $\mathbf{r}_i$ точки $i$ на значение $\mathbf{r}(t)$ интерполяционной функции. Поэтому желательно, чтобы коэффициенты $c_i(t)$ удовлетворяли условию $0 \leq c_i(t) \leq 1$. Если это свойство выполнено, то интерполяционная кривая лежит внутри выпуклой оболочки, натянутой на точки $\mathbf{r}_0, \mathbf{r}_1, \mathbf{r}_2,\dots, \mathbf{r}_n$ (см. обсуждение далее).

При построении сплайна, в отличие от интерполяционного многочлена, на каждом интервале $\Delta_i=[t_i,t_{i+1}]$ сетки параметра строится свой многочлен, но при этом накладываются условия, чтобы в точках сетки $t_i$ функция и какое-то число ее производных были непрерывны. Условий на гладкость выбирается столько, чтобы можно было однозначно
определить функции $c_i(t)$. Поскольку векторная функция $\mathbf{r}(t)$ эквивалентна трем скалярным функциям $(x(t), y(t), z(t))$, мы рассмотрим свойства сплайна для скалярной функции $f(t)$.

Пусть на отрезке $[a,b]$ задана сетка узлов $\{t_i\}$: $a=t_0< t_1<\dots <t_{n}= b$. Функция $S_{m}(t)$ называется сплайном степени $m$ с узлами на сетке $\{t_i\}$, если

  • на каждом отрезке $[t_i,t_{i+1}],\quad 0\le i <n$, функция $S_{m}(t)$ является многочленом степени $m$,
  • функция $S_{m}(t)$ имеет на отрезке $[a,b]$ непрерывные производные до $m-1$-го порядка включительно.

Если $m=1$, то сплайн представляет собой ломаную (на каждом интервале строится многогочлен первой степени), непрерывную в узлах сетки. При $m=2$ сплайн состоит из участков квадратичных парабол, непрерывных в точках сетки вместе с производной. При $m=3$ сплайн состоит из участков кубических парабол, а в точках сетки требуется непрерывность функции, ее первой и второй производных.

В общем случае сплайн $S_{m}(t)$ определяется на каждом сеточном интервале $m+1$ коэффициентами. Поскольку число интервалов равно $n$, то сплайн имеет $(m+1)n$ коэффициентов, для нахождения которых нужно столько же условий. Во внутренних узлах сетки $t_1, t_2,\dots, t_{n-1}$, согласно определению, непрерывна функция и $m-1$ призводных: $m$ условий в каждом узле, всего $m(n-1)$ условий. Если заданы значения функции в узлах (задача интерполяции), то получаем еще $n+1$ условие, т.е. всего их число $(m+1)n-(m-1)$. Недостающие $m-1$ условий могут быть заданы в в виде граничных условий. В частности, для сплайна первого порядка (ломаная линия, $m=1$) дополнительных условий не нужно, для сплайна второго порядка ($m=2$) требуется одно дополнительное условие, для кубического сплайна ($m=3$) нужно задать два граничных условия.

Аппроксимация. Сглаживающие кривые

Пусть, как и выше, в пространстве задано $n+1$ точек $P_0$, $P_1$, $P_2$, $\dots, P_n$, которые будем называть опорными. Задача аппроксимации состоит в построении кривой в каком-то смысле близкой к этим точкам, гладко аппроксимирующей ломаную, получающуюся при соединении точек отрезками в указанном порядке.

Будем строить сглаживающие кривые вида

(10)
\begin{align} {\bf r}(t)=\sum_{i=0}^n c_i(t){\bf r}_i \end{align}

где коэффициенты $c_i(t)$ не зависят от опорных точек. Выше мы уже обсуждали желательные свойства этих коэффициентов. При построении интерполяционной кривой условие $c_i(t_j)=\delta_{ij}$ не позволяет гладким коэффициентам $c_i(t)$ быть неотрицательными. Для сглаживающей кривой этого условия уже нет, так что мы будем требовать выполнения неравенства $c_i(t)\geq 0$. Как и выше, коэффициенты $c_i(t)$ будем строить в виде многочленов.

Далее будут описаны свойства двух типов сглаживающих кривых: кривой Безье и B-сплайна.

Интерполяция в OpenCascade

Процедуры класса GeomAPI_Interpolate обеспечивают построение интерполяционой кривой в виде B-сплайна по массиву точек. Инициализация данных выполняется созданием объекта класса, простейший конструктор которого имеет вид:

// Points       -- массив точек
// PeriodicFlag -- признак периодичности
// Tolerance    -- требуемая точность вычислений
   GeomAPI_Interpolate(
      const Handle(TColgp_HArray1OfPnt)& Points,
      const Standard_Boolean PeriodicFlag,
      const Standard_Real Tolerance);

Возможно также указание наклона кривой (вектора касательной) в отдельных точках.
Само построение интерполяционной кривой выполняется процедурой Perform(), проверить, прошло ли построение успешно, можно вызовом IsDone(), а получить ссылку на постоенную кривую можно посредством обращения к Curve():
// Построение BSpline-кривой, интерполирующей точки
   void Perform()

// Возвращается true, если построение выполнено успешно
   Standard_Boolean IsDone()

// Возвращается построенная BSpline-кривая
   Handle_Geom_BSplineCurve& Curve()

В качестве примера построим интерполяцию участка винтовой линии по пяти точкам.

// Определяем массив из пяти точек
   Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1,5);
// Точки винтовой линии с шагом по углу pi / 2
   double dphi = 3.14159/2.;
   double dz = 1. / 4.;
   int i;
   for( i=0; i<5; i++ )
      points->SetValue(i+1,gp_Pnt(cos(i*dphi),sin(i*dphi),i*dz ));
// Создаем объект класса
   GeomAPI_Interpolate interpolate(points,Standard_False,Precision::Approximation());
// Выполняем интерполяцию
   interpolate.Perform();
// Проверяем, все ли нормально
   if( ! interpolate.IsDone() ) {
   // Обработка ошибочной ситуации
   }
// Возвращается построенная BSpline-кривая
   Handle(Geom_BSplineCurve) сurve = interpolate.Curve();

Опросим параметры построенной кривой. Напомним, что B-сплайн кривая характеризуется

  • степенью;
  • числом узлов сетки параметра;
  • массивом значений узлов;
  • массвом значений кратностей узлов;
  • числом опорных точек;
  • массивом опорных точек.

Вот соответствующий участок программы:

//
// ----- Опрос параметров кривой
// Степень
   int degree = curve->Degree();   //   3
// Возвращает число узлов сетки параметра.
// Это длина массивов Knots, Multiplicities, Weights
    Standard_Integer nbKnots = curve->NbKnots();   // 5
// Возвращает массив значений узлов
    TColStd_Array1OfReal knots(1,nbKnots);
    curve->Knots(knots);   // ( 0.0000, 1.4361, 2.8722, 4.3084, 5.7445 )
// Возвращает массив кратностей узлов кривой
   TColStd_Array1OfInteger mults(1,nbKnots);
    curve->Multiplicities(mults);   // ( 4, 1, 1, 1, 4 )
// является рациональной
    Standard_Boolean isRational = curve->IsRational();   // false
// Возвращает число опорных точек.
// Это длина массива Poles
    Standard_Integer nbPoles = curve->NbPoles();   // 7
// Возвращает массив опорных точек кривой
    TColgp_Array1OfPnt poles(1,nbPoles);
    curve->Poles(poles);

Результат выполненного выше построения представлен на рисунке справа. Желтым цветом показан участок винтовой линии, красным — интерполяционная кривая. Точки интерполяции отмечены крестиком: в них обе кривые пересекаются. Зеленая ломаная соединяет опорные точки интерполяционного B-сплайна.

OCC_intrp01.jpg

Аппроксимация набора точек кривой в OpenCascade

Методы класса GeomAPI_PointsToBSpline позволяют построить B-слайн кривую, аппроксимирующую последовательность точек в пространстве. При создании объекта класса указываются параметры аппроксимации кривой, а именно:

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

Пример, аналогичный предыдущему

// Точки на винтовой линии
    TColgp_Array1OfPnt points(1,5);
    double dphi = 3.14159/2.;
    double dz = 1. / 4.;
    int i;
    for( i=0; i<5; i++ ) 
        points.SetValue(i+1,gp_Pnt(
            cos(i*dphi),sin(i*dphi),i*dz ));
// Аппроксмирующая кривая
// B-сплайн степени 3            
    Handle(Geom_BSplineCurve) curve = 
       GeomAPI_PointsToBSpline(points,3,3).Curve();
OCC_apprx01.jpg