Кривые и поверхности в компьютерной геометрии, I

Разбить на страницы

Показывать лекцию целиком

Введение: зачем и как нужно интерполировать

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

Рассмотрим простой пример. На первом курсе вам рассказывали, как интерполировать график функции одного вещественного переменного с помощью интерполяционного многочлена. Для этого, напомним, надо выбрать на графике (n + 1) точку и построить единственный многочлен степени n, график которого через эти точки проходит. При этом качество полученной интерполяции существенно зависит от расположения выбранных точек. Если выбрать точки на графике, например равномерно, то отклонение итогового многочлена от приближаемой функции может расти с ростом степени полинома:

In[1] := DynamicModule {RungePlot, г, Runge, x} ,
                 Manipulate[RungePlot[n],  
                   {{n, 15, "количество точек"}, 5, 50, 5, 
                             Appearance -> "Labeled" , ControlPlacement -> Bottom} , 
                         SaveDefinitions -> True] ,
                       Initialization: -> ( r[x_] := 1/(1+25x²);
                          Runge[n_] := Table[
                             {x, r[x]}, {x, -1., 1., 2. / (n-1)}];
                          RungePlot[n_] :=Module[{data, interp},
                              interp = Expand[
                                 InterpolatingPolynomial[Runge[n] , x] ] ;
                           Column[{Plot[{1/ (l + 25x^2), interp}, {x, -1, 1},
                               ImageSize -> {300} , Epilog -> {
                                    PointSize [Medium] , Red, Point/@ Runge [n] } ,
                                 PlotStyle-> {Orange, Blue}, PlotRange-> {-1, 1},
                                 Ticks -> {{-1, 1} , {-1, 1}} ,
                                 ImageSize -> 300, PlotRangePadding ->  .2] ,
                             ListLinePlot[Table[{x, Abs[N[N[interp] - r[x]]]},
                                {x, -1, 1, .01}] ,
                               ImageSize -> 300, PlotRange -> All,
                             Ticks -> { { -1, 1} , Automatic} ,
                             PlotLabel -> "Ошибка интерполяции" ,
                             PlotStyle -> Black, AspectRatio -> .25,
                             ImageSize -> 300] } , Center] ] ) ]

Проиллюстрированный здесь эффект носит название феномена Рунге.

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

В дословном переводе слово spline означает лекало или рейку используемую в строительных или морских конструкциях. Впервые понятие сплайна как кусочно-полиномиальной функции было введено Исааком Шенбергом (1903-1990, Румыния, США) в 1943-1946 годах при решении прикладных задач для американских ВВС, а затем и в теоретических исследованиях. Широкую известность и применение сплайны нашли после того, как в 60-е годы прошлого века французские математики-прикладники Пьер Безье, работавший на Рено, и Поль де Кастелье, работавший на Ситроен, независимо применили их для создания первых САПР в автомобилестроении. Конечно, мы не знаем точно, как выглядела их программа, но суть ее, вероятно, сводилась к следующей демо-версии:

In[2] :=DynamicModule [ {color, carbody, degrees, mesh}, 
               Manipulate[ 
                 Graphics3D[{
                    color, Specularity[White, 50] ,
                    BSplineSurf ace [carbody , SplineDegree -> Floor [degrees] ] , 
                    If[mesh, 
                      {Dashed, Blue, Line[carbody], Line[Transpose[carbody]], 
                        Red, PointSize [Medium] , Point /@ carbody} , { } ] } , 
                    Boxed -> False, Lighting -> "Neutral" , RotationAction -> "Clip" , 
                    Viewpoint -> {-2,-2,2}], 
                Row[ 
                    {Control [{{degrees, {2, 2}, "Порядок"}, {1, 1}, {14, 7}, 
                        ImageSize-> {140, 70}}], 
                  Dynamic[Floor[degrees]]}, Spacer[10]], 
               Row[ 
                   {Control[{{color, Yellow, "Цвет"}, 
                       ColorData["HTML", "GoldenRod"]}], 
                    Control[{{mesh, True, "Сетка"}, {False, True}}]}, Spacer[10]]  
                 ] , Initialization :-> (carbody = {
                    {{-.2, -2.25, 0}, {-.2, -2.25, 0}, {-.2, -2, 0}, {-.2, 0, 0},
                    {-.2, 2, 0}, {-.2, 2.25, 0}, {-.2, 2.25, 0}},  
                   (* \dots  еще тринадцать семерок точек  \dots  *) 
                     {{11.01, -2.25, 0}, {11.01, -2.25, 0}, {11.01, -2, 0},
                     {11.01,0,0}, {11.01,2,0}, {11.01,2.25,0}, {11.01,2.25,0}}
                  })

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

На первых нескольких занятиях этого семестра мы познакомимся с базовыми понятиями современной теории сплайнов. Разумеется, у нас нет времени сколько-нибудь углубляться здесь в подробности. Современная теория сплайнов представляет собой хорошо разработанную, разветвленную и сложную науку. Заинтересовавшиеся слушатели приглашаются на спецкурс кафедры Дифференциальной геометрии и приложений Компьютерная геометрия, который читают Д.П. Ильютко и Г.В. Носовский. См. также книгу Н. Н. Голованова, Д. П. Ильютко, Г. В. Носовского, А. Т. Фоменко, Компьютерная геометрия [1].

Мы начнем с одномерного случая.

Сплайны

Определение 5.2.1. Сплайном называется кривая, построенная по совокупности точек $$p_0,\dots ,p_n,$$ которые называются характеристическими или опорными точками сплайна. Требуется, чтобы сплайн r = r(t) проходил через заданные точки $$p_0, \dots ,p_n$$ при заданных значениях параметра $$t_0 \le \dots \le t_n: r(t_i) = p_i, 0 \le i \le n.$$ Значения параметра $$t_i (0 \le i \le п)$$ называются узловыми, а точки $$t_0, \dots ,t_n$$ на оси R - узлами сплайна.

Пример 5.2.1. Ломаная линия с вершинами в точках $$p_i$$:

$$r(t) = p_i(1-\omega)+p_i+1\omega,$$

где $$\omega=\frac{t-t_i}{t_{i+1}-t_i}$$ и $$t_i \le t \le t_{i+1}.$$ Параметр $$\omega$$ называется местным параметром на учacткe $$t_i \le t \le t_{i+1}.$$

Если $$t_i =i,$$ то параметризация ломаной называется равномерной. В этом случае параметрическая длина $$Т = t_n - t_0$$ ломаной равна $$п,$$ т. е. -числу опорных точек минус 1, если все опорные точки $$p_0, \dots, p_n$$ различны.

Если $$p_0 = p_n,$$ т. е. ломаная замкнута, ее параметрическая длина равна числу опорных точек (в предположении, что точки $$p_0, \dots , p_{n_1}$$ различны):

In[3] = DynamicModule [ {pts } , Manipulate [
                Show[Graphics[{Blue, Thick, Line[{pts}]}]], 
                 {{pts,   {{1,  1},   {2,  2.3},   {3, 2.7},   {4, 3},   {5,  4}}}, 
                  Locator}]]

Рассмотрим основные примеры сплайнов.

Сплайн Лагранжа

Пусть дана $$(n+ 1) $$ опорная точка $$p_i = (x_i,y_i) \in R^2, 0 \le i \le n.$$ Требуется найти полиномиальную по $$t$$ кривую $$r = r(t) = P^n(t) $$ такую, что $$r( t_i) = p_i, 0 \le i \le n,$$ при заданных значениях параметра $$t_0 \le \dots \le t_n.$$ Будем искать решение этой задачи в виде

$$r(t)=\sum_{i=0}^{n}L(t_i)p_i,$$

где

$$L_i(t_j)= \delta_{ij}= \begin{cases} 1, i=j,\\ 0, i \ne j, \end{cases} 0 \le i,j \le n$$

Таким свойствами обладают многочлены

$$L_i(t)=\frac{(t-t_0) \dots (t-t_{i-1})(t-t_{i+1})\dots (t-t_n)}{(t_i-t_0) \dots (t_i-t_{i-1})(t_i-t_{i+1}) \dots (t_i-t_n)}=\frac{\Pi_{j=0, i \ne j}^{n}(t-t_j)}{\Pi_{j=0, j \ne i}^{n}(t_i-t_j)}$$

С их помощью получаем сплайн Лагранжа:

$$r(t)=\sum_{i=0}^{n}\frac{\Pi_{j=0, i \ne j}^{n}(t-t_j)}{\Pi_{j=0, j \ne i}^{n}(t_i-t_j)}p_i$$

или

$$r(t)=\sum_{i=0}^{n}\frac{w_{0,n}(t)}{(t-t_0)w_{0,n}'(t_i)}p_i,$$

где

$$w_{0,n}(t) = (t-t_0) \dots (t-t_n),\\ w_{0 ,n}' (t_i ) = \frac{dw_{0,n}}{ dt }|_{t=t_i}= ( t_i - t_0) \dots ( t_i - t_{i-1})( t_i - t_{i+1}) \dots ( t_i - t_n).$$

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

Еще один существенный недостаток сплайнов Лагранжа состоит в том, что они могут выходить за пределы выпуклой оболочки своих опорных точек, причем довольно значительно.

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

Пример 5.2.2. Сплайн Лагранжа с возможностью управления опорными точками:

In[4] :=
         DynamicModule [ {pts , ptsO , L, t, p, n = 6, tt} , 
           ptsO ={{6.0, 0.0}, {0.0, 0.0}, {4.8, 5.2}, {0.3, 5.5}, {9.0, 3.0},
              {4.3, 8.2} , {10.0, 6.0}} ;
           tt= {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}; 
           p[i_, s_] : =
             Apply [Times , (s - tt [ [tt] ] )  /@ Drop [Range [n + 1] , {i + 1} ] ] ; 
           L = Table [p [i , t] / p [i , tt [ [i + 1] ] ] , {i , 0 , n} ] ; 
           Manipulate[ 
             Show[{ParametricPlot[L.pts, {t, tt[[l]], tt[[-l]]}, 
                  PlotRange -> {{-15, 15}, {-15, 15}}], 
               Graphics [{Red, Text [ToString [tt - 1] , pts [ [#] ] +{1.2, 0}] /@ 
                   Range[n + 1] , Text["Сплайн Лагранжа = сплайн Ньютона", 
                   {8, -14}]}]}], {{pts, ptsO}, Locator}]]

Сплайн Ньютона

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

А именно, сплайн Ньютона записывается в виде

$$r(t) = a_0 + a_1(t - t_0) + a_2(t - t_0)(t - t_1) + \dots + a_n(t - t_0) \dots (t - t_{n_1}),$$

где векторные коэффициенты $$a_0, \dots , a_n$$ ищутся из условий прохождения сплайна через точки $$p_0, \dots , p_n.$$ При этом точки $$p_0, \dots , p_k$$ полностью определяют коэффициенты $$a_0, \doys , a_k, 0 \le k \le n,$$ с помощью системы уравнений

$$\begin{cases} p_0 = a_0,\\ p_1=a_0+a_1(t_1-t_0),\\ p_2 = a_0 + a_1(t_2 - t_0) + a_2(t_2 - t_0)(t_2 - t_1),\\ ...........\\ p_n = a_0 + a_1(t_n - t-0) + a_2(t_n - t_0)(t_n - t_1) + \dots +a_n(t_n-t_0) \dots (t_n-t_{n_1}) \end{cases}$$

Эта система решается за один проход, так как матрица системы является нижнетреугольной. Из первых $$1 \le k \le n + 1$$ уравнений определяются первые $$k$$ неизвестные $$a_0, \dots , a_{k_1}.$$ Поэтому если мы добавим еще одну точку $$p_{n+1},$$ нам будет нужно лишь определить еще один коэффициент $$a_{n+1}$$ из еще одного, добавленного уравнения. Коэффициенты $$a_0, \dots ,a_n$$ при этом останутся прежними.

Задача 5.2.1. Написать программу для сплайна Ньютона с возможностью добавлять точку.

Сплайн Эрмита

Сплайны Эрмита предназначены для соединения двух заданных кривых новой кривой, полиномиальной по параметру, так, чтобы получившаяся совокупная кривая имела заданную степень гладкости в точках стыковки. Отсюда следует, что сплайн Эрмита задается двумя своими крайними точками и заданным числом m производных радиус-вектора в каждой из крайних точек. Сплайн Эрмита имеет вид полинома степени $$(2m + 1) $$ по $$t.$$ Степень полинома называется порядком сплайна. При этом $$(2m + 2) $$ векторных коэффициентов этого полинома определяются из системы $$(2m+ 2) $$ уравнений ( $$m + 1$$ условие для каждой из двух крайних точек).

Сплайн Эрмита может быть определен формулой

$$r(t)=\sum_{i=0}^m(g_i(t,m)r_0^i+h_i(t,m)r_1^i),$$

где $$r_0^i, r_1^i, 0 \le i \le m,$$ - заданные значения радиус-вектора кривой и его $$m$$ производных на концах, а $$g_i(t, m), h_i(t, m) $$ даются формулами

$$\begin{cases} g_i(t,m)=\frac{1}{i!}\sum_{j=0}^{m-i}C_{2m+1-i}^{j}t^{i+j}(1-t)^{2m+1-(i+j)},\\ h_i(t,m)=(-1)^i \frac{1}{i!}\sum_{j=0}^{m-i}C_{2m+1-i}^{j}t^{2m+1-(i+j)}(1-t)^{i+j} \end{cases}$$

Задача 5.2.2. Написать программу для сплайна Эрмита с возможностью добавлять точку. Запрограммировав сплайны Эрмита, проверить на Mathematica, что производные радиус-вектора на концах сплайна равны заданным значениям.

Составной сплайн Эрмита

Составной сплайн Эрмита предназначен для построения кривой заданной степени гладкости, проходящей через заданную последовательность опорных точек. Он получается путем соединения соседних опорных точек сплайнами Эрмита одного и того же порядка. При этом значения производных радиус-вектора в опорных точках (требуемые для построения составляющих сплайнов Эрмита) необходимо как-то доопределять. Это можно сделать по-разному, исходя из различных соображений. Простые способы доопределения производных приводят к простым вычислениям, однако при неравномерном расположении опорных точек такой сплайн может образовывать необоснованные петли и изгибы. Более сложные способы доопределения производных (псевдоупругие сплайны Эрмита, например) дают хорошие результаты, но приводят к большим вычислениям.

Пример 5.2.3. Составной сплайн Эрмита порядка 3 ( m = 1 ). Для построения такого сплайна необходимо предварительно доопределить первые производные радиус-вектора в опорных точках $$p_0, \dots ,p_n.$$ Для промежуточных и крайних опорных точек способ доопределения различен.

Промежуточные опорные точки. В промежуточных опорных точках можно, например, положить производные $$q_i = r'( t_i) $$ равными

$$q_i=\frac{p_{i+1}-p_{i-1}}{t_{i+1}-t_{i-1}}, 1 \le i \le n-1$$

т. е. рассмотреть аппроксимацию производных.

Однако при неравномерном расположении опорных точек такой способ доопределения производных может привести к появлению необоснованных петель. Поэтому обычно пользуются другими схемами доопределения $$q_i.$$ Например, можно положить

$$q_i=s_{i+1}\frac{p_1-p_{i-1}}{s_i+s_{i+1}}+s_i \frac{p_{i+1}-p_i}{s_i+s_{i+1}}$$

или

$$q_i=s_i \frac{p_i-p_{i-1}}{s_i+s_{i+1}}+s_{i+1} \frac{p_{i+1}-p_i}{s_i+s_{i+1}},\\ 1 \le i \le n-1$$

где $$s_i = ||p_i - p_{i-1}||$$ - расстояние между оседними опорными точками.

Крайние опорные точки. В двух крайних опорных точках $$p_0, p_n,$$ если в них не требуется гладко стыковать сплайн с другими, уже заданными кривыми, обычно пользуются условием равенства нулю третьих производных радиус-вектора на концах (условие свободных концов кривой): $$r^{(3)}(t_0) = 0, r^{(3)}(t_n) = 0.$$ Можно показать, что соответствующие значения первых производных на концах составного сплайна Эрмита третьего порядка даются формулами

$$q_0 = 2(p_1 - p_0) - q_1, \\ q_n = 2(p_n - p_{n-1}) - q_{n-1}.$$

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

Пример 5.2.4. Составной сплайн Эрмита третьего порядка с шестью управляемыми мышью опорными точками:

In[5]:=
       DynamicModule Г {pts , ptsO, n,  t, tt, i,  j, S, gO, gl, h0, hi, w, q, r} ,  
        pts0 = {{6.0,  -5.0},   {0.0,  0.0},   {4.8, 5.2},   {0.3,  7.5},   {6.8,  12.2},   {-2.3,  0.0}}; 
           tt= {1.0, 2.0,  3.0,  6.0,  7.0,  8.0}; 
           n = Length[ptsO] - 1; 
           Manipulate [ 
              S = Table [{r[i,  t, pts] ,  tt [ [i] ]  <=. t <== tt [ [i + 1] ] } ,   {i,  1, n} ] ; 
             Show[{ParametricPlot[
                 Piecewise[S] ,   {t,  tt[ [1] ] ,  tt [ [-1] ] } , PlotRange -> {{-10, 20},   {-10, 20}}
               ],
              Graphics [{Red, Text[ToString[* - 1] , pts [ [*] ] +{1.2, 0}]  /@Range [n + 1] , 
                       Text["Составной сплайн Эрмита 3 порядка",  {7, -9}]}]}], 
               {{pts, ptsO} , Locator}] ,  Initialization: -> (
                  g0 [x_]   : = 2 x³ - 3 x² + 1; gl [x_] : = x³ - 2 x² + x; h0 [x_]: = 3 x² - 2 x³ ; 
                  h1 [x_]   : = x³ - x² ;
                  w[i_,  t_]   :=  (t-tt[[i]]) / (tt[[i + l]] -tt[[i]]); 
                  q[i_, pts_] : = 
                     If[i==l, 2 (pts [[2]] -pts[[l]]) - (pts[[3]] -pts[[l]]) / (tt[[3]] -tt[[l]]), 
                        If[i==n + 1, 2 (pts [ [n + 1] ] - pts [ [n] ] ) -
                             (pts[[n + 1]] -pts[[n-l]]) / (tt[[n + l]] -tt[[n-l]]) , 
                           (pts[ [i + 1]] -pts[[i-l]]) / (tt[[i + l]] -tt[[i-l]])]]; 
                     r [i_,   t_, pts_]   : = gO [w[i,   t] ] * pts[ [i] ] + gl [w[i,   t] ] *q[i, pts] + 
                        h0[w[i,  t] ] *pts[ [i + 1] ] +hl[w[i,  t] ] *q[i + 1, pts] ) ]

Псевдоупругие сплайны Эрмита

При построении составного сплайна Эрмита третьего порядка неизвестные производные $$q_0, \dots ,q_n$$ в точках $$p_0, \dots ,p_n$$ могут быть определены также из следующих соображений.

Представим себе, что через точки $$p0, \dots ,p_n$$ должна пройти упругая нить, которая стремится занять положение, соответствующее минимуму запасенной энергии ее внутренних сил упругости. Эта энергия, как известно из физики, выражается интегралом от квадрата кривизны k(s) кривой, по которой пролегает нить:

$$E=\frac12 \int_{t_0}^{t_n}c(s)k^2(s)ds,$$

где s - натуральный параметр на кривой, а c(s) - коэффициент упругости нити. Мы будем считать, что c(s) = c = const.

Предположим, что искомый сплайн задается графиком функции y = y(x). Поскольку минимизация энергии упругих сил (5.3) является, вообще говоря, достаточно трудоемкой вычислительной задачей, мы заменим ее на более простую, но достаточную для наших целей задачу. А именно, заменим энергию E, задаваемую формулой (5.3), на псевдоупругую энергию, заменив кривизну k на модуль второй производной $$\left | \frac{d^2y}{dx^2} \right |.$$ Это приводит к задаче минимизации интеграла

$$E=\frac12 \sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}(y''(x))^2dx$$

за счет выбора величин $$y_0', \dots, y_n',$$ где $$y_i'=y'(x_i).$$ Условие стационарной точки целевой функции (5.4), получаемое приравниванием к нулю ее частных производных, имеет вид

$$\frac{\partial E}{\partial y_j'}=\sum_{i=0}^{n-1} \int_{x_i}^{x_{i+1}}y''(x) \frac{\partial y''(x)}{\partial y_j'}dx=0, j=0, \dots, n$$

Пусть точки $$p_0, \dots , p_n$$ на $$R^2$$ имеют координаты $$(x_0, y_0), \dots , (x_n, y_n) $$ соответственно. Будем считать, что $$x_i < x_{i+1}$$ для любого $$i, 0 \le i \le n - 1.$$ Введем обозначения: $$\Delta i = x_{i+1} -x_i.$$ Тогда местный параметр $$w$$ на участке $$(x_i,x_{i+1})$$ будет равен $$w = \frac{x- x_i}{\Delta _i}.$$ На отрезке $$[x_i,x_{i +1}] $$ обозначим $$y(x) $$ через $$y_{i, i+1}(x).$$

Согласно (5.1), на участке $$(x_i, x_{i+1}), i = 0, \dots , n - 1,$$ зависимость

$$\tilde y_{i,i+1}(w)=y_{i,i+1}(x(w))$$

имеет вид

$$\tilde y_{i, i+1}(w)=y_ig_0(w)+y_{i+1}h_0(w)+y_i' \Delta_ig_1(w)+y_{i+1}' \Delta_ih_1(w)$$

поскольку $$y_i'=\frac{dy_{i,j+1}}{dx} \left |_{x=x_i}= \frac{1}{\Delta_i} \frac{d \tilde y_{i, j+1}}{dw} \left |_{x=x_i}$$ и $$y_{i+1}'=\frac{dy_{i,j+1}}{dx} \left |_{x=x_{i+1}}=\frac{1}{\Delta_i} \frac{d \tilde y_{i,j+1}}{dw} \left |_{w=1}.$$

Подставляя (5.6) в (5.5), интегрируя по частям, пользуясь явным видом функций $$g_0( w ), h_0( w ), g_1( w ), h_1(w ) $$ из (5.2), а также тем, что $$\frac{d^4y}{dx^4},$$ приходим к следующей системе уравнений:

$$\begin{cases} y_{i-1,i}''(x_i)=y_{i,i+1}''(x_i), i=1, \dots, n-1\\ y_{0,1}''(x_0)=y_{n-1,n}''(x_n)=0 \end{cases}$$

которая обозначает, что вторые производные $$y??(x) ,$$ а следовательно, и кривизна кривой, должны быть непрерывны в точках стыковки $$p_1, \dots ,p_{n-1},$$ а на концах $$y?$$? обращается в нуль.

Подставим теперь (5.6) в (5.7) и снова воспользуемся явным видом функций $$g_0(w), h_0(w), g_1(w), h_1(w) $$ из (5.2). Используя равенства $$y_{i,j+1}''(x_i)=\frac{1}{\Delta_i^2} \frac{d \tilde y_{i,j+1}}{dw} \left |_{w=0} $$ и $$y_{i,j+1}''(x_{i+1})=\frac{1}{\Delta_i^2} \frac{d \tilde y_{i,j+1}}{dw} \left |_{w=1}$$ получим систему уравнений на $$y_i'$$ c ленточной 3-шаговой матрицей:

$$\begin{cases} y_{i-1}' \frac{1}{\Delta_{i-1}}+2y_i' \left ( \frac{1}{\Delta_{i-1}} + \frac {1}{\Delta_i} \right )+y_{i+1}'\frac{1}{\Delta_i}=3 \frac{y_i-y_{i-1}}{\Delta_{i-1}^2}+3 \frac{y_{i+1}-y_i}{\Delta_i^2},\\ 1 \le i \le n-1,\\ 2y_0' \Delta_0 + y_i' \Delta_0=3(y_1-y_0) (\mbox{граничное условие в точке} p_0\\ y_{n-1}' \Delta_{n-1}+2y_n' \Delta_{n-1}=3(y_n-y_{n-1})(\mbox{ граничное условие в} p_n) \end{cases}$$

Это - система из ( n + 1 )-го уравнения с ( n + 1 ) неизвестным $$y_0', \dots, y_n'.$$ Поскольку матрица системы имеет ленточный вид, то она допускает численное решение за линейное по n время.

Теперь рассмотрим более общий случай параметрически заданного сплайна r = r(t) (т. е. не обязательно в виде графика функции y = y( x) ).

На участке между точками $$p_i$$ и $$p_{i+1}$$ обозначим искомый сплайн Эрмита порядка 3 через $$r_i (w) (0 \le i \le n - 1) ,$$ где $$0 \le w \le 1$$ - местный параметр. Имеем, согласно (5.1),

$$r_i(w)=p_ig_0(w)+p_{i+1}h_0(w)+\dot r_i(0)g_1(w)+\dot r_i(1)h_1(w),\\ i=0, \dots, n-1$$

К этому надо добавить условия непрерывности первых производных в точках стыковки:

$$\dot r_{i-1}(1)=k_i \dot r_i(0), i=1, \dots, n-1$$

где $$k_i =\frac{t_i-t_{i-1}}{t_{i+1}-t_i}.$$ Будем считать, что $$\frac{|t_i-t_{i-1}|}{|t_{i+1}-t_i|}=\frac{|p_i-p_{i-1}|}{|p_{i+1}-p_i|}$$ т. е. параметр $$t$$ пропорционален натуральному. Тогда и местный параметр $$w$$ будет также пропорционален натуральному. Обозначим $$a_i = p_{i+1} - p_i, |a_i| = a_i.$$ Из пропорциональности параметра $$w$$ натуральному следует, что $$ds = a_idw$$ на участке $$(t_i,t_{i+1}).$$ Следовательно, кривизна $$k$$ пропорциональна $$\frac{\ddot r_i}{a_i^2},$$ а потому псевдоупругая энергия сплайна Эрмита может быть определена как

$$E=\sum_{i=0}^{n-1} \frac{1}{a_i^3}\int_0^1(\ddot r_i(w))^2dw)$$

При этом, согласно (5.8) и (5.9), $$\ddot r_i(w)$$ выражается через $$\dot r_i(0) $$ и $$\dot r_{i+1}(1) .$$ Условие стационарности (необходимое условие экстремума функции (5.10)) имеет вид

$$\frac{\partial E}{\partial \dot r_i(0)}=0, i=0, \dots, n-1,\\ \frac{\partial E}{\partial \dot r_{n-1}(1)}=0$$

Подставляя (5.8) и (5.9) в (5.11), приведем (5.11) к виду

$$\begin{cases} \ddot r_{i-1}(1)=k_i^2 \ddot r_i(0), i=1, \dots, n-1\\ \ddot r_0(0)= \ddot r_{n-1}(1)=0 \end{cases}$$

где $$k_i=\frac{a_{i-1}}{a_i}.$$ Подставим сюда (5.8). Получим после алгебраических преобразований следующее выражение:

$$\begin{cases} \dot r_{i-1}+2 \dot r_i(k_i+k_i^2)+\dot r_{i+1}k_i^2k_{i+1}=\\ =3(p_i-p_{i-1})+3(p_{i+1}-p_i)k_i^2, 1 \le i \le n-1,\\ 2 \dot r_0+\dot r_ik_1=3(p_1-p_0),\\ \dot r_{n-1}+2 \dot r_n=3(p_n-p_{n-1}), \end{cases}$$

где $$\dot r_i= \dot r_i(0), i=1, \dots, n-1, \dot r_n=\dot r_{n-1}(1).$$ Эта система состоит из $$(n + 1) $$ уравнений c $$(n+ 1) $$ неизвестными $$\dot r_0, \dots , \dot r_n.$$ Как и выше, матрица системы является ленточной, т. е. ненулевые элементы в ней расположены лишь в узкой окрестности главной диагонали. Такие системы уравнений допускают численное решение за линейное время с помощью хорошо разработанных алгоритмов.

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

Кубические сплайны

Пусть дана $$(n + 1) $$ точка $$p_0, \dots ,p_n \in R^2.$$ Требуется провести через эти точки кривую класса $$C^2,$$ т. е. имеющую непрерывные первые и вторые производные радиус-вектора по параметру: $$r = r(t), r(t_i) = p_i, i = 0, \dots , n,$$ и $$\dot r(t), \ddot r(t) $$ - непрерывны по $$t.$$ Решение этой задачи можно дать в виде кубического сплайна. На каждом участке $$[t_i,t_{i+1}]$$ такой сплайн описывается полиномом третьей степени по $$t.$$ При полиномиальной интерполяции степень полинома в данной задаче не может быть, вообще говоря, сделана меньше 3. Поэтому кубический сплайн имеет наименьшую степень в данной постановке задачи и определяется однозначно.

Пусть искомый сплайн задается параметрическим уравнением $$r = r(t) .$$ Обозначим $$s_i = \ddot r(t_i ), i = 0, \dots , n.$$ Поскольку сплайн является кубическим, вторая производная радиус-вектора будет линейной функцией на каждом отрезке $$[t_i,t_{i+1}]$$:

$$\frac{d^2r}{dt^2}=s_i \frac{t_{i+1}-t}{t_{i+1}-t_i}+s_{i+1}\frac{t-t_i}{t_{i+1}-t_i}.$$

Проинтегрируем (5.12) два раза. Получим

$$r(t)=s_i \frac{(t_{i+1}-t)^3}{6(t_{i+1}-t_i)}+s_{i+1}\frac{(t-t_i)^3}{6(t_{i+1}-t_i)}+c_1t+c_2$$

Постоянные $$c_1$$ и $$c_2$$ определим из условий, заданных на концах данного участка кривой:

$$\begin{cases} r(t_i)=p_i,\\ r(t_{i+1})=p_{i+1} \end{cases}$$

Подстановка (5.13) в (5.14) и вычисление $$c_1, c_2$$ приводит к следующему ответу:

$$r(t)=s_i \frac{(t_{i+1}-t)^3}{6(t_{i+1}-t_i)}+s_{i+1} \frac{(t-t_i)^3}{6(t_{i+1}-t_i)}+\\ +\left( \frac{p_i}{t_{i+1}-t_i}- s_i \frac{t_{i+1}-t_i}{6}\right)(t_{i+1}-t)+\\ + \left( \frac{p_{i+1}}{t_{i+1}-t_i}-s_{i+1} \frac{t_{i+1}-t_i}{6}\right) (t-t_i)$$

Правая часть (5.15) является кубическим полиномом по $$t$$ на отрезке $$[t_i,t_{i+1}].$$ Она содержит два неизвестных векторных параметра $$s_i$$ и $$s_{i+1}.$$ Их можно определить из условия согласованности первых производных на концах отрезков $$[t_i,t_{i+1}].$$ Пока мы еще не использовали этого условия (вторые производные уже согласованы в силу определения si). На левом конце отрезка $$[t_i,t_{i +1}] $$ мы имеем

$$\frac{dr}{dt} \left |_{t=t_i}=- \frac{(2s_i+s_{i+1})(t_{i+1}-t_i)}{6}+\frac{p_{i+1}-p_i}{t_{i+1}-t_i}$$

Аналогично, на правом конце отрезка $$[t_{i-1},t_i] ,$$ т.е. в той же точке $$t_i,$$ имеем

$$\frac{dr}{dt} \left |_{t=t_i}= \frac{(2s_i+s_{i-1})(t_i-t_{i-1})}{6}+\frac{p_i-p_i-1}{t_i-t_i-1}}$$

Приравнивая правые части (5.16) и (5.17), получим $$(n - 1) $$ уравнение для внутренних точек $$i = 1, \dots ,n-1$$:

$$s_{i-1}(t_i-t_{i-1})+2s_i(t_{i+1}-t_{i-1})+s_{i+1}(t_{i+1}-t_i)=6\frac{p_{i+1}-p_i}{t_{i+1}-t_i}-6 \frac{p_i-p_{i-1}}{t_i-t_{i-1}},\\ i=1, \dots, n-1$$

Система (5.18) представляет собой систему из $$(n - 1) $$ векторного уравнения (или, что то же самое, $$2(n - 1) $$ скалярного уравнения) с ленточной трехдиагональной матрицей и $$(n+1) $$ векторных неизвестных $$s_0, \dots , s_n$$ (что соответствует $$2(n+ 1) $$ скалярных неизвестных). Если кривая замкнута, то $$s_0 = s_n$$ и в систему (5.18) добавляется еще одно уравнение для $$i = n$$:

$$s_{n+1}=s_1.$$

Таким образом, для замкнутой кривой число неизвестных равно числу уравнений и система имеет единственное решение.

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

$$s_0 = s_n=0.$$

Еще один возможный вариант граничных условий:

$$s_0 = s_1, s_n = s_{n-1}$$

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

Определив из системы (5.18) и граничных условий вида (5.19) или (5.20) (или каких-либо других, следующих из постановки задачи) значения $$s_0, \dots , s_n,$$ получим формулу для искомого кубического сплайна на участке $$[t_i, t_{i+1}]$$:

$$r(t) = (1 - w)p_i + wp_{i+1} +((-2w + 3w^2 - w^3 )s_i + (-w + w^3 )s_{i+1} )\frac{(_{ti+1}-t_i)^2}{6},$$

где $$w =\frac{t-t_i}{t_{i+1}-t_i}$$ - местный параметр на участке $$[t_i,t_{i+1}].$$

Замечание 5.2.1.

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

Если взять равномерную параметризацию, т. е. положить $$t_i = i, i = 0, \dots ,n,$$ то уравнения (5.18) примут совсем простой вид:

$$s_{i-1} + 4s_i + s_{i+1} = 6(p_{i+1} - 2p_i + p_{i-1}), i = 1, \dots , n - 1.$$

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

$$\frac{|p_{i+1}-p_i|}{t_{i+1}-t_i} \approx \frac{|p_i-p_{i-1}|}{t_i-t_{i-1}}.$$

Сплайны в пакете Mathematica. Для создания сплайнов используется подпакет SplineFit, просто надо загрузить Splines' из Splines Package. Для этого используется функция SplineFit[data, type] , которая возвращает сплайн для данных data, используя сплайн-аппроксимацию типа type - по умолчанию это кубический сплайн Cube (другие типы - Bezier и CompositeBezier ).

Пример 5.2.5. Кубический сплайн с шестью управляемыми мышью опорными точками:

In[6] : = DynamicModule [ {pts , pts0, n = 6, spline},
                   Pts0= {{6.0, -5.0}, {0.0, 0.0}, {4.8, 5.2}, {20.3, 10.5}, {6.8, 12.2}, 
                         {20.0, 20.0}}; 
                   Manipulate[Show[{
                            ParametrioPlot[spline[pts][x], {x, 0, 6}, 
                              PlotStyle -> {Thick}] , Graphics [ {Red, 
                                Text[ToString[#- 1] ,
                                    pts[[#]] + {1.4, 0}]  /@ Range [n] }} } , 
                          PlotRange ->{{-l,23}, {-6, 21}}], 
                      {{pts, ptsO}, Locator}], 
                   Initialization : -> ( 
                        Needs["Splines'"]; 
                         spline[pfcs_] := SplineFit[pts, Cubic])]

Кривые Безье

Обычные кривые Безье

Поставим следующую задачу. Для любого набора из (n+ 1) опорных точек $$p_0,\dots, p_n \in R^3$$ требуется найти такие функции $$f_0(t), \dots , f_n(t) ,$$ не зависящие от $$p_0, \dots , p_n$$ и имеющие вид полиномов по $$t,$$ что кривая

$$r(t)=p_of_o(t)+ \dots + p_nf_n(t), t \in [0,1]$$

обладает следующими свойствами.

$$r(0)=p_0,r(1)=p_n.$$

Кривая (5.21) не зависит от выбора декартовой системы координат в $$R^3$$ или, другими словами, формула (5.21) инвариантна относительно движений пространства.

Если все опорные точки $$p_0, \dots , p_n$$ лежат в одной плоскости, то кривая (5.21) плоская и лежит внутри выпуклой оболочки опорных точек $$p_0, \dots ,p_n.$$

Теорема 5.1. Условия 1 - 3 равносильны условиям

$$\sum_{i=0}^{n} f_i (t) = 1$$ для любой точки $$t \in [0,1] $$ ;

$$f_i(t) \ge 0$$ для любого $$t \in [0,1] $$ ;

$$f_0(0) = 1, f_i (0) = 0$$ для любого $$1 \le i \le n$$ и $$f_n(1) = 1, f_i(1) = 0$$ для любого $$0 \le i \le n-1.$$

Определение 5.3.1. Фиксируем целое $$n > 0.$$ Базисом Бернштейна на отрезке [0,1] называется система из (n+1) полинома степени n, задаваемых формулами

$$B_i^n(t)=\frac{n!}{i!(n-i)!}t^i(1-t)^{n-i}=C_n^it^i(1-t)^{n-i}, 0 \le i \le n$$

Следующая программа рисует полиномы Бернштейна для $$0 \le n \le 20$$ и для всех $$k \le n.$$

Пример 5.3.1. Полиномы Бернштейна:

In[7]:=
       DynamicModule[{i, k, n, DivZ, PowerZ, Bernstein},
          DivZ[s_, 0] := 0; 
       DivZ[s_, t_] := s/ t;
          PowerZ[(x_)?((#1 ==0) ), (y_)?((#2 ==0) )] := 1; 
       PowerZ[x_,   y_]    :=  x^y; 
         Bernstein [k_,   n_,    t_] : = Binomial [n, k]  * PowerZ [1 - t, n - k]  * PowerZ [t, k] ; 
        Manipulate [Module [ {ps , f} , If[i > n, i = n] ;
            ps = Table [ {Hue [. 67 , .6, . 6] , Opacity!. 1 (8 + 2 Boole [i == k]) ] ,
                Thickness [0.001 (1 + 4 Boole [k == i]) ] } , {k, 0 , n} ] ; 
            f[t_] := Table [Bernstein [k, n, t] , {k, 0, n}];
            Show[piot [Evaluate@f [x] , {x, 0, 1}, PlotStyle -> ps , PlotPoints -> 25 , 
               PlotRange ->  {0 , 1.1},
               Filling -> If [fill, {i + 1 ->  {Axis, {Opacity [. 15] , Green}}}, None], 
               PlotLabel -> Style [B^i,n[t] == Bernstein [i , n, t] , 14]], 
           Graphics [ If [n ==0, {} , With[{j = i + .3 Boole [i == 0] - . 3 Boole [i == n] } , 
                   Arrow[{{.5, 1.1} , {DivZ[j, n] , Bernstein[i, n, DivZ[j, n]]}}]]]] , 
           ImageSize -> {450, 300}, ImagePadding -> {{25, 25}, {25, 38}}]], 
        {{n, 12, "Степень полинома n"} , 0, 20, 1}, 
        { {i, 5 , " Номер полинома i"}, 0, n, 1}, 
        {{fill, True, "Заполнение"}, {True, False}, ControlPlacement -> Bottom}]
]

Базис Бернштейна удовлетворяет свойствам a)-c) , сформулированным в теореме, и поэтому в качестве функций $$f_0(t), \dots ,f_n(t) $$ можно взять базис Бернштейна.

Определение 5.3.2. Пусть дана (n + 1) опорная точка $$p_0, \dots , p_n.$$ Тогда кривой Безье, определяемой точками $$p_0, \dots , p_n,$$ называется кривая

$$r_B(t)=\sum_{i=1}^nB_i^n(t)p_i, t \in [0,1].$$

Число $$n$$ называется степенью кривой Безье.

Замечание 5.3.1.

Кривая Безье является полиномиальной, а следовательно, бесконечно гладкой на всей области определения параметра [0,1].

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

Из определения кривой Безье непосредственно следует, что она проходит через свои крайние опорные точки при значениях параметра $$t = 0$$ и $$t = 1$$ соответственно. Однако через все промежуточные опорные точки кривая Безье, вообще говоря, не проходит. Таким образом, опорные точки являются лишь средствами управления кривой Безье.

В своих крайних опорных точках кривая Безье касается звеньев своей опорной ломаной (ломаной, построенной по опорным точкам).

Кривые Безье в пакете Mathematica. Mathematica имеет встроенную функцию построения кривой Безье BezierCurve[{p1, p2, p3, p4}], однако эта функция строит кривую Безье только по 4 опорным точкам. Если в качестве ее аргументов задать больше, чем 4 точки, то вместо кривой Безье будет построена составная кривая Безье, имеющая изломы в каждой четвертой опорной точке (попробуйте). Для того чтобы в Mathematica построить кривую Безье по более чем 4 опорным точкам, надо пользоваться функцией BezierFunction[pts] , где [pts] - список координат опорных точек. Эта функция задает зависимость радиус-вектора кривой Безье от параметра. Кривую Безье можно изобразить, применив к ней оператор ParametricPlot или ParametricPlot3D.

Пример 5.3.2. Кривая Безье, построенная по четырем опорным точкам (точки управляемы):

In[8] : = DynamicModule [{pts , pts0 , n = 4} ,
                 pts0 = {{-3, 0} , {0, 4}, {2, 4} , {3, 2}} ; 
                 Manipulate[ 
                    Show [ {Graphics[{BezierCurve[pts] , Green, Line[pts] ,
                           Red, Text[ToString[#- 1] , pts [ [#]] +{0.4, 0}]  /@Range [4]}] } , 
                      PlotRange -> {{-4, 5} , {-2,5}}], {{pts , ptsO} , Locator} ] 
              ]

Следующий пример демонстрирует пункт 4 замечания.

Пример 5.3.3. Кривая Безье, построенная по пяти опорным точкам, с управляемыми второй и предпоследней точками. Цель программы - показать, что кривая Безье всегда касается звеньев своей опорной ломаной в крайних точках.

In[9]:=
        DynamicModule[{u, v, t, pts}, 
          Manipulate [pts = {{-3, 0}, и, {0, 0}, v, {3, 2}}; 
             BezierF = BezierFunction[pts] ;
             Show [ Graphics [ {RGBColor [. 12 , .34, . 65] , Line [ { {-3 , 0} , u} ] , 
               Line[{v, {3, 2}}] , AbsolutePointSize[6] , Point[{u, v} ] ,
               RGBColor[l, .21, 0] , PointSize [Large] , Point[{{-3, 0}, {0, 0}, {3, 2}}]}, 
            Axes -> True] , 
          Graphics [{Red, Text [ToString[# - 1] , pts [ [#] ] +{0.3, 0}] /@Range[5],
              Text["Составной сплайн Эрмита 3 порядка", {7, -9}]}], 
            ParametricPlot[BezierF[t], {t, 0, 1},
              PlotStyle -> {{Red, AbsoluteThickness [2] } } ] , PlotRange -> 6, 
           ImageSize -> 6 x 72 , Background -> GrayLevel [. 98] ] , 
        {{u, {2, -3}}, {-6, -6}, {6, 6}}, {{v, {-6, -3}}, {-6, -6}, {6, 6}},  
         ControlType -> Locator] ]

Пример 5.3.4. Пример пространственной кривой Безье, построенной по 8 опорным точкам. Все точки управляемы покоординатно с помощью ползунков.

In[10]: =
       DynamicModule [ {pts , р0 , pi, р2 , рЗ , р4 , р5 , р6 , р7 , t, g} , 
         Manipulate [pts := {р0 , pi, p2 , p3 , p4 , p5 , р6 , p7} ; 
             Show[ { ParametricElot3D[g[pts] [t] , {t, 0, 1}, PlotStyle -> {Green, Thick}, 
                 Axes -> True, AxesEdge -> Automatic, AxesOrigin -> {0, 0, 0}, 
                 PlotRange -> {{-15, IS}, {-15, 15}, {-15, 15}}], 
              Graphics3D [ {Red, PointSize [Large] , Point [pts [[#]] ] /@Range[8], 
                  Text[ToString[# - 1] , pts [ [*f] ] + {0.5, 0.5, 0.5}]  /@ Range [8] }] } ] , 
           {{p0, {0.0, 0.0, 0.0}}, {-1, -1, -1}, {1, 1, 1}}, 
           {{p1, {2.0, 1.0, 1.0}}, {1, 0, 0}, {3, 2, 2}}, 
           {{p2, {1.8, 5.2, 2.4}}, {0.8, 4.2, 1.4}, {2.8, 6.2, 3.4}}, 
           {{p3, {-2.3, 6.5, 3.0}}, {-3.3, 4.5, 2.0}, {-1.3, 7.5, 4.0}}, 
           {{p4, {-4.0, 4.0, 5.0}}, {-5.0, 3.0, 4.0}, {-3.0, 5.0, 6.0}}, 
           {{p5, {-7.0, 1.0, 6.0}}, {-8.0, 0.0, 5.0}, {-6.0, 2.0, 7.0}}, 
           {{рб, {4.8, 2.2, 8.4}}, {3.8, 1.2, 7.4}, {5.8, 3.2, 9.4}}, 
           {{p7, {9.3, 6.5, 12.0}}, {8.3, 5.5, 11.0}, {10.3, 7.5, 13.0}}, 
           ControlType -> Slider 
        ], Initialization: -> (g[pts_] : = BezierFunction [pts] ) ]

Для построения кривых Безье существуют разные алгоритмы. Рассмотрим один из наиболее популярных алгоритмов.

Алгоритм де Кастелье

Алгоритм де Кастелье позволяет построить кривую Безье с помощью рекуррентного соотношения

$$B_i^n=tB_{i-1}^{n-1}+(1-t)B_i^{n-1}.$$

Пусть даны опорные точки $$p_0, \dots , p_n.$$ Обозначим через $$r_i(t; k) $$ кривую Безье, построенную по точкам $$p_i, \dots , p_{i+k} (i + k \le n) .$$

Тогда искомая кривая Безье - это кривая $$r_0(t;n) .$$ Имеем $$r_i(t;0) = p_i, 0 \le i \le n.$$ Выведем рекуррентные формулы для определения $$r_0(t; n) = r(t) $$ исходя из начальных данных: $$r_i(t; 0) \equiv p_i (0 \le i \le n) .$$ Имеем

$$r(t)=p_0(1-t)^n+p_nt^n+\sum_{i=1}^{n-1}(tB_{i-1}^{n-1}+(1-t)B_i^{n-1}p_i=\\ =(1-t)(p_0(1-t)^{n-1}+\sum_{i=0}^{n-1}B_i^{n-1}p_i)+t(p_nt^{n-1}+\sum_{i=0}^{n-1}B_{i-1}^{n-1}p_i)=\\ =(1-t) \sum_{i=0}^{n-1}B_i^{n-1}p_i+t \sum_{i=0}^{n-1}B_i^{n-1}p_{i+1}=(1-t)r_0(t; n-1)+tr_1(t; n-1).$$

Отсюда получаем следующую рекуррентную формулу де Кастелье для вычисления кривой Безье:

$$r_i(t;k)=(1-t)r_i(t; k-1)+tr_{i+1}(t;k-1),$$

где $$i + k \le n.$$

По формуле (5.22), начиная от значений $$r_i(t;0) = p_i,$$ получаем последовательно значения $$r_i (t;k), 1 \le k \le n.$$ На последнем шаге вычисляем $$r_0(t;n) = r(t) .$$ Это и есть алгоритм де Кастелье.

Пример 5.3.5. Алгоритм де Кастелье:

In[11]: =
        DynamicModule [ {s, r, b, t,   p} , 
             b[t_, i_, r_]   := (1 - t) b[t, i, r- 1] + tb[t, i + 1, r - 1] ; 
             b[t_, 0, 1] =  (1-t) p[0] + tp[l];   b[t_, 1, 1] =  (1-t) p[l] +tp[2]; 
             b[t_, 2,  1]  =  (1-t) p[2] + tp[3] ; p[0]  = {0.0,  0.5}; 
             p[l]  = {0.5, 2.5} ; p[2]  = {3.0, 2.0}; p[3]  = {4.0,  0.0}; 
             Manipulate[ 
                Module[{pit}, 
                   plt=ParametricPlot[b[t, 0, 3] ,   {t, 0, s} ,
                              PlotStyle -> {Blue, Thickness [0.01] }] [ [1] ] ; 
                   Which [r == 1, 
                        Graphics [{{Lighter [Gray] , Line[{p[0] ,p[l],p[2],p[3]}]}, 
                          {PointSize[0.015] , Point[{p[0] , p [1] , p [2] , p[3]}]}, 
                          {Purple, PointSize[0.015],
                            Point [{b[s,  0,  l],b[s, 1, l],b[s, 2, 1]}]}, 
                          {Purple, Thickness[0.005], Line[{p[0], b[s, 0, 1]}], 
                            Line[{p[l] , b[s, 1, 1]}] , Line [{p [2], b[s, 2, 1]}]}}, 
                            PlotRange -> {{-0.1, 4.1},   {-0.1, 2.6}}, ImageSize -> {300, 275}] , 
                       r == 2, Graphics [{{Lighter [Gray] , Line[{p[0] ,p[l],p[2],p[3]}]}, 
                          {PointSize [0.015] , Point [{p[0] ,p[l],p[2],p[3]}]}, 
                          {Lighter[Gray], 
                            Line[{{b[s, 0, l],b[s, 1, 1]},   {b[s, 1, l],b[s, 2, 1]}}]}, 
                          {Purple, PointSize[0.015], 
                            Point [{b[s, 0, 1] , b[s, 1, 1], b[s, 2, 1], b[s, 0, 2] , b[s, 1, 2]}]}, 
                          {Purple,  Thickness[0.005], 
                            Line[{{p[0], Ь[з,  0, 1]},   {p[l], b[з, 1,  1]},   {p[2], b[s, 2,  1]}, 
                             {b[s,  0,  1], b[s, 0, 2]},   {b[s,  1,  1], b[s,  1, 2]}}]}}, 
                     PlotRange -> {{-0.1,  4.1},   {-0.1, 2.6}}, ImageSize -> {300, 275}] , 
                  r == 3, Graphics [ {pit,   {Lighter [Gray] , Line[{p[0] , p [1] , p[2] , p[3] } ]} , 
                     {PointSize [0.015] , Point [{p[0] ,p[l],p[2],p[3]}]}, 
                     {Lighter[Gray], 
                       Line[{{b[s,  0, 1] , b[s, 1,  1]},   {b[s, 1,  l],b[s, 2, 1]}, 
                          {b[s, 0, 2], b[s,  1, 2]}}]}, 
                     {Purple, PointSize[0.015], 
                       Point[{b[s, 0, 1] , b[s, 1, 1] , b[s, 2, 1] , b[s, 0, 2] , b[s, 1, 2] ,
                            b[s, 0, 3]}]},   {Purple, Thickness[0.005], 
                      Line[{{p[0], b[s, 0, 1]},   {p[l], b[s, 1, 1]},   {p[2], b[s, 2, 1]}, 
                          {b[s,  0,  1], b[s,  0, 2]},   {b[s,  1,  1], b[s,  1, 2]}, 
                          {b[s,  0, 2], b[s, 0,  3]}}]}}, 
               PlotRange -> {{-0.1,  4.1},   {-0.1, 2.6}},  ImageSize -> {300, 275}]] ] , 
          {{r,  3,  "Шаг рекурсии."},   {1, 2,  3}}, 
         {{s, 0.5, "t   (параметр на кривой)"}, 0.01, 1}] ]

Пример 5.3.6. Следующий пример строит последовательность кривых Безье в форме "цветка", соседние кривые имеют только одну общую крайнюю точку. Все опорные точки расположены на двух окружностях: первая имеет радиус 1, а радиус второй управляется ползунком.

In[12]: =
             Manipulate[ 
               DynamicModule[{ptlists, img, pts2}, 
                 pts2 = 
                  Table [ If [ EvenQ[i] , 1. , r] *
                       {Sin [2 Pii *iMul/ iDiv] , Cos [2 Pi i * iMul / iDiv] } ,
                   {i,0, iLoops / 2 * LCM [iDiv, iDeg, iMul] } ] ;
                ptlists = Partition [ pts2 , iDeg + 1, iDeg]; 
                img = Graphics[{
                   Table[{If[bColor, Hue[2 *k/Length[ptlists]] , Black] , 
                       BezierCurve [ptlists [ [k] ] , SplineDegree -* iDeg] } , 
                    {k, Length[ptlists]}] }, ImageSize -" {375, 375}, Axes 4 True]] , 
          {{bColor, True, "Расцветка кривых Безье"}, {True, False}}, 
          { {iDiv, 23,
             "Количество равностоящих точек окружностей, в которых (возможно 
                 не во всех) будут располагаться все опорные точки"}, 2, 
           40, 1, ImageSize -> Tiny, Appearance -> "Labeled"} , 
         {{iDeg, 9, "Степень кривых Безье"} , 1, 11, 1, ImageSize -> Tiny, 
              Appearance -> "Labeled" } , 
          {{iMul, 7, "Шаг, с помощью которого выбираются опорные точки"},
              1, 10, 1, ImageSize->" Tiny, Appearance -> "Labeled"}, 
          {{r, .5, "Радиус второй окружности"}, 0, 1.5, .01, ImageSize -> Tiny, 
              Appearance -> "Labeled" } , 
         {{iLoops, 1, "Добавление опорных точек"}, 1, 4,1, ImageSize -> Tiny, 
              Appearance -> "Labeled" } , 
            Control Placement -> Left]

Операторная форма кривой Безье

Пусть даны опорные точки $$p_0, \dots ,p_n$$ и параметр $$t \in [0,1] .$$ Определим оператор правого сдвига E, действующий на кривых вида $$r(t) =\sum_{i=0}^{n-1}f_i(t)p_i,$$ где $$f_i (t), i = 0, \dots , n - 1,$$ - гладкие функции от $$t \in [0,1] ,$$ по формуле

$$Er(t)=\sum_{i=0}^{n-1}f_i(t)p_{i+1}$$

Аналогично, на кривых вида $$r(t) = \sum_{i=0}^{n-1}g_i(t)p_i, t \in [0,1] ,$$ определен оператор левого сдвига L:

$$Lr(t)=\sum_{i=1}^{n}g_i(t)p_{i-1}$$

Тогда на множестве кривых, для которых определен оператор $$L \circ E$$ или $$E \circ L,$$ имеем $$L \circ E = Id, E \circ L = Id.$$ Другими словами, $$L = E^{-1}, E = L^{-1}.$$ Мы будем пользоваться также степенями операторов $$L$$ и $$E: L^k = L \circ \dots \circ L, E^k = E \circ \dots \circ E$$ (где $$L$$ и $$E$$ повторяются $$k$$ раз), применяя их к тем кривым, для которых эти степени определены.

Кривую $$r(t) = p_i$$ мы будем обозначать просто $$p_i.$$

Лемма 5.1. Кривая Безье, построенная по точкам $$p_0, \dots ,p_n,$$ может быть задана в виде

$$r(t) = ( (1-t)L + t )^n p_n,$$

а также в виде

$$r(t) = ( (1-t)+tE )^n p_0$$

В операторных обозначениях алгоритм де Кастелье имеет вид

$$r_i(t; k) = ( (1 -t)+ tE ) r_i(t; k-1), (i + k \le n),$$

поскольку $$Er_i(t; k - 1) = r_{i+1}(t; k - 1) .$$

В частности

$$r_0(t;n)= ( (1-(t)+tE ) r_0(t;n-1),\\ r_0(t; n - 1) = (1 -t) + tE ) r_0(t; n - 2),\\ ............................................................\\ r_0(t;1) = p_0,$$

что опять дает нам формулу (5.24):

$$r(t) = r_0(t; n) = ( (1 -t)+ tE )^n p_0.$$

Лемма 5.2. Для любого значения $$t \in [0,1] $$ соответствующая точка $$r(t) $$ на кривой Безье $$r(t) = r0(t; n) $$ делит отрезок, соединяющий точки $$r_0(t; n-1) $$ и $$r_1(t;n-1) ,$$ в отношении $$t : (1-t) .$$ Сама кривая $$r(t) = r_0(t;n) $$ касается этого отрезка в точке $$r(t) .$$

Рассмотрим две операции на кривых Безье.

Деление кривой Безье на две кривые Безье того же порядка в отношении t* : (1 - *)

Пусть даны опорные точки $$p_0, \dots , p_n$$ и кривая Безье $$r = r(t) ,$$ построенная по ним. Требуется для некоторого заданного $$t^* \in [0,1] $$ построить опорные точки $$p_0^a, \dots , p_n^a$$ и $$p_0^b, \dots , p_п^b$$ так, чтобы кривые Безье для этих новых точек совпали бы с двумя дугами исходной кривой Безье, на которые она разбивается точкой $$r(t*).$$

Данная задача имеет важное значение в теории кривых Безье. Необходимость деления кривых Безье возникает в тех случаях, когда требуется изменить какую-то часть кривой, оставляя некоторую другую ее часть неизменной. Поскольку смещение любой опорной точки меняет всю кривую Безье в целом, в таких случаях приходится разбивать кривую Безье на несколько частей. Это дает возможность менять только одну из получившихся частей исходной кривой Безье (ставших уже самостоятельными кривыми Безье), не трогая остальные части исходной кривой.

Для решения поставленной задачи введем новый параметр $$\tau =\frac {t}{t^*}.$$ Имеем $$t = \tau t^*.$$ Запишем формулу (5.24) в виде

$$r(t) = r_0(t; n) = (1 - \tau t^* + \tau t^*E)^np_0 = ( (1-\tau ) + (1-t^*+t^*E) \tau )^n p_0$$

Введем новый оператор сдвига $$E_a = (1-t^*+t^*E) $$ и перепишем (5.25) в виде (5.24), положив $$p_0^a = p_0, p_i^a = E_i^a p_0^a$$:

$$r_0(t; n) = r_0^a(t; n) = ( (1 -\tau) + \tau E_a )^n p_0^a,$$

где $$\tau \in [0,1], t = \tau t^*.$$

Таким образом, если мы построим кривую Безье по $$(n + 1) $$ точкам $$p_0^a, \dots ,p_n^a,$$ определенным по формуле $$p_i^a = (1 - t^* + t^*E)^ip_0,$$ эта новая кривая Безье совпадает с дугой исходной кривой Безье $$r(t), 0 \le t \le t^*.$$

Аналогично, с помощью нового параметра $$\xi = \frac{t-t^*}{1-t^*}$$ и нового оператора левого сдвига $$L_b = (1-t^*)L+t^* (L = E^{-1})$$ определим $$(n+1) $$ новую опорную точку $$p_0^b, \dots ,p_n^b$$ по формуле $$p_{n_i}^b = L_b^i p_n^b, i = 0, \dots ,n, p_n^b = p_n.$$ Тогда дуга исходной кривой Безье на участке $$[t, 1] $$ может быть задана формулой (5.23), которую можно записать в виде

$$r(t) = ( (1-t)L + t )^n p_n =\\ = (( 1 - (1 - t^*) \xi - t*) L+(1 -t^*) \xi+t^*)^n p_n =\\ = ( ((1-t^*)L + t^*) (1-\xi ) + \xi )^n p_n= ( (1-\xi ) L_b +\xi )^n p_n.$$

Следовательно, эта дуга будет совпадать с кривой Безье, построенной по точкам $$p_0^b,\dots , p_n^b.$$

Процедура построения точек $$p_0^a, \dots ,p_n^a$$ и $$p_0^b, \dots ,p_n^b$$ может быть описана более наглядно следующим образом.

Исходный многоугольник с вершинами в опорных точках назовем многоугольником нулевого уровня $$p_0^{(0)} = p_0, \dots ,p_n^{(0)} = p_n.$$

Строим многоугольники уровней $$1,2, \dots ,n$$ так, что вершины многоугольника следующего уровня делят стороны многоугольника предыдущего уровня в отношении $$t^* : (1 - t^*)$$:

$$p_i^{(k)} = (1 - t^*)p_i^{(k-1)} + t^*p_{i+1}^{(k-1)}, i = 0, \dots ,n - k.$$

Многоугольник $$k$$ -го уровня будет иметь $$(n-k+1) $$ вершин: $$p_0^{(k)}, \dots ,p_{n-k}^{(k)} .$$ При этом, согласно построению, сдвигая все исходные вершины вперед, мы сдвигаем вперед и все вершины многоугольника любого уровня, т. е. $$Ep_i^{(k)} = p_{i+1}^{(k)} .$$ Тогда

$$p_i^a =p_0^{(i)}.$$

В самом деле, имеем по построению,

$$p_0^{(0)} =p_0,\\ p_0^{(1)} = (1 - t^*)p_0 + t^*p_1 = (1 - t^* + t^*E)p_0 = E_ap_0^{(0)},\\ …………………………………………….\\ p_0^{(k)} = (1 - t^*)p_0^{(k-1)} + t^*p_1^{(k-1)} = (1 - t^* + t^*E)p_0^{(k-1)} = E_ap_0^{(k-1)},$$

так как $$Ep_0^{(k-1)} = p_1^{(k-1)} .$$ Следовательно, $$p_0^{(k)} = E_a^kp_0^{(0)} = E_a^kp_0 = p_k^a ,$$ где $$E_a = (1 - t^*) + t*E .$$

Аналогично строятся точки $$p_0^b, \dots ,p_n^b ,$$ только надо брать последовательность последних, а не первых вершин многоугольников $$k $$ -го уровня.

Условия сохранения гладкости сопряжения при делении кривой Безье

Деление кривой Безье, описанное выше, применяется в частности тогда, когда надо управлять формой этой кривой, зафиксировав некоторую точку на ней (в этой точке мы и делим кривую). Однако при последующем управлении двумя частями кривой гладкость сопряжения в точке деления может падать. Если мы хотим, чтобы эта гладкость была не менее $$C^k ,$$ то мы можем добиться этого, оставив на месте точки

$$p_n^a, \dots, p_{n-k}^a,\\ p_o^b, \dots, p_k^b.$$

Данное утверждение следует из того, что первые $$k $$ производных на конце кривой Безье определяются соответствующими крайними $$(k+1) $$ опорными точками. Поскольку исходная кривая была бесконечно гладкой, то, оставляя на месте точки (5.27), мы обеспечим сопряжение производных до порядка $$k $$ включительно в точке деления.

Увеличение числа опорных точек без изменения формы кривой Безье

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

Для кривой Безье с опорными точками $$p_0, \dots , p_n $$ имеем

$$r(t) = r_0(t;n) = (1-t + Et)^np_0 $$

или

$$r(t)=\sum_{i=0}^{n}C_n^i(1-t)^{n-i}t^ip_i$$

Произведем следующие преобразования параметрического представления кривой Безье с целью увеличения числа опорных точек на единицу без изменения формы исходной кривой:

$$r(t)=(1-t+t)\sum_{i=0}^{n}C_n^i(1-t)^{n-i}t^ip_i=\\ =\sum_{i=0}^{n}C_n^i(1-t)^{n-i+1}t^ip_i+\sum_{i=0}^{n}C_n^i(1-t)^{n-i}t^{i+1}p_i=\\ =\sum_{i=0}^{n+1}C_{n+1}^i(1-t)^{n+1-i}t^i \left ( \frac{C_n^i}{C_{n+1}^i}(1 - \delta_{(n+1)i})p_i+ \frac{C_n^{i-1}}{C_{n+1}^i}(1- \delta_{0i})p_{i-1} \right ),$$

где $$\delta_{ij} $$ - символ Кронекера-Капелли:

$$\delta_{ij}= \begin{cases} 1 \mbox{ при} i=j,\\ 0 \mbox{ при} i \ne j \end{cases} $$

Положим теперь

$$q_i=\frac{n+1-i}{n+1}p_i(1-\delta_{(n+1)i})+\frac{i}{n+1}p_{i-1}(1-\delta_{0i}), i=0, \dotsm n+1$$

Отметим, что $$p_{n+1} $$ и $$p - 1 $$ входят в последнее равенство с нулевыми коэффициентами, поэтому доопределять эти (не определенные пока) точки нет необходимости.

Рациональные кривые Безье

Обычные кривые Безье получаются с помощью итерирования операции деления отрезка в отношение $$t : (1 -t) .$$ Так, например: $$p_t = (1 -t)p_0 +tp_i = (1 -t+tE)p_0 = r_0(t; 1) $$ - кривая Безье первого порядка (здесь отрезок $$p_{0pl} $$ разделен в отношении $$t : (1 - t)$$ ). Аналогично, $$r_0(t; n) = (1 - t + tE)^np_0 $$ -кривая Безье порядка $$n .$$ Здесь уже операция деления отрезка в отношении $$t:(1-t) $$ проитерированаn раз.

Рациональные кривые Безье получаются таким же образом с помощью операции рационального деления отрезка с весами $$\omega_0 $$ и $$\omega_х .$$

Определение 5.4.1. Рациональной кривой Безье порядка $$n ,$$ построенной по заданным опорным точкам $$p_0, \dots , p_n $$ с заданными весами $$\omega_0, \dots ,\omega_n ( \omega_i > 0) ,$$ называется кривая

$$r_0(t;n)=\frac{(1-t+tE)^n(\omega_0p_0)}{(1-t+tE)^n\omega_0}=\frac{(1-t-tE)^n(\omega_0p_0)}{\omega_0(t;n)},$$

где $$t \in [0,1], \omega_0(t; n) = (1 -t + tE)^n\omega_0 .$$ Здесь, как и выше, оператор $$E $$ действует одновременно и на последовательности опорных точек $$p_0, p_1 \dots , p_n ,$$ и на последовательности весов $$\omega_0,\omega_1, \dots , \omega_n .$$

Согласно биному Ньютона, рациональная кривая Безье $$r_0(t;n) $$ может быть представлена в виде

$$r_0(t;n)=\sum_{i=0}^nC_n^it^i(1-t)^{n-i}\frac{\omega_i}{\omega_0(t;n)}p_i$$

В случае равных весов $$\omega_0 = \dots = \omega_n $$ знаменатель $$\omega_0(t; n) $$ превращается в постоянную, а рациональная кривая Безье - в обычную кривую Безье.

Замечание 5.4.1.

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

Так же как и для обычной кривой Безье, в своих крайних опорных точках рациональная кривая Безье касается звеньев своей опорной ломаной (ломаной, построенной по опорным точкам).

Так же как и для обычной кривой Безье, $$i $$ -тая производная на конце рациональной кривой Безье определяется $$(i + 1) $$ крайней опорной точкой, а также весами этих точек.

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

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

При смещении опорной точки $$p_i $$ на вектор $$\xi_i $$ все точки рациональной кривой Безье смещаются в направлении вектора $$\xi_i $$ (на разные величины). Поскольку обычная кривая Безье является частным случаем рациональной кривой Безье (в случае равенства весов), то же верно и для обычной кривой Безье.

При увеличении веса $$\omega_i $$ опорной точки pi все точки рациональной кривой Безье испытывают смещение в направлении опорной точки pi (на разные величины). Это смещение может оказаться нулевым (например, концы кривой не смещаются при изменении весов), но вектор смещения всегда направлен от исходного положения точки кривой к опорной точке $$p_i ,$$ если вес этой опорной точки увеличивается, и в противоположную сторону, если он уменьшается.

Рациональные кривые Безье в пакете Mathematica. Mathematica не имеет встроенной поддержки рациональных кривых Безье. В частности, функция BezierFunction[pts] не работает, если в качестве ее аргумента взять последовательность вещественных чисел (весов), а не координат точек пространства размерности 2 или больше. Поэтому при вычислении знаменателя формулы радиус-вектора рациональной кривой Безье с помощью Mathematica следует дополнить список весов еще одной координатой до двумерных векторов, а затем использовать лишь первую компоненту получившейся двумерной вектор-функции.

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

In[13]: =
           DynamicModule [ {pts0 , w0, wl, w2, w3; w4 , w5 , w6, w7 , ww, n} , 
              Pts0 = {{-10.0, 0.0}, {-3.0, 7.0}, {1.1, -7.2}, {4.3, -8.5},
                   {-16.0, 12.0}, {-2.1, 1.7}, {-8.2, 1.4}, {12.2, 5.4}}; 
                n = Length [pts0] ; 
                Manipulate[ 
                ww = {w0 , wl, w2 , w3 , w4 , w5 , w6 , w7 } ; 
                 Show[ 
                   {Graphics [{Red, Text [ToString [# - 1] , pts [ [#] ] +{1.2, 0.2}] /@Range[n], 
                      Text[ 
                          "Рациональная кривая Безье с весами $$\omega^i$$ в сравнении с 
                               обычной кривой 
                              Безье", {-7.2, 13}]}], 
                        ParametricPlot[{fBez [pts, ww] [t] / fwBez [ww] [t] [ [1] ] , f0Bez [pts] [t] } , 
                              {t, 0, 1}] 
                          } , Axes -> True, PlotRange -> {{-30, 15} , {-10, 15}}] , 
                        {{pts, pts0}, Locator},
                        {{w0, 1}, 1, 100}, {{wl, 1}, 1, 100}, {{w2, 1}, 1, 100},
                        {{w3, 1}, 1, 100}, {{w4, 1}, 1, 100}, {{w5, 1}, 1, 100}, {{w6, 1}, 1, 100}, 
                        {{w7, 1} , 1, 100} , 
                        Initialization : -> (
                               f0Bez [JGT;LS_] : = BezierFunction [pts] ;
                               fwBez [ww_]   := BezierFunction[Table [{w[ [i] ] , 1}, {i, 1, n} ] ] ;
                               fBez [pts_f ww_]   : = BezierFunction [pts ww] ; )]]

Деление рациональной кривой Безье

Для заданного $$t^* $$ из интервала $$(0,1) $$ требуется построить два новых набора опорных точек и весов $$\{p_0^a, \dots , p_n^a \}, \{\omega_0^a, \dots \omega_n^a\} $$ и $$\{p_0^b, \dots , p_п^b\},\{\omega_0^b, \dots ,\omega_n^b\} $$ так, чтобы кривые Безье, построенные по этих новым наборам, совпали бы (как геометрическое место точек) соответственно с двумя дугами исходной кривой Безье, на которые она разбивается точкой $$r(t^*) .$$

Данная задача имеет важное значение в теории кривых Безье. Необходимость деления кривых Безье возникает в тех случаях, когда требуется изменить какую-то часть кривой, оставляя некоторую другую ее часть неизменной. Поскольку смещение любой опорной точки кривой Безье меняет всю ее в целом, возникает необходимость разбивать кривую Безье на несколько частей. Это дает возможность менять только одну из получившихся частей исходной кривой Безье (ставших уже самостоятельными кривыми Безье), не трогая остальные части исходной кривой.

Деление рациональной кривой Безье сводится к делению обычной кривой Безье, которая порождает данную рациональную кривую.

Рассмотрим рациональную кривую Безье, построенную по опорным точкам $$p_0, \dots , p_n $$ с весами $$\omega_0 > 0, \dots , \omega_n > 0 .$$ Обозначим ее через $$r_0(t; n), t \in [0,1] .$$ Как мы знаем, эта кривая порождается обычной кривой Безье, построенной по вспомогательным опорным точкам $$\{\hat p_0, \dots ,\hat p_n\} = \{\omega_0p_0, \dots , \omega_np_n\} ,$$ которые зависят от выбора начала координат. При этом начало координат необходимо выбирать так, чтобы оно не совпадало ни с одной из опорных точек $$p_0, \dots , p_n .$$ Обозначим построенную таким образом вспомогательную обычную кривую Безье через $$\hat r_0(t;n), t \in [0,1] .$$ Разделим вспомогательную кривую $$\hat r_0(t;n) $$ в точке $$t = t^* $$ в отношении $$t^* : (1 - t^*) ,$$ как это делалось выше. Таким образом, возникает три вспомогательные (обычные) кривые Безье:

$$\hat r_0(t;n) $$ - это кривая Безье, построенная по опорным точкам $$\hat p_0, \dots ,\hat p_n , t \in [0,1],$$

$$\hat r_0^a(s;n) $$ - это кривая Безье, построенная по опорным точкам $$\hat p_0^a, \dots ,\hat p_n^a, s \in [0,1],$$

$$\hat r_0^b(q;n) $$ - это кривая Безье, построенная по опорным точкам p $$\hat p_0^b, \dots , \hat p_n^b, q \in [0,1].$$ Здесь

$$\hat p_i^a=\hat r_0(t^*;i)=(1-t^*+Et^*)^i \hat p_0,\\ \hat p_{n-i}^b=\hat r_{n-i}(t^*;i)=((1-t^*)L+t^*)^i \hat p_n, 0 \le i \le n$$

Кривые Безье $$\hat r_0^a 0(s; n) $$ и $$\hat r_0^b(q; n) $$ порождают две рациональные кривые Безье $$r_0^a(s;n) $$ и $$r_0^b(q;n) $$ по следующему правилу:

$$r_0^a(s;n) $$ - рациональная кривая Безье, построенная по опорным точкам $$p_0^a, \dots ,p_n^a $$ с весами $$\omega_0^a, \dots ,\omega_n^a ,$$

$$r_0^b(q;n) $$ - рациональная кривая Безье, построенная по опорным точкам $$p_0^b, \dots ,p_n^b $$ с весами $$\omega_0^b, \dots ,\omega_b^b ,$$ где

$$p_i^a=r_0(t^*;i)=\frac{(1-t^*+Et*)^i\omega_0p_0}{(1-t^*+Et^*)^i\omega_0},\\ p_{n-i}^b=r_{n-i}(t^*;i)=\frac{((1-t^*)L+t^*)^i\omega_np_n}{((1-t^*)L+t^*)^i \omega_n},$$ $$\omega_i^a=(1-t^*+Et^*)^i \omega_0,\\ \omega_{n-i}^b=((1-t^*)L+t^*)^I \omega_n$$

Из формул (5.28), (5.29) и (5.30) следует, что

$$\hat p_i^a=\omega_i^ap_i^a,\\ \hat p_i^b=\omega_i^bp_i^b, i=0, \dots, n$$

Равенства (5.31) означают, что кривые Безье $$\hat r_0^a $$ и $$\hat r_0^b $$ соответствуют рациональным кривым Безье $$r_0^a $$ и $$r_0^b $$ с весами $$\{\omega_i^a\} $$ и $$\{\omega_i^b\} $$ соответственно.

Лемма 5.3. Кривые $$r_0^a(s; n) $$ и $$r_0^b(q; n) $$ делят кривую $$r_0(t; n) $$ в точке $$t = t* $$ в отношении $$t* :(1-t) .$$

Доказательство см. в [1].

Увеличение числа опорных точек рациональной кривой Безье

Операция увеличения числа опорных точек рациональной кривой Безье сводится к уже известной нам операции увеличения числа опорных точек обычной кривой Безье в пространстве $$R_+^{k+1} $$ и применению отображения $$\varphi : R_+^{k+1} \to R^k, \varphi (x, \omega) =\frac{x}{\omega} .$$ Более подробно, алгоритм состоит в следующем.

По данной рациональной кривой Безье $$r_0(t; n) ,$$ определяемой опорными точками $$p_0, \dots , p_n $$ и весами $$\omega_0, \dots , \omega_n ,$$ строим обычную кривую Безье $$r_0(t; n) $$ в пространстве $$R_+^{k+1} ,$$ определяемую опорными точками $$\{ (\omega_0p_0, \omega_0), \dots , (\omega_n p_n, \omega_n) \} , \tilde p_i = (\omega_i p_i,\omegai) \in R_+^{k+1} .$$ Тогда $$\varphi: \tilde r_0(t;n) \to r_0(t;n) .$$

Увеличиваем число опорных точек обычной кривой Безье $$r_0(t;n) $$ на единицу способом, описанным в разделе 3.6. А именно, строим новые опорные точки $$\tilde q_i \in R^{k+1}, i = 0, \dots , n + 1 ,$$ по формуле

$$\tilde q_i = \frac{n+1-i}{n+1} \tilde p_ i (1 -\delta_{(n+1)i}) +\frac{i}{n+1} \tilde p_{i-1}(1 -\delta_{0 i}), i = 0, \dots ,n + 1.$$

Точки $$\tilde p_{-1}$$ и $$\tilde p_{n+1}$$ можно определить произвольно, т. к. они входят в правую часть последней формулы с нулевыми коэффициентами.

Применяем к новым опорным точкам $$\tilde q_i$$ преобразование $$\varphi : q_i = \varphi (\tilde q_i) .$$ Поскольку кривая Безье $$\tilde r_0(t;n + 1) ,$$ определяемая опорными точками $$\{\tilde q_0, \dots ,\tilde q_{n+1}\},$$ - это та же самая кривая, что и $$\tilde r_0(t;n) ,$$ то $$\varphi (\tilde r_0(t;n+1) ) =r_0(t;n) .$$ Но одновременно по определению рациональной кривой Безье имеем: $$\varphi (\tilde r_0(t;n+1) ) = r_0(t;n+ 1) ,$$ где $$r_0(t;n+ 1) $$ - рациональная кривая Безье, определяемая опорными точками и весами $$\{ (q_0,\omega_0?), \dots , (q_{n+1}, \omega_{n+1}') \},$$ где $$\omega_i? =\tilde q_i^{ k+1} -(k + 1) $$ -е компоненты точек $$\tilde q_i$$ в $$R^{k+1}$$

Окончательные формулы имеют вид

$$\begin{cases} q_i=\frac{i \omega_{i-1}p_{i-1}(n+1-i)\omega_ip_i}{i \omega_{i-1}+(n+1-i) \omega_i}, \omega_i'=\frac{i \omega_{i-1}+(n_1-i)\omega_i}{n+1},\\ \end{cases},\\ i=0, \dots, n+1$$

где $$\{ (p_i , \omega_i) \}_{i =0}^n$$ - исходные опорные точки и веса данной рациональной кривой Безье i а $$\{(q_i, \omega_i')\}_{i=0}^{n+1}$$ - новые опорные точки и веса той же кривой, но - с увеличенным на единицу количеством опорных точек (и, соответственно, с большим на единицу количеством весов).

Страницы: