.3. B-КРИВЫЕ

Основные недостатки кривых Безье следующие:

1) смещение любой опорной точки  влияет на всю кривую, за исключением ее крайних точек;

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

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

Поставим задачу следующим образом. Пусть даны опорные точки (Subscript[p, 1]),...,(Subscript[p, n]) с весами (Subscript[ω, 1]),...,(Subscript[ω, n]). Будем искать  кривую r(t) в виде:    r(t)=(((∑_(i=1)^n)(Subscript[N, i])(t)(Subscript[ω, i])(Subscript[p, i]))/((∑_(i=1)^n)(Subscript[N, i])(t)(Subscript[ω, i]))), (Subscript[t, min]) ≤ t≤ (Subscript[t, max]),
где (Subscript[N, i])(t) — набор кусочно-полиномиальных функций таких, что

1) (Subscript[N, i])(t)=0 при t∈[(Subscript[t, min]),(Subscript[t, max])]\((Subscript[a, i]),(Subscript[b, i])) при некоторых (Subscript[t, min])<(Subscript[a, i]) и (Subscript[b, i])<(Subscript[t, max]).

2) для любой гладкой функции f(t) существует линейная комбинация ∑(Subscript[c, i])(Subscript[N, i])(t), интерполирующая функцию f(t) в заданных узлах (Subscript[t, i]): f((Subscript[t, i]))=(∑_(j=1)^n)(Subscript[c, j])(Subscript[N, j])((Subscript[t, i])), i=1,...,n;

3)   (∑_(i=1)^n)(Subscript[N, i])(t)=1 для каждого t∈[(Subscript[t, min]),(Subscript[t, max])].

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

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

.3.1. РАЗДЕЛЕННЫЕ РАЗНОСТИ

Пусть имеется гладкая функция f(t) и фиксированы некоторые значени (Subscript[t, 1]),...,(Subscript[t, m+1]) переменной t (не обязательно различные и необязательно упорядоченные по возрастанию или убыванию t), которые мы будем называть УЗЛАМИ. Мы также будем рассматривать переменное значение t, которое будем считать самым первым узлом. Поскольку он предшествует узлу (Subscript[t, 1]), припишем ему индекс 0 (хотя в большинстве случаев будем обозначать его просто t без индекса): (Subscript[t, 0]) = t,(Subscript[t, 1]),...,(Subscript[t, m+1]).  Обозначим (Subscript[f, i])=f((Subscript[t, i])), (Subscript[f, t])=f(t)=f((Subscript[t, 0]))=(Subscript[f, 0]).

Определение. РАЗДЕЛЕННОЙ РАЗНОСТЬЮ функции f(t) порядка k≤m-i+1 в узле (Subscript[t, i]) называется величина (R^((k)))(Subscript[f, i])=(((R^((k-1)))(Subscript[f, i+1])-(R^((k-1)))(Subscript[f, i]))/((Subscript[t, i+k])-(Subscript[t, i]))); 
(R^((0)))(Subscript[f, i])=(Subscript[f, i]).

В случае совпадения узлов (Subscript[t, i])=(Subscript[t, i+1])=...=(Subscript[t, i+k]) в правой части формулы для разделенной разности появляется отношение вида (0/0). Эта неопределенность может быть раскрыта путем предельного перехода по возмущенным узлам, что дает производную функции f(t) порядка k в точке (Subscript[t, i])=(Subscript[t, i+1])=...=(Subscript[t, i+k]):

(R^((k)))(Subscript[f, i])=((1)/(k!))(((d^k))/(dt^k))f((Subscript[t, i])).

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

ПРОГРАММА 1. Разделенная разность пятого порядка функции Cos(t) как функция первого узла.

ЗАМЕЧАНИЕ. Программа не вычисляет производные при совпадении узлов.

DynamicModule[{n=5,tt,F,h,L,t,t2,t3,t4,t5,f,p,S,z},
f[t_]:=Cos[t];
Manipulate[
tt=Join[{t},ss];
p[i_,s_]:=Apply[Times,(s-tt[[#]])&/@Drop[Range[n],{i}]];
S=Flatten[Table[Table[ss[[j]]-ss[[i]],{i,1,j-1}],{j,2,Length[ss]}]];
z=Apply[Times,S];
If[z≠0,
L=Table[1/p[i,tt[[i]]],{i,1,n}];
F=Table[f[tt[[i]]],{i,1,n}];
h=F.L/.t->0;
Show[
Plot[{f[t],F.L/h},{t,-20,20},PlotRange->{{-25,25},{-10,10}}],
Graphics[{Red,PointSize[Large],Point[Table[{ss[[i]],0},{i,1,Length[tt]-1}]],Text[ToString[#],{ss[[#]]+1.2,0}]&/@Range[4],Axes->True}]
],
Show[Graphics[Text["Узлы совпадают"],PlotRange->{{-25,25},{-20,20}}]]
],
{{ss,{-7.2,0.7,5.2,12.0},"Узлы"},
With[{r=Range[1,Length[tt]-1]},
Grid[{Spacer[3],Subscript[Style["t",Italic],#-1],Slider[Dynamic@ss[[#1]],{-20,20,.01},Appearance->"UpArrow",ImageSize->{200,12}],Dynamic@ss[[#1]]}&/@r]]&,ControlPlacement->Bottom},TrackedSymbols ->{ss},SaveDefinitions->True
]
]

УПРАЖНЕНИЕ. Обратите внимание, что разделенная разность быстро убывает при удалении вправо и влево от крайних узлов. Каков порядок этого убывания?

УПРАЖНЕНИЕ. Написать программу, которая бы учитывала совпадение узлов.

.3.2. УСЕЧЕННАЯ СТЕПЕННАЯ ФУНКЦИЯ И ЕЕ РАЗДЕЛЕННЫЕ РАЗНОСТИ

Определение. Пусть 0≤m  — целое число. Рассмотрим функцию (Subscript[σ, m])(z,t)=((SubPlus[(z-t)])^m), где индексом ``+'' мы обозначаем положительную часть функции: (SubPlus[f])(t)=max{f(t),0}. При фиксированном t функция (Subscript[σ, m])(z,t), рассматриваемая как функция z, называется УСЕЧЕННОЙ СТЕПЕННОЙ ФУНКЦИЕЙ.

ПРОГРАММА 2. Разделенная разность пятого порядка усеченной кубической функции. Разделенная разность вычисляется от  (Subscript[σ, m])(z,t)=((SubPlus[(z-t)])^m), рассматриваемой как функция z при фиксированном t. Программа демонстрирует равенство нулю разделенной разности при всех t, лежащих левее и правее всех узлов.

f[t_,z_]:=If[t<z,(z-t)^3,0];
Plot[f[-6,z],{z,-20,20}]

DynamicModule[{h,z,t,st,M,G,i,x,y},
h[x_,y_]:=If[x<y,(y-x)^3,0];
Manipulate[
p[i_,s_]:=Apply[Times,(s-st[[#]])&/@Drop[Range[Length[st]],{i}]];
z=Apply[Times,(st[[#]]-st[[#+1]])&/@Range[Length[st]-1]];
If[z≠0,
M=Table[1/p[i,st[[i]]],{i,1,Length[st]}];
G=Table[h[t,st[[i]]],{i,1,Length[st]}];
Show[
If[
OrderedQ[st],
Plot[G.M,{t,-20,20}],
Do[If[st[[i-1]]>st[[i]],st[[i-1]]=st[[i]]],{i,2,Length[st]}];
Plot[G.M,{t,-20,20}]],
Graphics[{Red,PointSize[Large],Point[Table[{st[[i]],0},{i,1,Length[st]}]],Text[ToString[#-1],{st[[#]]+1.2,0}]&/@Range[5],Axes->True}]],
Show[Graphics[Text["Узлы совпадают"],PlotRange->{{-25,25},{-20,20}}]]
],
{{st,{-13.2,-6.0,0.7,5.2,12.0},"Узлы"},
With[{r=Range[1,Length[st]]},
Grid[{Spacer[3],Subscript[Style["t",Italic],#-1],Slider[Dynamic@st[[#1]],{-15,13,.01},Appearance->"UpArrow",ImageSize->{200,12}],Dynamic@st[[#1]]}&/@r]]&,ControlPlacement->Bottom},TrackedSymbols ->{st},SaveDefinitions->True
]
]

УПРАЖНЕНИЕ. Написать программу, которая бы учитывала совпадение узлов.

.3.3. B-СПЛАЙНЫ

На основе разделенных разностей усеченной степенной функции определим функциональные коэффициенты (Subscript[N, i])(t) .

Пусть  (Subscript[t, 0])≤(Subscript[t, 1])≤...≤(Subscript[t, m]), тогда (Subscript[σ, m-1])[(Subscript[t, 0]),...,(Subscript[t, m])](t) – разделенная разность (Subscript[σ, m-1])(z,t), рассматриваемой как функция z при фиксированном t.

Определения. НОРМИРОВАННЫМ B-СПЛАЙНОМ m-го ПОРЯДКА (Subscript[N, i,m])(t), привязанным справа к первому узлу (Subscript[t, i]) неубывающей последовательности узлов (Subscript[t, i])≤(Subscript[t, i+1])≤...≤(Subscript[t, i+m]), называется функция (Subscript[N, i,m])(t)=((Subscript[t, i+m])-(Subscript[t, i]))(Subscript[σ, m-1])[(Subscript[t, i]),...,(Subscript[t, i+m])](t).

НЕНОРМИРОВАННЫМ B-СПЛАЙНОМ m-го ПОРЯДКА (Subscript[M, i,m])(t), привязанным справа к первому узлу (Subscript[t, i]) неубывающей последовательности узлов (Subscript[t, i])≤(Subscript[t, i+1])≤...≤(Subscript[t, i+m]), называется функция (Subscript[M, i,m])(t)=(Subscript[σ, m-1])[(Subscript[t, i]),...,(Subscript[t, i+m])](t).

НОРМИРОВАННЫМ B-СПЛАЙНОМ m-го ПОРЯДКА (Subscript[N, m,i])(t), привязанным слева к последнему узлу (Subscript[t, i]) неубывающей последовательности узлов (Subscript[t, i-m])≤(Subscript[t, i-m+1])≤...≤(Subscript[t, i]), называется функция (Subscript[N, i,m])(t)=((Subscript[t, i])-(Subscript[t, i-m]))(Subscript[σ, m-1])[(Subscript[t, i-m]),...,(Subscript[t, i])](t).

НЕНОРМИРОВАННЫМ B-СПЛАЙНОМ m-го ПОРЯДКА (Subscript[N, m,i])(t), привязанным слева к последнему узлу (Subscript[t, i]) неубывающей последовательности узлов (Subscript[t, i-m])≤(Subscript[t, i-m+1])≤...≤(Subscript[t, i]), называется функция (Subscript[M, i,m])(t)=(Subscript[σ, m-1])[(Subscript[t, i-m]),...,(Subscript[t, i])](t)

ЗАМЕЧАНИЕ. Иногда порядок B-сплайна определяют не числом узлов минус единица, а степенью усеченной функции (как многочлена t). При таком подходе введенные выше B-сплайны будут иметь (по определению) порядок не m, а (m-1). В этом случае для них применяются другие обозначения: (Subscript[N, i-]^(m-1))(t)=(Subscript[N, i,m])(t),(Subscript[M, i-]^(m-1))(t)=(Subscript[M, i,m])(t),(Subscript[N, -i]^(m-1))(t)=(Subscript[N, m,i])(t),(Subscript[M, -i]^(m-1))(t)=(Subscript[M, m,i])(t).  И те и другие обозначения присутствуют в литературе по B-кривым. Мы будем пользоваться первыми. В Mathematica используются вторые.

ТЕОРЕМА КОКСА-ДЕ БУРА.  Для ненормированных B-сплайнов при любых t,1≤m, α, β=α+m справедливы следующие формулы:

(Subscript[M, α,m])(t)=((((Subscript[t, β])-t)(Subscript[M, α+1,m-1])(t)+(t-(Subscript[t, α]))(Subscript[M, α,m-1])(t))/((Subscript[t, β])-(Subscript[t, α])));

(Subscript[M, m,β])(t)=((((Subscript[t, β])-t)(Subscript[M, m-1,β])(t)+(t-(Subscript[t, α]))(Subscript[M, m-1,β-1])(t))/((Subscript[t, β])-(Subscript[t, α]))).
При этом, ненормированные B-сплайны первого порядка по определению полагаются равными разделенным разностям функции (Subscript[σ, 0])(z)=((SubPlus[(z-t)])^0):

(Subscript[M, i,1])(t)=((1)/((Subscript[t, i+1])-(Subscript[t, i]))), если t∈[(Subscript[t, i]),(Subscript[t, i+1])), и 0, иначе;

(Subscript[M, 1,i])(t)=((1)/((Subscript[t, i])-(Subscript[t, i-1]))), если t∈[(Subscript[t, i-1]),(Subscript[t, i])), и 0, иначе.

Для нормированных B-сплайнов при любых t,1≤m, α, β=α+m справедливы следующие формулы:

(Subscript[N, α,m])(t)=(((Subscript[t, β])-t)/((Subscript[t, β])-(Subscript[t, α+1])))(Subscript[N, α+1,m-1])(t)+((t-(Subscript[t, α]))/((Subscript[t, β-1])-(Subscript[t, α])))(Subscript[N, α,m-1])(t);

(Subscript[N, m,β])(t)=(((Subscript[t, β])-t)/((Subscript[t, β])-(Subscript[t, α+1])))(Subscript[N, m-1,β])(t)+((t-(Subscript[t, α]))/((Subscript[t, β-1])-(Subscript[t, α])))(Subscript[N, m-1,β-1])(t);

(Subscript[N, i,1])(t)=1, если t∈[(Subscript[t, i]),(Subscript[t, i+1])), и 0, иначе;

(Subscript[N, 1,i])(t)=1, если t∈[(Subscript[t, i-1]),(Subscript[t, i])), и 0, иначе.

Теорема Кокса-Де Бура позволяет вычислять B-сплайны рекуррентно.

ЗАМЕЧАНИЕ. Для сокращения объема вычислений следует сначала  вычислить ненормированные B-сплайны нужного порядка по формулам Кокса-Де Бура, а затем нормировать их.

УПРАЖНЕНИЕ. С помощью встроенной функции Mathematica BSplineBasis[{m, {(Subscript[t, 0]), (Subscript[t, 1]), ...,(Subscript[t, n])}},i, t], вычисляющей тот из n-m нормированныx B-сплайнов (Subscript[N, i,m+1])(t) порядка m+1 (или m при втором способе определения порядка, см. выше), который привязан справа к i-тому (0≤i≤n-(m+1)) узлу из расширенного множества неубывающих узло (Subscript[t, 0])≤(Subscript[t, 1])≤...≤(Subscript[t, n]) , наглядно изучить ненормированные и нормированные B-сплайны, привязанные к узлам справа. Обратите внимание, что в Mathematica узлы в расширенном множестве нумеруются начиная с нуля. Соответственно, из (n-m) В-сплайнов в базисе нумеруются в Mathematica с нулевого до n-(m+1)-го.

ЗАМЕЧАНИЕ.  При работе с нижеследующей программой обратите внимание, что при совпадении первых m+1 узлов, значение нулевого (самого первого) нормированного В-сплайна из базиса в нулевом (самом первом) узле становится равно 1, тогда как все остальные В-сплайны базиса обращаются в этом узле в ноль. Это значит, что когда такой базис из нормированных В-сплайнов используется для вычисления радиус-вектора В-кривой (аналогично базису Бернштейна для кривых Безье), получившаяся кривая будет начинаться в своей крайней левой опорной точке. Аналогично, при совпадении последних m+1 узлов расширенного множества, последний сплайн базиса в последнем узле примет значение 1, а остальные будут равны в нем нулю, и, следовательно, правый конец соответствующей кривой будет находится в ее последней опорной точке. См.
Утверждение ниже, в формулировке которого порядок В-сплайна m соответствует значению m+1 в обозначениях Mathematica (как уже говорилось, Mathematica использует второй способ определения порядка В-сплайна, а значит, уменьшает его на единицу по сравнению с тем способом, который используем мы).

УПРАЖНЕНИЕ. Написать программу на Mathematica, вычисляющую B-сплайны (Subscript[N, m,i])(t), привязанные к узлам слева.

ПРОГРАММА 3. Следующие две программы строят графики В-сплайнов, привязанных к узлам справа. В данном примере n=9, то есть расширенное множество узлов состоит из 10 точек, начиная с нулевой, а m=3, то есть порядок В-сплайна при принятом нами определении равен 4.

DynamicModule[{m=3,i,tt,knots},
For[i=1,i≤m+1,i++,tt[i]=0]; (*m+1 начальных узло равны 0*)
tt[5]=4;
tt[6]=4;
tt[7]=6;
tt[8]=9; (*задали m+1 конечных узлов*)
Manipulate[
knots=Join[Table[tt[i],{i,1,m+1}],ss,Table[tt[j],{j,m+2,2m+2}]];
Show[
If[OrderedQ[ss],
Plot[BSplineBasis[{m,knots},0,t],{t,knots[[1]],1.1},
PlotRange->{{0,1.1},{0,1.2}},PlotLabel->Style[(Subscript[N, 0,m+1])[t],14]],
Do[If[ss[[i-1]]>ss[[i]],ss[[i-1]]=ss[[i]]],{i,2,Length[ss]}];
Plot[BSplineBasis[{m,knots},0,t],{t,knots[[1]],knots[[Length[knots]]]+1},
PlotRange->{{0,1.1},{0,1.2}}],PlotLabel->Style[(Subscript[N, 0,m+1])[t],14]
]],
{{ss,{0.5,0.7},"Узлы"},
With[{r=Range[1,Length[ss]]},
Grid[{Spacer[3],Subscript[Style["t",Italic],#+m],Slider[Dynamic@ss[[#1]],{0,1,.01},Appearance->"UpArrow",ImageSize->{200,12}],Dynamic@ss[[#1]]}&/@r]]&,ControlPlacement->Bottom},TrackedSymbols ->{ss}, SaveDefinitions->True
]
]

DynamicModule[{m=3,n=9,i,j,tt,knots,a},
For[i=1,i≤m+1,i++,tt[i]=0]; (*m+1 начальных узло равны 0*)
For[j=1,j≤m+1,j++,tt[n-j]=7]; (*задали m+1 конечных узлов*)
Manipulate[
knots=Join[Table[tt[i],{i,1,m+1}],ss,Table[tt[j],{j,m+2,2m+2}]];
a=(knots[[1]]+knots[[n]])/2;
(*f[t_]:=Table[BSplineBasis[{m,knots},l,t],{l,0,n-m-1}];*)
Show[
If[OrderedQ[ss],
Plot[Evaluate[Table[BSplineBasis[{m,knots},l,t],{l,0,n-m-1}]],{t,knots[[1]],8},PlotRange->{{knots[[1]],8},{0,1}},Filling->If[fill,{i+1->{Axis,{Opacity[.15],Green}}},None],PlotLabel->Style[(Subscript[N, i,m+1])[t],14]],
Do[If[ss[[i-1]]>ss[[i]],ss[[i-1]]=ss[[i]]],{i,2,Length[ss]}];
Plot[Evaluate[Table[BSplineBasis[{m,knots},l,t],{l,0,n-m-1}]],{t,knots[[1]],8},PlotRange->{{knots[[1]],8},{0,1}},Filling->If[fill,{i+1->{Axis,{Opacity[.15],Green}}},None],PlotLabel->Style[(Subscript[N, i,m+1])[t],14]]],
Graphics[Arrow[
{{2.5,1.1},{i+0.1,BSplineBasis[{m,knots},i,i+0.1]}}]]
],
{{i,0,"Номер полинома i"},0,n-m-1,1},
Delimiter,
{{ss,{0.5,1.7},"Узлы"},
With[{r=Range[1,Length[ss]]},
Grid[{Spacer[3],Subscript[Style["t",Italic],#+m],Slider[Dynamic@ss[[#1]],{0,2,.01},Appearance->"UpArrow",ImageSize->{200,12}],Dynamic@ss[[#1]]}&/@r]]&,ControlPlacement->Bottom},
{{fill,True,"Заполнение"},{True,False},ControlPlacement->Bottom},SaveDefinitions->True
]
]

.3.4. B-КРИВЫЕ

Определение. B-КРИВОЙ порядка m, построенной по опорным точкам (Subscript[p, 1]),...,(Subscript[p, n]), m≤n, с весами (Subscript[ω, 1]),...,(Subscript[ω, n]), называется кривая
r(t)=(((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i])(Subscript[p, i]))/((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i]))),
определяемая опорными точками ((Subscript[{(Subscript[p, i])}, i=1])^n), весами ((Subscript[{(Subscript[ω, i])}, i=1])^n) и неубывающей последовательностью узлов (Subscript[t, 1]),...,(Subscript[t, n+m]).

ЗАМЕЧАНИЯ.

1) Управляющими параметрами B-кривой являются:
    –  порядок кривой m,
    – ее опорные точки ((Subscript[{(Subscript[p, i])}, i=1])^n),
    – веса опорных точек ((Subscript[{(Subscript[ω, i])}, i=1])^n),
    – неубывающая последовательность вещественных чисел (Subscript[t, 1]),...,(Subscript[t, n+m]), называемая РАСШИРЕННЫМ МНОЖЕСТВОМ УЗЛОВ B-КРИВОЙ.

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

2) В отличие от кривых Безье, B-кривые, вообще говоря, не проходят через свои крайние опорные точки. Однако такого прохождения нетрудно добиться с помощью следующего утверждения.

УТВЕРЖДЕНИE. Если первые m узлов расширенного множества совпадают, то B-кривая начинается в своей первой опорной точке (Subscript[p, 1]). Аналогично, если последние m узлов  расширенного множества совпадают, то B-кривая заканчивается в своей последней опорной точке (Subscript[p, n]).

Mathematica имеет встроенные функции BSplineFunction и BSplineCurve, которые позволяют строить В-кривые по заданному множеству опорных точек пространства (R^n), n>1, автоматическим или полуавтоматическим подбором параметров. Расширенное множество узлов по умолчанию подбирается автоматически, однако его  можно задать самому с помощью параметра SplineKnots. Например, SplineKnots -> Automatic означает автоматический подбор расширенного множества узлов, а SplineKnots -> {0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 2} задает расширенное множество узлов непосредственно.

Так же дело обстоит и с весами опорных точек. Веса можно задавать с помощью SplineWeights. Например, SplineWeights -> Automatic означает автоматический выбор весов, а SplineWeights -> {1, 1, 10, 1, 1, 1} задает конкретные веса для шести опорных точек.

Порядок сплайнов можно задавать с помощью SplineDegree, но при этом должно быть SplineKnots  -> Automatic .

Вычисление порядка происходит так. Число сплайнов в базисе, равное n-m (см. выше), должно равняться числу опорных точек. Поэтому порядок сплайнов может быть вычислен как разность числа узлов и числа опорных точек. Однако для Mathematica эта разность должна быть подобрана так, чтобы она равнялась 4, иначе В-кривая может быть построена в Mathematica неправильно.

Чтобы воспользоваться возможностью полноценного управления В-кривой, в том числе и весов, нужно запрограммировать В-кривую самостоятельно с помощью встроенной функции BSplineBasis.

УПРАЖНЕНИЕ. С помощью встроенных функций Mathematica SplineFunction и BSplineCurve построить В-кривую по заданным опорным точкам с автоматическим выбором узлов на плоскоси и в трехмерном пространстве.

ПРОГРАММА 4. В-кривая на плоскости с автоматическим выбором узлов, построенная двумя способами: с помощью функций BSplineFunction и BSplineCurve соответственно. На втором графике для сравнения наприсована также простая кривая Безье, построенная по тем же опорным точкам.

DynamicModule[{pts0,f,g,n,pts1,pts2,t},
n=Length[pts0];
Row[{
Manipulate[f=BSplineFunction[pts1];
 Show[{
Graphics[{Green,Line[pts1],Black,Text["В-кривая, построенная через BSplineFunction",{7,10.0}],Text[ToString[#-1],pts1[[#]]+{0.7,0}]&/@Range[n]},PlotRange->{{-4,20},{-4,10}},Axes->True],ParametricPlot[f[t],{t,0,1}]
}],{{pts1,pts0},Locator}],
Manipulate[g=BezierFunction[pts2];
Show[{
Graphics[{Red,BSplineCurve[pts2,SplineKnots->Automatic],Black,Text["В-кривая, построенная через BSplineCurve (красная)",{11,10.0}],
Text["Криая Безье (голубая)",{10,9.0}],Text[ToString[#-1],pts2[[#]]+{0.7,0}]&/@Range[n]},PlotRange->{{-4,20},{-4,10}},Axes->True],
ParametricPlot[g[t],{t,0,1},PlotStyle->{Blue}]
}],{{pts2,pts0},Locator}
]
},"    "],Initialization:>(
pts0={{-1,0.0},{0,0},{0.5,1.5},{2.0,3.0},{4,0},{6,3.0},{8.0,2.0},{9.0,0},{1.0,1.0},{0.0,3.0},{2.0,4.},{4.0,1.5},{6.0,6.0},{11,3},{12.0,0.0},{14.5,1.5},{16,6},{17.0,3.0},{15.0,0.0},{12.0,2.0},{12.0,4.0},{14.0,7.0}}
)]

ПРОГРАММА 5. В-кривая в пространстве с автоматическим выбором узлов в сравнении с кривой Безье, построенной по тем же опорным точкам. Обратите внимание насколько более послушно ведет себя В-кривая при управлении опорными точками, чем кривая Безье. Кривая Безье по-другому  — более глобально — усредняет влияние опорных точек на кривую, что делает ее локальное управление за счет смещения отдельных опорных точек гораздо более вязким, чем в случае В-кривой).

DynamicModule[{pts0,f,pts,n},
Manipulate[
f=BezierFunction[pts];
Show[{
Graphics3D[{Cyan,Thick,BSplineCurve[pts,SplineKnots->Automatic],Red,PointSize[Large],Point[pts],Text[ToString[#-1],pts[[#]]+{0.4,0.4,0.4}]&/@Range[n],Black,Text["X",{18,-0.3,-0.3}],Black,Text["Y",{-0.3,18,-0.3}],Black,Text["Z",{-0.3,-0.3,18}]},Axes->True,AxesEdge->Automatic,AxesOrigin->{0,0,0}],ParametricPlot3D[f[t],{t,0,1},PlotStyle->{Thick}]},Boxed->False,PlotRange->{{-18,18},{-18,18},{-18,18}}],{
{pts,pts0,"Узлы"},
With[{r=Range[1,n]},
Grid[{Spacer[3],Subscript[Style["(t^1)",Italic],#-1],
Slider[Dynamic@pts[[#1]][[1]],{-15,15,.01},Appearance->"Left",ImageSize->{200,12}],Dynamic@pts[[#1]][[1]]}&/@r]]&,ControlPlacement->Left},
Delimiter,
{
{pts,pts0,"Узлы"},
With[{r=Range[1,n]},
Grid[{Spacer[3],Subscript[Style["(t^2) и (t^3)",Italic],#-1],
Slider2D[Dynamic@{pts[[#1]][[2]],pts[[#1]][[3]]},{{-15,-15},{15,15}},Appearance->"Left",ImageSize->{50,30}],Dynamic@{pts[[#1]][[2]],pts[[#1]][[3]]}}&/@r]]&,ControlPlacement->Left}],Initialization:>(
pts0={{0.0,0.0,0.0},{2.0,1.0,1.0},{1.8,5.2,2.4},{-2.3,6.5,3.0},{-4.0,4.0,5.0},{-7.0,1.0,6.0},{4.8,2.2,8.4},{9.3,6.5,12.0},{12.3,4.5,13.0},{11.0,3.0,15.0},{0.0,-1.0,13.0},{4.8,-2.2,14.4},{1.3,3.5,12.0}};
n=Length[pts0]
)]

ПРОГРАММА 6. В-кривая на плоскости с автоматическим выбором порядка, но с управляемыми мышью опорными точками и узлами (узлы должны быть неубывающими).

DynamicModule[{pts0,tt0,n,m,f,tt,i,pts},
pts0={{-9.0,0.0},{-4.0,5.0},{1.1,7.2},{4.3,8.5},{6.0,12.0},{12.1,14.7},{7.2,9.4},{11.5,5.9},{7.0,7.0},{14.1,4.7},{11.2,-1.4},{-5.5,10.9}};
tt0=Table[i-1,{i,1,16}];
n=Length[pts0];
m=Length[tt0];
Manipulate[
f=BSplineCurve[pts,SplineKnots->tt];
Show[{If[OrderedQ[tt],
Graphics[f,Axes->True,PlotRange->{{-15,20},{-10,20}}],
Do[If[tt[[i-1]]>tt[[i]],tt[[i-1]]=tt[[i]]],{i,2,m}];
Graphics[f,Axes->True]
],
Graphics[{Red,Point[pts],Text[ToString[#-1],pts[[#]]+{0.4,0.4}]&/@Range[n]},PlotRange->{{-15,20},{-10,20}}]}],
{{pts,pts0},Locator},
{{tt,tt0,"Узлы"},
With[{r=Range[1,m]},
Grid[{Spacer[3],Subscript[Style["t",Italic],#-1],Slider[Dynamic@tt[[#1]],{0,15,.01},Appearance->"UpArrow",ImageSize->{200,12}],Dynamic@tt[[#1]]}&/@r]]&,ControlPlacement->Left}],Initialization:>(
pts0={{-9.0,0.0},{-4.0,5.0},{1.1,7.2},{4.3,8.5},{6.0,12.0},{12.1,14.7},{7.2,9.4},{11.5,5.9},{7.0,7.0},{14.1,4.7},{11.2,-1.4},{-5.5,10.9}};
tt0=Table[i-1,{i,1,16}];
n=Length[pts0];
m=Length[tt0]
)]

ПРОГРАММА 7. В-кривая на плоскости с автоматическим выбором порядка и узлов, но с управляемыми мышью опорными точками и весами. На том же графике изображена также простая кривая Безье, управляемая теми же опорными точками.

DynamicModule[{pts0,n,w0,f,g,d,curve,w},
pts0={{-9.0,0.0},{-4.0,5.0},{1.1,7.2},{4.3,8.5},{6.0,12.0},{12.1,14.7},{10.2,9.4},{7.5,7.9},{2.0,4.0},{1.1,1.2},{-17.2,18.4},{-16.5,-15.9}};
n=Length[pts0];
w0=Table[1,{i,1,n}];
Manipulate[
f=BezierFunction[pts];
g=BSplineCurve[pts,SplineWeights->w,SplineDegree->d];Show[Graphics[{Red,Point[pts],Text[ToString[#-1],pts[[#]]+{0.4,0.4}]&/@Range[n]},Axes->True],If[curve===1,Graphics[g],ParametricPlot[f[t],{t,0,1},PlotStyle->{Red}]]],
{{d,2,"Степень B-кривой"},Dynamic@RadioButtonBar[Dynamic[d],Range[0,7]]&},
{{curve,1,"Кривая"},{1->"B-кривая",2->"Кривая Безье"}},
{{pts,pts0},Locator},
Delimiter,
{{w,w0,"Узлы"},
With[{r=Range[1,n]},
Grid[{Spacer[3],Subscript[Style["w",Italic],#-1],
Slider[Dynamic@w[[#1]],{0,4,0.01},Appearance->"Left",ImageSize->{200,12}],Dynamic@w[[#1]]}&/@r]]&,ControlPlacement->Left}],Initialization:>(
pts0={{-9.0,0.0},{-4.0,5.0},{1.1,7.2},{4.3,8.5},{6.0,12.0},{12.1,14.7},{10.2,9.4},{7.5,7.9},{2.0,4.0},{1.1,1.2},{-17.2,18.4},{-16.5,-15.9}};
n=Length[pts0];
w0=Table[1,{i,1,n}]
)
]

УПРАЖНЕНИЕ. Самостоятельно запрограммировать В-кривую с помощью встроенной функции BSplineBasis, так, чтобы можно было воспользоваться возможностью полноценного управления В-кривой. Обеспечить возможность произвольного допустимого порядка В-сплайнов.

.3.5. АЛГОРИТМ ВЫЧИСЛЕНИЯ РАДИУС-ВЕКТОРА  B-КРИВОЙ

Пусть даны опорные точки (Subscript[p, 1]),...,(Subscript[p, n]) и их веса (Subscript[ω, 1]),...,(Subscript[ω, n]), а также расширенное множество неубывающих узлов {(Subscript[t, 1]),...,(Subscript[t, n+m]):(Subscript[t, 1])≤...≤(Subscript[t, n+m])}. На основании вышесказанного опишем алгоритм вычисления радиус-вектора B-кривой

r(t)=(((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i])(Subscript[p, i]))/((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i]))).

1) Фиксируем t. Положим (Subscript[t, 0]) равным минус бесконечность, а (Subscript[t, n+m+1]) – плюс бесконечность. Тогда существует единственное (Subscript[i, 0]) такое, что t∈[(Subscript[t, (Subscript[i, 0])]),(Subscript[t, (Subscript[i, 0])+1])),  где 0≤(Subscript[i, 0])≤n+m;

2) Для определенного выше индекса (Subscript[i, 0]), если 1≤(Subscript[i, 0])≤n+m, вычисляем единственное ненулевое значение ненормированного B-сплайна первого уровня (m=1): (Subscript[M, (Subscript[i, 0]),1])(t)=((1)/((Subscript[t, (Subscript[i, 0])+1])-(Subscript[t, (Subscript[i, 0])]))).
Для индекса (Subscript[i, 0]), находящегося вне множества {1,...,n}, значение (Subscript[M, (Subscript[i, 0,1])]) не вычисляется.

3) С помощью соотношений Кокса-Де Бура вычисляем все отличные от нуля в точке t ненормированные сплайны m-го порядка (Subscript[M, j,m])(t) при  1≤j≤n:
(Subscript[M, j,m])(t)=((((Subscript[t, j+m])-t)(Subscript[M, j+1,m-1])(t)+(t-(Subscript[t, j]))(Subscript[M, j,m-1])(t))/((Subscript[t, j+m])-(Subscript[t, j]))).

4) Вычисляем нормированные сплайны (Subscript[N, j,m])(t) для каждого 1≤j≤n по формуле: (Subscript[N, j,m])(t)=((Subscript[t, j+m])-(Subscript[t, j]))(Subscript[M, j,m])(t).

5) Окончательно вычисляем r(t) по формуле  
r(t)=(((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i])(Subscript[p, i]))/((∑_(i=1)^n)(Subscript[N, i,m])(t)(Subscript[ω, i]))).

2.1. ПОВЕРХНОСТИ, ОПРЕДЕЛЯЕМЫЕ МАТРИЦАМИ ОПОРНЫХ ТОЧЕК И ВЕСОВ

2.1.1. ПОВЕРХНОСТИ БЕЗЬЕ

Определение. Пусть дано (m+1)(n+1) точек в пространстве (R^3), образующих прямоугольную матрицу – сетку опорных точек поверхности Безье (Subscript[P, m n])=(Subscript[((Subscript[p, i j])), i=0,m;j=0,n]).  Соответствующая этой сетке ПОВЕРХНОСТЬ БЕЗЬЕ порядка m x n определяется  двупараметрическим выражением, где оба параметра u, v меняются  на отрезке [0,1]:

S(u,v;m,n)=F(u,m)P (F^T)(v,n),

где F(u,m)=((Subscript[f, 0])(u,m),...,(Subscript[f, m])(u,m)), F(v,n)=((Subscript[f, 0])(v,n),...,(Subscript[f, n])(v,n)), (Subscript[f, i])(ξ,k)=((((Subscript[C, k])^i)(1-ξ))^(k-i))(ξ^i).

РАЦИОНАЛЬНАЯ ПОВЕРХНОСТЬ БЕЗЬЕ определяется аналогично рациональной кривой Безье по сетке опорных точек (Subscript[P, m n]) и матрице их весов (Subscript[W, m n])=(Subscript[((Subscript[w, ij])), i=0,m;j=0,n]): элементами которой являются вещественные положительные числа (Subscript[w, i j]) – веса опорных точек (Subscript[p, i j]). Радиус-вектор рациональной поверхности Безье определяется формулой:

(Subscript[S, рац])(u,v;m,n)=((F(u,m)Q (F^T)(v,n))/(F(u,m)W (F^T)(v,n))),

где матрица Q=((Subscript[q, i j])) состоит из элементов  (Subscript[q, i j])=(Subscript[w, i j])(Subscript[p, i j]).

УПРАЖНЕНИЕ. Наглядно изучить влияние опорных точек и их весов на рациональную поверхность Безье, меняя исходные данные pts (опорные точки) и  w (их веса) в нижеследующей программе. В матрице весов стоят двумерные вектора, но используется только их первая компонента. Вторая компонента фиксирована. Это связано с неспособность Mathematica применять функцию BezierFunction к матрицам из скаляров.

ПРОГРАММА 1. Рациональная поверхность Безье с возможностью непосредственного управления весами с помощью движков.

DynamicModule[{pts,a,w0,pw,w,g,f,w6,w7,w10,w11,i,j,out,ins,u,v,n,pts0},
pts={{{0,0,0},{0,1,0},{0,2,0},{0,3,0}},{{1,0,0},{1,1,1},{1,2,1},{1,3,0}},{{2,0,0},{2,1,1},{2,2,1},{2,3,0}},{{3,0,0},{3,1,0},{3,2,0},{3,3,0}}};
pts0=Flatten[pts,1];
n=Length[pts0];
Row[{
Manipulate[
w0={{{2,a},{3,a},{4,a},{5,a}},{{3,a},{w5,a},{w6,a},{6,a}},{{4,a},{w9,a},{w10,a},{7,a}},{{9,a},{12,a},{14,a},{32,a}}};
w=Table[w0[[i,j,1]],{i,1,4},{j,1,4}];
pw=Table[w0[[i,j,1]] pts[[i,j]],{i,1,4},{j,1,4}];
g=BezierFunction[w0];
f=BezierFunction[pw];
Show[{
Graphics3D[{PointSize[Large],Red,Map[Point,pts],Green,Text[ToString[ #-1],pts0[[#]]+{0.01,0.04,0.04}]&/@Range[n]}],
Graphics3D[{Gray,Dashed,Line[pts],Line[Transpose[pts]]}],
ParametricPlot3D[f[u,v]/g[u,v][[1]],{u,0,1},{v,0,1},PlotStyle->FaceForm[out,ins]]
}],
Column[{
Control[{{ins,Green,"Внутренний цвет"},Green}],
Control[{{out,Red,"Внешний цвет"},Red}]
},Right],
{{w5,10},1,100},{{w6,10},1,100},{{w9,10},1,100},{{w10,10},1,100}
]
MatrixForm[pts]
},"   "],
Initialization:>(
pts={{{0,0,0},{0,1,0},{0,2,0},{0,3,0}},{{1,0,0},{1,1,1},{1,2,1},{1,3,0}},{{2,0,0},{2,1,1},{2,2,1},{2,3,0}},{{3,0,0},{3,1,0},{3,2,0},{3,3,0}}};
pts0=Flatten[pts,1];
a=1;
n=Length[pts0]
)
]

2.1.2. ГЕОМЕТРИЧЕСКИЙ СМЫСЛ ПОВЕРХНОСТИ БЕЗЬЕ

Геометрически, поверхность Безье получается следующим образом:

1. На первом шаге строим (n+1) кривую Безье по опорным точкам, расположенным в столбцах матрицы (Subscript[P, m n]). Будем считать параметр на всех этих кривых согласованным и обозначим его через u.

2. Далее, для каждого фиксированного значения параметра u∈[0,1] строим кривую Безье, имеющую опорные точки, расположенные на построенных выше (n+1) кривых так, что все они соответствуют выбранному значению параметра u. Эти кривые будут заметать некую гладкую поверхность. Она и будет искомой поверхностью Безье.

Можно показать, что если начать строить кривые не со столбцов матрицы (Subscript[P, m n]), а с ее строк, то поверхность в итоге получится та же самая.

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

2.1.3. ДЕЛЕНИЕ ПОВЕРХНОСТИ БЕЗЬЕ

Деление поверхностей Безье целиком основано на делении кривых Безье. Как мы уже знаем, каждая из кривых Безье (обычных или рациональных), построенных по столбцам сетки опорных точек (Subscript[P, m n]), может быть разделена в параметрической точке u=(SuperStar[u]) ∈ [0,1] без изменения ее формы. При таком делении получаются две поверхности Безье (составляющие вкупе исходную поверхность) и соответствующие им две матрицы опорных точек Безье, каждая из которых имеет тот же порядок, что и исходная матрица: (Subscript[P, m n]^a) и (Subscript[P, m n]^b). Затем, каждую из получившихся двух поверхностей Безье надо поделить в параметрической точке v=(SuperStar[v]) (с помощью деления кривых Безье, построенных по строкам соответствующих матриц). В итоге получится четыре поверхности Безье и соответствующие им четыре матрицы опорных точек: (Subscript[P, m n]^(a c)),  (Subscript[P, m n]^(a d)), (Subscript[P, m n]^(b c)) и (Subscript[P, m n]^(b d)). Эти четыре поверхности Безье будут лежать на исходной поверхности Безье, делить ее на четыре части и в совокупности составлять исходную поверхность.

2.1.4. ДОБАВЛЕНИЕ НОВЫХ СТРОК И СТОЛБЦОВ ОПОРНЫХ ТОЧЕК (ИЗМЕЛЬЧЕНИЕ СЕТКИ БЕЗЬЕ) БЕЗ ИЗМЕНЕНИЯ ПОВЕРХНОСТИ БЕЗЬЕ

Добавление новых опорных точек в сетку Безье без изменения самой поверхности происходит по формулам, полностью аналогичным рассмотренным выше формулам добавления новых опорных точек в кривую Безье без ее изменения.  Если дана матрица опорных точек поверхности Безье (Subscript[P, m n])=(Subscript[((Subscript[p, i j])), i=0,m;j=0,n]), то новая, увеличенная на одну строку и на один столбец матрица опорных точек (Subscript[Q, m+1 n+1]) для той же самой поверхности Безье строится по формулам:

(Subscript[q, i j])=(1/(m+1))(1/(n+1))((m+1-i)(1-(Subscript[δ, (m+1)i]))+i(1-(Subscript[δ, 0i]))(E^(-1))((n+1-j)(1-(Subscript[δ, (n+1)j]))+j(1-(Subscript[δ, 0 j]))(F^(-1)))(Subscript[p, i j]),

где E и F – операторы сдвига вперед по первому и второму индексу соответственно: E (Subscript[p, i j])=(Subscript[p, (i+1)j]) и F (Subscript[p, i j])=(Subscript[p, i(j+1)]).
При этом угловые точки новой сетки совпадут с углами старой: (Subscript[q, 00])=(Subscript[p, 00]), (Subscript[q, (m+1)0])=(Subscript[p, m 0]), (Subscript[q, 0(n+1)])=(Subscript[p, 0n]), (Subscript[q, (m+1)(n+1)])=(Subscript[p, m n]).

2.1.5.  ВЫЧИСЛЕНИЕ ДИФФЕРЕНЦИАЛЬНО-ГЕОМЕТРИЧЕСКИХ ХАРАКТЕРИСТИК
     ПОВЕРХНОСТИ БЕЗЬЕ В ПРОИЗВОЛЬНОЙ ТОЧКЕ

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

Таким образом, достаточны следующие формулы.

Касательные вектора к поверхности Безье в угловой точке (Subscript[p, 00]) могут быть вычислены по формулам: (Subscript[T, u])(0)=m((Subscript[p, 10])-(Subscript[p, 00])) и (Subscript[T, v])(0)=n((Subscript[p, 01])-(Subscript[p, 00])).

Нормаль к поверхности Безье в угловой точке (Subscript[p, 00]) дается формулой: n=(((Subscript[T, u])(0) X (Subscript[T, v])(0))/(|(Subscript[T, u])(0) X (Subscript[T, v])(0)|)).

Вторые производные радиус-вектора поверхности Безье по параметрическим координатам в точке (Subscript[p, 00]) равны:

((∂^2)/(∂(u^2)))(Subscript[|, u=0,v=0])S(u,v;m,n)=m(m-1)((Subscript[p, 20])-2(Subscript[p, 10])+(Subscript[p, 00])),

((∂^2)/(∂(v^2)))(Subscript[|, u=0,v=0])S(u,v;m,n)=n(n-1)((Subscript[p, 02])-2(Subscript[p, 01])+(Subscript[p, 00])),

((∂^2)/(∂u∂v))(Subscript[|, u=0,v=0])S(u,v;m,n)=m n((Subscript[p, 11])-(Subscript[p, 10])-(Subscript[p, 01])+(Subscript[p, 00])).

Вычислив, согласно приведенным формулам, первые и вторые производные радиус-вектора поверхности в угловой точке, а также нормаль к поверхности в этой точке, можно по формулам дифференциальной геометрии вычислить гауссову и среднюю кривизну поверхности Безье в этой точки, а также главные направления и главные кривизны. Все эти величины будут выражаться через опорные точки (Subscript[p, 00]), (Subscript[p, 01]), (Subscript[p, 02]), (Subscript[p, 10]) и (Subscript[p, 20]).

Для произвольной точки S((SuperStar[u]),(SuperStar[v]);m,n) поступаем так:

1) делим поверхность в этой точке и вычисляем узлы (Subscript[p, 00]^(b d)), (Subscript[p, 01]^(b d)), (Subscript[p, 02]^(b d)), (Subscript[p, 10]^(b d)), (Subscript[p, 20]^(b d)).

2) вычисляем геометрические характеристики в данной точке через эти узлы.

ПРОГРАММА 2. Чайник нарисованный с помощью поверхности Безьу. Чайник состоит из заплат.

DynamicModule[{potPoints,pp2,potPatches,t,i,n,Bernstein,BezierSurface},potPoints={{1.4,0,2.4},{1.4,-0.784,2.4},{0.784,-1.4,2.4},{0,-1.4,2.4},{1.3375,0,2.53125},{1.3375,-0.749,2.53125},{0.749,-1.3375,2.53125},{0,-1.3375,2.53125},{1.4375,0,2.53125},{1.4375,-0.805,2.53125},{0.805,-1.4375,2.53125},{0,-1.4375,2.53125},{1.5,0,2.4},{1.5,-0.84,2.4},{0.84,-1.5,2.4},{0,-1.5,2.4},{-0.784,-1.4,2.4},{-1.4,-0.784,2.4},{-1.4,0,2.4},{-0.749,-1.3375,2.53125},{-1.3375,-0.749,2.53125},{-1.3375,0,2.53125},{-0.805,-1.4375,2.53125},{-1.4375,-0.805,2.53125},{-1.4375,0,2.53125},{-0.84,-1.5,2.4},{-1.5,-0.84,2.4},{-1.5,0,2.4},{-1.4,0.784,2.4},{-0.784,1.4,2.4},{0,1.4,2.4},{-1.3375,0.749,2.53125},{-0.749,1.3375,2.53125},{0,1.3375,2.53125},{-1.4375,0.805,2.53125},{-0.805,1.4375,2.53125},{0,1.4375,2.53125},{-1.5,0.84,2.4},{-0.84,1.5,2.4},{0,1.5,2.4},{0.784,1.4,2.4},{1.4,0.784,2.4},{0.749,1.3375,2.53125},{1.3375,0.749,2.53125},{0.805,1.4375,2.53125},{1.4375,0.805,2.53125},{0.84,1.5,2.4},{1.5,0.84,2.4},{1.75,0,1.875},{1.75,-0.98,1.875},{0.98,-1.75,1.875},{0,-1.75,1.875},{2.,0,1.35},{2.,-1.12,1.35},{1.12,-2.,1.35},{0,-2.,1.35},{2.,0,0.9},{2.,-1.12,0.9},{1.12,-2.,0.9},{0,-2.,0.9},{-0.98,-1.75,1.875},{-1.75,-0.98,1.875},{-1.75,0,1.875},{-1.12,-2.,1.35},{-2.,-1.12,1.35},{-2.,0,1.35},{-1.12,-2.,0.9},{-2.,-1.12,0.9},{-2.,0,0.9},{-1.75,0.98,1.875},{-0.98,1.75,1.875},{0,1.75,1.875},{-2.,1.12,1.35},{-1.12,2.,1.35},{0,2.,1.35},{-2.,1.12,0.9},{-1.12,2.,0.9},{0,2.,0.9},{0.98,1.75,1.875},{1.75,0.98,1.875},{1.12,2.,1.35},{2.,1.12,1.35},{1.12,2.,0.9},{2.,1.12,0.9},{2.,0,0.45},{2.,-1.12,0.45},{1.12,-2.,0.45},{0,-2.,0.45},{1.5,0,0.225},{1.5,-0.84,0.225},{0.84,-1.5,0.225},{0,-1.5,0.225},{1.5,0,0.15},{1.5,-0.84,0.15},{0.84,-1.5,0.15},{0,-1.5,0.15},{-1.12,-2.,0.45},{-2.,-1.12,0.45},{-2.,0,0.45},{-0.84,-1.5,0.225},{-1.5,-0.84,0.225},{-1.5,0,0.225},{-0.84,-1.5,0.15},{-1.5,-0.84,0.15},{-1.5,0,0.15},{-2.,1.12,0.45},{-1.12,2.,0.45},{0,2.,0.45},{-1.5,0.84,0.225},{-0.84,1.5,0.225},{0,1.5,0.225},{-1.5,0.84,0.15},{-0.84,1.5,0.15},{0,1.5,0.15},{1.12,2.,0.45},{2.,1.12,0.45},{0.84,1.5,0.225},{1.5,0.84,0.225},{0.84,1.5,0.15},{1.5,0.84,0.15},{-1.6,0,2.025},{-1.6,-0.3,2.025},{-1.5,-0.3,2.25},{-1.5,0,2.25},{-2.3,0,2.025},{-2.3,-0.3,2.025},{-2.5,-0.3,2.25},{-2.5,0,2.25},{-2.7,0,2.025},{-2.7,-0.3,2.025},{-3.,-0.3,2.25},{-3.,0,2.25},{-2.7,0,1.8},{-2.7,-0.3,1.8},{-3.,-0.3,1.8},{-3.,0,1.8},{-1.5,0.3,2.25},{-1.6,0.3,2.025},{-2.5,0.3,2.25},{-2.3,0.3,2.025},{-3.,0.3,2.25},{-2.7,0.3,2.025},{-3.,0.3,1.8},{-2.7,0.3,1.8},{-2.7,0,1.575},{-2.7,-0.3,1.575},{-3.,-0.3,1.35},{-3.,0,1.35},{-2.5,0,1.125},{-2.5,-0.3,1.125},{-2.65,-0.3,0.9375},{-2.65,0,0.9375},{-2.,-0.3,0.9},{-1.9,-0.3,0.6},{-1.9,0,0.6},{-3.,0.3,1.35},{-2.7,0.3,1.575},{-2.65,0.3,0.9375},{-2.5,0.3,1.125},{-1.9,0.3,0.6},{-2.,0.3,0.9},{1.7,0,1.425},{1.7,-0.66,1.425},{1.7,-0.66,0.6},{1.7,0,0.6},{2.6,0,1.425},{2.6,-0.66,1.425},{3.1,-0.66,0.825},{3.1,0,0.825},{2.3,0,2.1},{2.3,-0.25,2.1},{2.4,-0.25,2.025},{2.4,0,2.025},{2.7,0,2.4},{2.7,-0.25,2.4},{3.3,-0.25,2.4},{3.3,0,2.4},{1.7,0.66,0.6},{1.7,0.66,1.425},{3.1,0.66,0.825},{2.6,0.66,1.425},{2.4,0.25,2.025},{2.3,0.25,2.1},{3.3,0.25,2.4},{2.7,0.25,2.4},{2.8,0,2.475},{2.8,-0.25,2.475},{3.525,-0.25,2.49375},{3.525,0,2.49375},{2.9,0,2.475},{2.9,-0.15,2.475},{3.45,-0.15,2.5125},{3.45,0,2.5125},{2.8,0,2.4},{2.8,-0.15,2.4},{3.2,-0.15,2.4},{3.2,0,2.4},{3.525,0.25,2.49375},{2.8,0.25,2.475},{3.45,0.15,2.5125},{2.9,0.15,2.475},{3.2,0.15,2.4},{2.8,0.15,2.4},{0,0,3.15},{0,-0.002,3.15},{0.002,0,3.15},{0.8,0,3.15},{0.8,-0.45,3.15},{0.45,-0.8,3.15},{0,-0.8,3.15},{0,0,2.85},{0.2,0,2.7},{0.2,-0.112,2.7},{0.112,-0.2,2.7},{0,-0.2,2.7},{-0.002,0,3.15},{-0.45,-0.8,3.15},{-0.8,-0.45,3.15},{-0.8,0,3.15},{-0.112,-0.2,2.7},{-0.2,-0.112,2.7},{-0.2,0,2.7},{0,0.002,3.15},{-0.8,0.45,3.15},{-0.45,0.8,3.15},{0,0.8,3.15},{-0.2,0.112,2.7},{-0.112,0.2,2.7},{0,0.2,2.7},{0.45,0.8,3.15},{0.8,0.45,3.15},{0.112,0.2,2.7},{0.2,0.112,2.7},{0.4,0,2.55},{0.4,-0.224,2.55},{0.224,-0.4,2.55},{0,-0.4,2.55},{1.3,0,2.55},{1.3,-0.728,2.55},{0.728,-1.3,2.55},{0,-1.3,2.55},{1.3,0,2.4},{1.3,-0.728,2.4},{0.728,-1.3,2.4},{0,-1.3,2.4},{-0.224,-0.4,2.55},{-0.4,-0.224,2.55},{-0.4,0,2.55},{-0.728,-1.3,2.55},{-1.3,-0.728,2.55},{-1.3,0,2.55},{-0.728,-1.3,2.4},{-1.3,-0.728,2.4},{-1.3,0,2.4},{-0.4,0.224,2.55},{-0.224,0.4,2.55},{0,0.4,2.55},{-1.3,0.728,2.55},{-0.728,1.3,2.55},{0,1.3,2.55},{-1.3,0.728,2.4},{-0.728,1.3,2.4},{0,1.3,2.4},{0.224,0.4,2.55},{0.4,0.224,2.55},{0.728,1.3,2.55},{1.3,0.728,2.55},{0.728,1.3,2.4},{1.3,0.728,2.4},{0,0,0},{1.5,0,0.15},{1.5,0.84,0.15},{0.84,1.5,0.15},{0,1.5,0.15},{1.5,0,0.075},{1.5,0.84,0.075},{0.84,1.5,0.075},{0,1.5,0.075},{1.425,0,0},{1.425,0.798,0},{0.798,1.425,0},{0,1.425,0},{-0.84,1.5,0.15},{-1.5,0.84,0.15},{-1.5,0,0.15},{-0.84,1.5,0.075},{-1.5,0.84,0.075},{-1.5,0,0.075},{-0.798,1.425,0},{-1.425,0.798,0},{-1.425,0,0},{-1.5,-0.84,0.15},{-0.84,-1.5,0.15},{0,-1.5,0.15},{-1.5,-0.84,0.075},{-0.84,-1.5,0.075},{0,-1.5,0.075},{-1.425,-0.798,0},{-0.798,-1.425,0},{0,-1.425,0},{0.84,-1.5,0.15},{1.5,-0.84,0.15},{0.84,-1.5,0.075},{1.5,-0.84,0.075},{0.798,-1.425,0},{1.425,-0.798,0}}; (*это точки каркаса чайника, их 306*)
pp2=With[{avg=Mean[potPoints]},Map[(#-avg)&,potPoints]]; (*усреднили все точки и составили из них список*)
potPatches={{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16},{4,17,18,19,8,20,21,22,12,23,24,25,16,26,27,28},{19,29,30,31,22,32,33,34,25,35,36,37,28,38,39,40},{31,41,42,1,34,43,44,5,37,45,46,9,40,47,48,13},{13,14,15,16,49,50,51,52,53,54,55,56,57,58,59,60},{16,26,27,28,52,61,62,63,56,64,65,66,60,67,68,69},{28,38,39,40,63,70,71,72,66,73,74,75,69,76,77,78},{40,47,48,13,72,79,80,49,75,81,82,53,78,83,84,57},{57,58,59,60,85,86,87,88,89,90,91,92,93,94,95,96},{60,67,68,69,88,97,98,99,92,100,101,102,96,103,104,105},{69,76,77,78,99,106,107,108,102,109,110,111,105,112,113,114},{78,83,84,57,108,115,116,85,111,117,118,89,114,119,120,93},{121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136},{124,137,138,121,128,139,140,125,132,141,142,129,136,143,144,133},{133,134,135,136,145,146,147,148,149,150,151,152,69,153,154,155},{136,143,144,133,148,156,157,145,152,158,159,149,155,160,161,69},{162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177},{165,178,179,162,169,180,181,166,173,182,183,170,177,184,185,174},{174,175,176,177,186,187,188,189,190,191,192,193,194,195,196,197},{177,184,185,174,189,198,199,186,193,200,201,190,197,202,203,194},{204,204,204,204,207,208,209,210,211,211,211,211,212,213,214,215},{204,204,204,204,210,217,218,219,211,211,211,211,215,220,221,222},{204,204,204,204,219,224,225,226,211,211,211,211,222,227,228,229},{204,204,204,204,226,230,231,207,211,211,211,211,229,232,233,212},{212,213,214,215,234,235,236,237,238,239,240,241,242,243,244,245},{215,220,221,222,237,246,247,248,241,249,250,251,245,252,253,254},{222,227,228,229,248,255,256,257,251,258,259,260,254,261,262,263},{229,232,233,212,257,264,265,234,260,266,267,238,263,268,269,242},{270,270,270,270,279,280,281,282,275,276,277,278,271,272,273,274},{270,270,270,270,282,289,290,291,278,286,287,288,274,283,284,285},{270,270,270,270,291,298,299,300,288,295,296,297,285,292,293,294},{270,270,270,270,300,305,306,279,297,303,304,275,294,301,302,271}};(*чайник сосотоит из заплат, их 32 штуки*)
Bernstein[i_Integer,n_Integer,t_Symbol]:=Binomial[n,i] t^i (1-t)^(n-i); (*полиномы Бернштейна*)
Bernstein[n_Integer,t_Symbol]:=Map[Bernstein[#,n,t]&,Range[0,n]]; (*базис -- список всех полиномов Бернштейна*)
BezierSurface[pts_,u_,v_,m_]:=Module[{dim=Dimensions[pts],avg=Mean[Flatten[pts,1]]},Bernstein[dim[[1]]-1,u].Map[Bernstein[dim[[2]]-1,v].#&,Map[(#)&,pts]]+m*avg] ;(*формула поверхности Безье*)
Manipulate[ParametricPlot3D[Evaluate[BezierSurface[#,s,t,m]&/@(Transpose[Partition[Part[pp2,#],4]]&/@potPatches)],(*делаем матрицу опорных точек и для них рисуем поверхность Безье*){s,0,1},{t,0,1},Axes->None,Boxed->False,Evaluated->False,PlotPoints->ControlActive[{3,3},{8,8}],PlotStyle->Directive[Specularity[White,50],Yellow],MaxRecursion->0,Mesh->None,PerformanceGoal->"Speed",ImageSize->{450,350}],{{m,0,"explode"},0,1},SaveDefinitions->True]
]

2.1.1. В-ПОВЕРХНОСТИ

Пусть задана матрица опорных точек ((Subscript[p, i j])), 1≤i≤n, 1≤j≤m и соответствующая матрица весов ((Subscript[w, i j])). Фиксируем числа p и q – порядки B-сплайнов по двум координатным направлениям u и v соответственно.  Тогда B-ПОВЕРХНОСТЬ порядка (p,q), построенная по матрице опорных точек ((Subscript[p, i j])) с весами  ((Subscript[w, i j])) с помощью двух расширенных множества узлов в направлениях u и v соответственно, определяется формулой

r (u, v)=(((∑_(i=1)^n)(∑_(j=1)^m)(Subscript[N, i,q])(v)(Subscript[N, j,p])(u)(Subscript[w, i j])(Subscript[p, i j]))/((∑_(i=1)^n)(∑_(j=1)^m)(Subscript[N, i,q])(v)(Subscript[N, j,p])(u)(Subscript[w, i j]))).

Если ввести новые обозначения для входящих в эту формулу произведений B-сплайнов: (Subscript[N, i j,p q])(u,v)=(Subscript[N, i,q])(v)(Subscript[N, j,p])(u), то вышеприведенная формула параметрического задания В-поверхности перепишется в виде:

r (u, v)=(((∑_(i=1)^n)(∑_(j=1)^m)(Subscript[N, i j,p q])(u,v)(Subscript[w, i j])(Subscript[p, i j]))/((∑_(i=1)^n)(∑_(j=1)^m)(Subscript[N, i j,p q])(u,v)(Subscript[w, i j]))).

Если расширенные множества узлов в u и v  – направлениях имеют вид:
(Subscript[u, 1])=...=(Subscript[u, p]); (Subscript[u, p+1]),...,(Subscript[u, n]); (Subscript[u, n+1])=...=(Subscript[u, n+k;])
(Subscript[v, 1])=...=(Subscript[v, q]); (Subscript[v, q+1]),...,(Subscript[v, n]); (Subscript[v, n+1])=...=(Subscript[v, n+q,])
то область изменения параметров (u,v) на B-поверхности будет следующей:
 (Subscript[u, p]) ≤ u ≤ (Subscript[u, n+1]);  (Subscript[v, q])≤v≤(Subscript[v, n+1]).

УПРАЖНЕНИЕ. Наглядно изучить влияние опорных точек и расширенного множества узлов на В-поверхность, порождаемую МАТЕМАТИКОЙ с помощью следующей программы. Поменять в программе автоматический выбор расширенных множеств узлов на прямое задание этих узлов. Изучить влияние расширенных множеств узлов на В-поверхность.

ПРОГРАММА 3. В-поверхность. Для обеих поверхностей задан одинаковый остов. В первом случае поверхность ``замыкается'' только по оси x, во втором случаем по обеим осям.

pts=Table[{Cos[2 Pi u/6] Cos[v],Sin[2 Pi u/6] Cos[v],v},{u,6},{v,-1,1,1/2}];
f=BSplineFunction[pts,SplineClosed->{True,False},SplineKnots->Automatic];
g=BSplineFunction[pts,SplineClosed->{True,True},SplineKnots->Automatic];

Show[Graphics3D[{PointSize[Large],Red,Map[Point,pts]}],Graphics3D[{Gray,Dashed,Line[pts],Gray,Line[Transpose[pts]]}],ParametricPlot3D[f[u,v],{u,0,1},{v,0,1}]]

Show[Graphics3D[{PointSize[Large],Red,Map[Point,pts]}],Graphics3D[{Gray,Dashed,Line[pts],Gray,Line[Transpose[pts]]}],ParametricPlot3D[g[u,v],{u,0,1},{v,0,1}]]

ПРОГРАММА 4. Моделирование кузова автомобиля с помощью B-поверхности.

DynamicModule[{colour,carbody,deg,mesh},
Manipulate[
Graphics3D[{
colour,Specularity[White,50],
BSplineSurface[carbody,SplineDegree->Floor[deg]],
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[{{deg,{1,1},"Порядок"},{1,1},{14,7},ImageSize->{140,70}}],Dynamic[Floor[deg]]},Spacer[10]],
Row[{Control[{{colour,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}},
{{0,-2.25,0},{0,-2.25,1.75},{0,-2,1.75},{0,0,1.75},{0,2,1.75},{0,2.25,1.75},{0,2.25,0}},
{{1,-2.25,0},{1,-2.25,1.75},{1,-2,1.75},{1,0,1.75},{1,2,1.75},{1,2.25,1.75},{1,2.25,0}},
{{1.01,-2.25,1.5},{1.01,-2.25,1.75},{1.01,-2,1.75},{1.01,0,1.75},{1.01,2,1.75},{1.01,2.25,1.75},{1.01,2.25,1.5}},
{{2.99,-2.25,1.5},{2.99,-2.25,1.75},{2.99,-2,1.75},{2.99,0,1.75},{2.99,2,1.75},{2.99,2.25,1.75},{2.99,2.25,1.5}},
{{3,-2.25,0},{3,-2.25,1.75},{3,-2,1.75},{3,0,1.75},{3,2,1.75},{3,2.25,1.75},{3,2.25,0}},
{{4,-2.25,0},{4,-2.25,1.75},{4,-2,1.75},{4,0,1.75},{4,2,1.75},{4,2.25,1.75},{4,2.25,0}},
{{5,-2.25,0},{5,-2.25,1.75},{5,-2,3.5},{5,0,3.5},{5,2,3.5},{5,2.25,1.75},{5,2.25,0}},
{{8,-2.25,0},{8,-2.25,1.75},{8,-2,3.5},{8,0,3.5},{8,2,3.5},{8,2.25,1.75},{8,2.25,0}},
{{8.01,-2.25,1.5},{8.01,-2.25,1.75},{8.01,-2,3.5},{8.01,0,3.5},{8.01,2,3.5},{8.01,2.25,1.75},{8.01,2.25,1.5}},
{{9,-2.25,1.5},{9,-2.25,2},{9,-2,2},{9,0,2},{9,2,2},{9,2.25,2},{9,2.25,1.5}},
{{9.99,-2.25,1.5},{9.99,-2.25,2},{9.99,-2,2},{9.99,0,2},{9.99,2,2},{9.99,2.25,2},{9.99,2.25,1.5}},
{{10,-2.25,0},{10,-2.25,2},{10,-2,2},{10,0,2},{10,2,2},{10,2.25,2},{10,2.25,0}},
{{11,-2.25,0},{11,-2.25,2},{11,-2,2},{11,0,2},{11,2,2},{11,2.25,2},{11,2.25,0}},
{{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}}
})
]

ПРОГРАММА 5. Интерполирование заданного рисунка с помощью поверхности Безье и B-поверхности.

DynamicModule[{f,colfun,prevn=5,prevapprox=True,start,locators,grid,order,approx,images},
images={,};
Manipulate[
(* Создали список, который задает рисунок, и для этого списка построили B-функцию, т.е. мы ``синтерполировали'' рисунок*)
colfun=BSplineFunction[ImageData[images[[ii]]],SplineDegree->1];
(* Первоначальная сетка *)
If[n=!=prevn,pts=Flatten[Table[{i,j},{i,0,1,1/(n-1)},{j,0,1,1/(n-1)}],1];If[approx,order=n-1,order=2];prevn=n];
(* Выбираем порядок интерполяции *)
If[approx=!=prevapprox,If[approx,order=n-1,order=2];prevapprox=approx];
If[approx,
f=BSplineFunction[Partition[pts,n],SplineDegree->order,Method->{"Extrapolation"->"Clamp"}],f=Interpolation[Flatten[Table[{{N[1/(n-1)]i,N[1/(n-1)]j},pts[[n i+j+1]]},{i,0,n-1},{j,0,n-1}],1],Method->"Spline",InterpolationOrder->order]
];
(*Рисуем картинки*)
ParametricPlot[f[u,v],{u,0,1},{v,0,1},ColorFunction->(colfun[1-#4,#3]&),MaxRecursion->0,PlotPoints->ControlActive[20,100],Mesh->None,Epilog->{If[grid,{Gray,Opacity[.7],Line[Partition[pts,n]],Line[Transpose[Partition[pts,n]]]},{}],If[locators,Locator/@pts,{}]},ImageSize->400,PlotRange->{{0,1},{0,1}},PlotRangePadding->Scaled[.03],FrameTicks->None],
(*элементы управления*)
{{pts,Flatten[Table[{i,j},{i,0,1,.25},{j,0,1,.25}],1]},{0,0},{1,1},Locator,Appearance->None},
{{ii,1,"Объкт интерполирования"},{1->"Медведь",2->"Шахматная доска"},ControlType->RadioButton},
Delimiter,
{{n,4,"Размер сетки"},4,8,1,ControlType->RadioButton},
{{approx,True,"Метод"},{True->"Поверхность Безье",False->"B-поверхность"},ControlType->RadioButton},
{{order,3,"Порядок B-поверхности"},Dynamic[RadioButtonBar[Dynamic[order],Range[2,n-1,1],Enabled->Not[approx]]]&},
Delimiter,
Row[{Control[{{grid,True,"Показать сетку"},{True,False}}],
Control[{{locators,True,"Показать локаторы"},{True,False}}],
Control[{{start,"",""},Button["Первоначальная сетка",pts=Flatten[Table[{i,j},{i,0,1,1/(n-1)},{j,0,1,1/(n-1)}],1]]&}]},Spacer[10]],
(* Код, соответствующий второму рисунку,следующий: Image[Table[Boole[Xor[EvenQ[Quotient[i-1,20]],EvenQ[Quotient[j-1,20]]]],{i,200},{j,200},{k,3}]] *)
Initialization:>(images={,};)]
]

2.2. ПОВЕРХНОСТИ, ЗАТЯГИВАЮЩИЕ ЗАДАННЫЕ ГРАНИЧНЫЕ КОНТУРЫ

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

2.2.1. ЛИНЕЙЧАТЫЕ ПОВЕРХНОСТИ

Определение. ЛИНЕЙЧАТАЯ ПОВЕРХНОСТЬ – геометрическое место точек отрезков прямых, соединяющих соответствующие друг другу точки двух заданных  кривых.  Пусть (Subscript[ρ, 1])(t), (Subscript[t, 0])≤t≤(Subscript[t, 1]), и (Subscript[ρ, 2])(w), (Subscript[w, 0])≤w≤(Subscript[w, 1]), – две заданные кривые. Тогда определяемая ими линейчатая поверхность имеет вид
ρ(u,v)=(1-v)(Subscript[ρ, 1])(t)+(Subscript[vρ, 2])(w), 
где t=(Subscript[t, 0])(1-u)+(Subscript[t, 1])u, w=(Subscript[w, 0])(1-u)+(Subscript[w, 1])u, 0≤u,v≤1.

Если обе кривые (Subscript[ρ, 1]) и (Subscript[ρ, 2]) замкнуты, то линейчатая поверхность замкнута по u. По параметру v линейчатая поверхность всегда незамкнута.

Если обе кривые (Subscript[ρ, 1]) и (Subscript[ρ, 2]) – отрезки прямых, то линейчатая поверхность называется БИЛИНЕЙНОЙ поверхностью и определяется только своими четырьмя угловыми точками (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]), (Subscript[p, 4]):

r(u,v)=(1-u)(1-v)(Subscript[p, 1])+u(1-v)(Subscript[p, 2])+(1-u) v (Subscript[p, 3])+u v (Subscript[p, 4]).

2.2.2. СЕКТОРИАЛЬНЫЕ И ТРЕУГОЛЬНЫЕ ПОВЕРХНОСТИ

Определение. СЕКТОРИАЛЬНАЯ ПОВЕРХНОСТЬ – частный случай линейчатой поверхности, когда одна из кривых (Subscript[ρ, 1]) и (Subscript[ρ, 2]) вырождается в точку:
r(u,v)=(1-v)ρ(u)+v p,
где ρ(u) – кривая, p – точка, (Subscript[u, 1])≤u≤(Subscript[u, 2]), 0≤v≤1.

Секториальная билинейная поверхность называется ТРЕУГОЛЬНОЙ ПОВЕРХНОСТБЮ.

2.2.3. ПОВЕРХНОСТИ КУНСА

2.2.3.1. ЛИНЕЙНЫЕ ПОВЕРХНОСТИ КУНСА

Определение. Пусть заданы четыре кривые, образующие криволинейный четырехугольник с вершинами в точках (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]), (Subscript[p, 4]):
(Subscript[γ, a])=r(u,0), 0≤u≤1,r(0,0)=(Subscript[p, 1]), r(1,0)=(Subscript[p, 2]);
(Subscript[γ, b])=r(0,v), 0≤v≤1,r(0,0)=(Subscript[p, 1]), r(0,1)=(Subscript[p, 3]);
(Subscript[γ, c])=r(1,v), 0≤v≤1,r(1,0)=(Subscript[p, 2]), r(1,1)=(Subscript[p, 4]);
(Subscript[γ, d])=r(u,1),0≤u≤1,r(0,1)=(Subscript[p, 3]), r(1,1)=(Subscript[p, 4]).
Обозначим: (Subscript[α, 0])(t)=1-t, (Subscript[α, 1])(t)=t. Функции (Subscript[α, 0])(t), (Subscript[α, 1])(t) будем называть ФУНКЦИЯМИ СМЕЩЕНИЯ линейной поверхности Кунса.
Сама ЛИНЕЙНАЯ ПОВЕРХНОСТЬ КУНСА определяется следующей формулой:
r(u,v)=(Subscript[α, 0])(v)r(u,0)+(Subscript[α, 0])(u)r(0,v)+(Subscript[α, 1])(v)r(u,1)+(Subscript[α, 1])(u)r(1,v)-(Subscript[α, 0])(u)(Subscript[α, 0])(v)r(0,0)-(Subscript[α, 1])(u)(Subscript[α, 0])(v)r(1,0)-(Subscript[α, 0])(u)(Subscript[α, 1])(v)r(0,1)-(Subscript[α, 1])(u)(Subscript[α, 1])(v)r(1,1).

ЗАМЕЧАНИЯ.

1) Линейная поверхность Кунса, вообще говоря, не является линейчатой поверхностью. Название ``линейная'' обусловлено тем, что в ее определении используются линейные функции смещения (Subscript[a, 0]) и (Subscript[α, 1]).

2) Параметрическая область поверхности Кунса может иметь вид прямоугольника, не обязательно совпадающего с квадратом [0,1]x[0,1].

ПРИМЕР. Зададим с помощью поверхности Кунса область D плоскости, ограниченную эллипсом ((x^2)/(a^2))+((y^2)/(b^2))=1 с полуосями a>0, b>0. Представим границу области D в виде четырех кривых – дуг эллипса, стыкующихся в точках (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]), (Subscript[p, 4]). Соответствующая линейная поверхность Кунса имеет вид
r(u,v)=a(cos u+cos v-1+((2u(1-cos v-sin v))/π)+((2v(1-cos u-sin u))/π))(Subscript[e, 1])+b(sin u-sin v-((2u(1-cos v-sin v))/π)+((2v(1-cos u-sin u))/π))(Subscript[e, 2]), 0≤u,v≤(π/2).

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

ПРОГРАММА 6. Поверхность Кунса, построенная по кривым, лежащим на заданных поверхностях. Кривые раскрашены красным, зеленым, желтым и оранжевым.

DynamicModule[{r1=30., r2=30., r3=30.,d=50.,u,v,m1,m2,m3,m4,f1,f2,f3,f4,rfel,α,MV,MU,lam1,lam2,rfk},
r1=30.; r2=30.; r3=30.; d=50.;
lam1=2.; lam2=1.;
f1[m1_]:={r1-r1*Cos[u*Pi],lam1*d*v/2,m1*(Sin[u*Pi]*Sin[v*Pi])^4} ;
f2[m2_]:={r2-r2*Cos[u*Pi],d-lam2*d/2+lam2*d*v/2,m2*(Sin[u*Pi])^2};
f3[m3_]:={0.8*r3*u,d*v,m3*(Sin[v*3*Pi])^2};
f4[m4_]:={1.6*r1+0.4*r1*u+(v^4-2*v^3+v^2)*100*m4*u,d*v,0};
Manipulate[
Show[{ParametricPlot3D[Evaluate@f1[m1],{u,0,1},{v,-1,0},Mesh->{24,24},PlotRange->All],
ParametricPlot3D[Evaluate@{f1[m1]/.v->0},{u,0,1},Mesh->{24,24},PlotRange->All,PlotStyle->{Red}],
         ParametricPlot3D[Evaluate@f2[m2],{u,0,1},{v,1,2} ,Mesh->{24,2},PlotRange->All],
ParametricPlot3D[Evaluate@{f2[m2]/.v->1},{u,0,1},Mesh->{24,24},PlotRange->All,PlotStyle->{Green}],
ParametricPlot3D[Evaluate@f3[m3],{u,-1,0},{v,0,1}, Mesh->{2,24},PlotRange->All],
ParametricPlot3D[Evaluate@{f3[m3]/.u->0},{v,0,1},Mesh->{24,24},PlotRange->All,PlotStyle->{Yellow}],
ParametricPlot3D[Evaluate@f4[m4],{u,1,2},{v,0,1},Mesh->{2,24},PlotRange->All],
ParametricPlot3D[Evaluate@{f4[m4]/.u->1},{v,0,1},Mesh->{24,24},PlotRange->All,PlotStyle->{Orange}],
ParametricPlot3D[Evaluate@rfel[m1,m2,m3,m4],{u,0,1},{v,0,1},Mesh->{24,24}, PlotRange->All]},
ImageSize->{500,380}, Axes->None, Boxed->False,PlotRange->{{-30,135},{-50,100},{-30,30}}
],
{{m1,10,"``Ухабистая'' поверхность"},-15., 15.,Appearance-> "Labeled"},
{{m2,20,"Цилиндрическая поверхность"},-30, 30,Appearance-> "Labeled"},
{{m3,10,"Волновая поверхность"},-10, 10, Appearance-> "Labeled"},
{{m4,1,"Плоская поверхность"},-5, 5,Appearance-> "Labeled"}, ControlPlacement->Top,
Initialization:>(
r1=30.; r2=30.; r3=30.; d=50.;
lam1=2.; lam2=1.;
f1[m1_]:={r1-r1*Cos[u*Pi],lam1*d*v/2,m1*(Sin[u*Pi]*Sin[v*Pi])^4} ;
f2[m2_]:={r2-r2*Cos[u*Pi],d-lam2*d/2+lam2*d*v/2,m2*(Sin[u*Pi])^2};
f3[m3_]:={0.8*r3*u,d*v,m3*(Sin[v*3*Pi])^2};
f4[m4_]:={1.6*r1+0.4*r1*u+(v^4-2*v^3+v^2)*100*m4*u,d*v,0};
α[t_]:=t;
MV[v_]:={1-α[v],α[v]}; 
MU[u_]:={1-α[u],α[u]}; (* наши функции α *)
Table[rfk[m1,m2,m3,m4][i],{i,3}];
Do[rfk[m1_,m2_,m3_,m4_][i]=
MU[u].{f3[m3][[i]]/.u->0,f4[m4][[i]]/.u->1}+
MV[v].{f1[m1][[i]]/.v->0,f2[m2][[i]]/.v->1}-
MU[u].{{f3[m3][[i]]/.{u->0,v->0},f3[m3][[i]]/.{u->0,v->1}},{f4[m4][[i]]/.{u->1,v->0},f4[m4][[i]]/.{u->1,v->1}}}.MV[v]//Together;
,{i,3}];
rfel[m1_,m2_,m3_,m4_]={rfk[m1,m2,m3,m4][1],rfk[m1,m2,m3,m4][2],rfk[m1,m2,m3,m4][3]};
)
]
]

2.2.3.2. ОБОБЩЕННЫЕ ПОВЕРХНОСТИ КУНСА

При стыковке линейных поверхностей Кунса по граничным кривым производные в направлении, ортогональном границе, претерпевают разрыв. Чтобы гладко склеивать поверхности Кунса необходимо строить их таким образом, чтобы они имели на границе заданные (согласованные) производные до k-го порядка включительно.  Производные в касательном направлении определяются уравнениями самих граничных кривых, а производные в трансверсальном к границе направлении необходимо задавать дополнительно.

Определение. ОБОБЩЕННАЯ ПОВЕРХНОСТЬ КУНСА n-ГО ПОРЯДКА называется поверхность

r(u,v)=(∑_(i=0)^1)(∑_(k=0)^n)(Subscript[α, k,i])(u)(Subscript[r, u]^((k)))(i,v)+(∑_(j=0)^1)(∑_(k=0)^n)(Subscript[α, k,j])(v)(Subscript[r, v]^((k)))(u,j)-(∑_(i=0)^1)(∑_(j=0)^1)(∑_(k=0)^n)(∑_(l=0)^n)(Subscript[α, k,i])(u)(Subscript[α, l,j])(v)(Subscript[r, u v]^((k,l)))(i,j), 0≤u,v≤1, 
задаваемая четырьмя граничными кривыми r(u,0), r(0,v), r(1,v),r(u,1), стыкующимися в точках (Subscript[p, 1])=r(0,0), (Subscript[p, 2])=r(1,0),
(Subscript[p, 3])=r(0,1) и  (Subscript[p, 4])=r (1, 1) , а также чистыми частными производными в трансверсальном направлении на граничных кривых до n-го порядка включительно и смешанными частными производными в угловых точках (Subscript[p, 1]), (Subscript[p, 2]), (Subscript[p, 3]), (Subscript[p, 4]) до 2n-го порядка включительно, а также – ОБОБЩЕННЫМИ ФУНКЦИЯМИ СМЕЩЕНИЯ (Subscript[α, k,i])(t), где 0≤k≤n, i=0,1, в качестве которых могут быть взяты произвольные гладкие функции, удовлетворяющие соотношениям:

(((d^m)(Subscript[α, k,i])(t))/(d (t^m)))(Subscript[|, t=j])=(Subscript[δ, i]^j)(Subscript[δ, k]^m). Здесь 0≤k,m≤n.

ЗАМЕЧАНИЕ. В качестве обобщенных функций смещения всегда можно выбрать многочлены.

ПРИМЕР. КУБИЧЕСКАЯ ПОВЕРХНОСТЬ КУНСА – это обобщенная поверхность Кунса, построенная для

(Subscript[α, 0,0])(t)=1-3(t^2)+2(t^3), (Subscript[α, 0,1])(t)=3(t^2)-2(t^3), (Subscript[α, 1,0])(t)=t-2(t^2)+(t^3), (Subscript[α, 1,1])(t)=-(t^2)+(t^3) и k,l=1.

УПРАЖНЕНИЕ. Написать на МАТЕМАТИКЕ программу построения кубической поверхности Кунса.

2.2.4. ПОВЕРХНОСТИ, ПОСТРОЕННЫЕ ПО ОСТОВУ ИЗ КРИВЫХ

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

2.2.4.1. ПОВЕРХНОСТИ ЭРМИТА

Аналогично сплайнам Эрмита строятся  ПОВЕРХНОСТИ ЭРМИТА.

Пусть имеется семейство кривых: (Subscript[r, 0])(u)=r(u,0),..., (Subscript[r, i])(u)=r(u,i),..., (Subscript[r, n])(u)=r(u,n).
Без ограничения общности будем считать, что параметр u для всех кривых пробегает один и тот же отрезок: (Subscript[u, 0])≤u≤(Subscript[u, 1]). Тогда поверхность Эрмита, затягивающая остов кривых r(u,i), 0≤i≤n, определяется формулой аналогичной формуле для сплайна Эрмита:
r(u,v)=r(u,j)(Subscript[α, 0])(w)+r(u,j+1)(Subscript[α, 1])(w)+(Subscript[r, v])(u,j)(Subscript[β, 0])(w)+(Subscript[r, v])(u,j+1)(Subscript[β, 1])(w), если j≤v≤j+1. Здесь w=v-j – местный параметр на данном участке поверхности (0 ≤ w ≤ 1), (Subscript[r, v])(u,j) и (Subscript[r, v])(u,j+1) – производные в трансверсальном направлении на заданном остове кривых, (Subscript[α, 0])(w)=1-3(w^2)+2(w^3), (Subscript[α, 1])(w)=3(w^2)-2(w^3), (Subscript[β, 0])(w)=w-2(w^2)+(w^3), (Subscript[β, 1])(w)=-(w^2)+(w^3). Параметрическая область поверхности Эрмита имеет
вид (Subscript[u, 0])≤u≤(Subscript[u, 1]),0 ≤ v ≤ n.

Если трансверсальные производные (Subscript[r, v])(u,i) на линиях остова не заданы, то их, аналогично случаю сплайна Эрмита, можно определить, например так:
(Subscript[r, v])(u,j)=((r(u,j+1)-r(u,j-1))/2),1 ≤j≤ n-1,
(Subscript[r, v])(u,0)=2r(u,1)-2r(u,0)-(Subscript[r, v])(u,1),
(Subscript[r, v])(u,n)=2r(u,n)-2r(u,n-1)-(Subscript[r, v])(u,n-1).

При таком определении третьи производные на кривых (Subscript[r, 0]) и (Subscript[r, n]) в трансверсальном направлении обратятся в ноль.

Сами исходные кривые (Subscript[r, 0]),...,(Subscript[r, n]), образующие остов поверхности Эрмита, могут быть построены как кривые Безье или B-кривые. При этом получаются поверхности класса (C^1), имеющие на двух участках границы управляемые производные в трансверсальном направлении (а именно, на
участках r(u,0) и r(u,n)). На участках границы r((Subscript[u, 0]),v) и r((Subscript[u, 1]),v) производные, вообще говоря, управляются лишь приблизительно – с помощью выбора кривых остова (Subscript[r, i])(u).

2.2.4.2. ПРИМЕНЕНИЕ ПОВЕРХНОСТЕЙ ЭРМИТА: ПОВЕРХНОСТЬ ПЕРЕХОДА

Примером применения поверхностей Эрмита может служить задача объединения двух поверхностей в одну общую поверхность путем соединения их краев с помощью ПОВЕРХНОСТИ ПЕРЕХОДА.

Пусть даны две поверхности M и N, которые надо сопрячь, соединив поверхностью перехода край r=(Subscript[r, M])(t,(Subscript[u, 1]^M)) первой поверхности и край r=(Subscript[r, N])(t,(Subscript[u, 0]^N)) второй поверхности. Будем считать, что оба сопрягаемых края отнесены к одному и тому же параметру t, изменяющемуся в пределах
(Subscript[t, 0])≤t≤(Subscript[t, 1]). Второй параметр на поверхностях M и N обозначим (в обоих случаях) через u. Пусть параметрические пространства поверхностей имеют вид: (Subscript[t, 0])≤t≤(Subscript[t, 1]), (Subscript[u, 0]^M)≤u≤(Subscript[u, 1]^M)для поверхности M:r=(Subscript[r, M])(t,u) и (Subscript[t, 0])≤t≤(Subscript[t, 1]), (Subscript[u, 0]^N)≤u≤(Subscript[u, 1]^N)для поверхности N:r=(Subscript[r, N])(t,u). Без ограничения общности будем считать, что (Subscript[u, 1]^M)<(Subscript[u, 0]^N).  Тогда поверхность перехода от M к N может быть описана формулой поверхности Эрмита:

r(t,w)=(Subscript[r, M])(t,(Subscript[u, 1]^M))(1-3(w^2)+2(w^3))+(Subscript[r, N])(t,(Subscript[u, 0]^N))(3(w^2)-2(w^3))+(∂/(∂w))(Subscript[|, u=(Subscript[u, 1]^M)])(Subscript[r, M])(t,u)(w-2(w^2)+(w^3))+(∂/(∂w))(Subscript[|, u=(Subscript[u, 0]^N)])(Subscript[r, N])(t,u)(-(w^2)+(w^3)),

где w=((u-(Subscript[u, 1]^M))/((Subscript[u, 0]^N)-(Subscript[u, 1]^M))) – местный параметр по направлению u, (Subscript[t, 0])≤t≤(Subscript[t, 1]), 0≤w≤1. Данная поверхность Эрмита гладко сопрягает поверхности M и N по заданным краям. Класс гладкости сопряжения – (C^1).

ПРОГРАММА 7. Приклейка ручки с помощью поверхности перехода.

DynamicModule[{v,r1,r2,s1,s2,dr1,dr2,ds1,ds2,r,s},
r1[u_,v_]:={-3+(7+1.5Cos[2π u])Sin[2π v+π],1.5Sin[2π u],(7+1.5Cos[2π u])Cos[2π v+π]}; (*тор, u и v в пределах [0,1]*)
r2[u_,v_]:={5Cos[2π u]*Sin[π (v-1)],5Sin[2π u]*Sin[π (v-1)],5Cos[π (v-1)]}; (*сфера, u в пределах [0,1], v в пределах [1,2]*)

dr1[t_]:=D[r1[t,v], v]*(7/6-(1/2-1/10))/.v->1/2-1/10;
dr2[t_]:=D[r2[t,v], v]*(7/6-(1/2-1/10))/.v->7/6;
r[t_,w_]:=r1[t,1/2-1/10](1-3(w^2)+2(w^3))+r2[t,7/6](3(w^2)-2(w^3))+dr1[t](w-2(w^2)+(w^3))+dr2[t](-(w^2)+(w^3));

s1[u_,v_]:=r2[u,v+2]; (* u в пределах [0,1], v в пределах [-1,0]*)
s2[u_,v_]:=r1[u,v];
ds1[t_]:=D[s1[t,v], v]*(1/10+1/6)/.v->-1/6;
ds2[t_]:=D[s2[t,v], v]*(1/10+1/6)/.v->1/10;
s[t_,w_]:=s1[t,-1/6](1-3(w^2)+2(w^3))+s2[t,1/10](3(w^2)-2(w^3))+ds1[t](w-2(w^2)+(w^3))+ds2[t](-(w^2)+(w^3));

Manipulate[
Show[{
ParametricPlot3D[r1[u,v],{u,0,1},{v,1/10,1/2-1/10},PlotStyle-> FaceForm[Orange,Blue]],
ParametricPlot3D[r2[u,v],{u,0,1},{v,1+1/6,1+5/6},PlotStyle-> FaceForm[Orange,Blue]],
ParametricPlot3D[r[t,w],{t,0,1},{w,0,w1},PlotStyle-> FaceForm[Orange,Blue]],
ParametricPlot3D[s[t,w],{t,0,1},{w,0,w2},PlotStyle-> FaceForm[Orange,Blue]]
},PlotRange->{{-12,10},{-10,10},{-10,10}},Axes->False,Boxed->False
],{{w1,0.001},0.001,1},{{w2,0.001},0.001,1},TrackedSymbols->{w1,w2},SaveDefinitions->True]
]

2.2.4.3. ПОВЕРХНОСТИ ЛАГРАНЖА

Аналогично интерполяционной кривой Лагранжа, можно построить ПОВЕРХНОСТЬ ЛАГРАНЖА, затягивающую заданный остов кривых, см. Поверхности Эрмита. Эта поверхность описывается радиус-вектором:
r(u,v)=(Subscript[Σ, i=0]^n) (Subscript[L, i])(v) r(u,i), (Subscript[u, 0])≤u≤(Subscript[u, 1]),0 ≤ v ≤ n, где

(Subscript[L, i])(v)=(((Subscript[Π, j=0, j≠i]^n)(v-j))/((Subscript[Π, j=0, j≠i]^n)(i-j))).
Поверхность Лагранжа имеет те же недостатки, что и интерполяционный полином Лагранжа – необоснованные изгибы при неравномерном расположении кривых остова в пространстве. По параметру v поверхность будет бесконечно гладкой и представляет собой многочлен.  Гладкость по параметру u зависит от гладкости заданных кривых остова.

2.2.4.4. ПОВЕРХНОСТИ ГОРДОНА

Пусть задан остов поверхности, состоящий из двух семейств кривых в обоих параметрических направлениях, которые образуют 	СЕТКУ – т.е. каждая кривая одного семейства пресекает в одной и только одной точке каждую кривую другого семейства. Различные кривые одного и того же
семейства не пересекаются. Обозначим заданную таким образом сетку кривых через

r=r(u,0), ...,r=r(u,n), (первое семейство);
r=r(0,v), ...,r=r(m,v), (второе семейство).

Здесь 0≤u≤m,0 ≤ v ≤ n. Тогда радиус-вектор точки пересечения i-ой кривой первого семейства с j-ой кривой второго семейства будет равен r(i,j), 0≤i≤m,0 ≤ j ≤ n.

Опрделение. ПОВЕРХНОСТЬЮ ГОРДОНА, затягивающей сетку кривых, называется поверхность:

r(u,v)=(Subscript[Σ, j=0]^n)(Subscript[L, j]^n)(v)r(u,j)+(Subscript[Σ, i=0]^m)(Subscript[L, i]^m)(u)r(i,v)-(Subscript[Σ, j=0]^n)(Subscript[Σ, i=0]^m)(Subscript[L, i]^m)(u)(Subscript[L, j]^n)(v)r(i,j), 0 ≤ u ≤ m, 0 ≤ v ≤ n, где
(Subscript[L, i]^m)(u)=(((Subscript[Π, r=0, r≠i]^m)(u-r))/((Subscript[Π, r=0, r≠i]^m)(i-r))), (Subscript[L, j]^n)(v)=(((Subscript[Π, s=0, s≠j]^n)(v-s))/((Subscript[Π, s=0, s≠j]^n)(j-s))) – коэффициенты Лагранжа.
ОБОБЩЕННАЯ ПОВЕРХНОСТЬ ГОРДОНА задается формулой, участвующей в определении поверхности Гордона, в которой вместо коэффициентов Лагранжа используются произвольные (C^1)-функции (Subscript[α, k])(w), удовлетворяющие условию (Subscript[α, k])(w)=1, если w=k, и (Subscript[α, k])(w)=0, если w≠k.

Например, в качестве (Subscript[α, k])(w) можно взять кубические сплайны. Гладкость обобщенной поверхности Гордона зависит от гладкости используемых функций (Subscript[α, k])(w). Поверхности Гордона, построенные по коэффициентам Лагранжа, имеют гладкость (C^∞), но могут образовывать необоснованные изгибы.

2.2.4.5. ПОВЕРХНОСТИ, ЗАТЯГИВАЮЩИЕ  СЕТКУ КРИВЫХ ЗАПЛАТАМИ КУНСА

Остов, состоящий из сетки кривых можно затянуть поверхностью, состоящей из заплат, затягивающих каждую отдельную ячейку сетки и гладко сопрягающихся между собой. В качестве таких заплат могут выступать поверхности Кунса. Получается следующая поверхность
r(u,v)=(({{(Subscript[α, 0])(t), (Subscript[α, 1])(t), (Subscript[β, 0])(t), (Subscript[β, 1])(t)}}))(({{r(i,v)}, {r(i+1,v)}, {(Subscript[r, u])(i,v)}, {(Subscript[r, u])(i+1,v)}}))+(({{(Subscript[α, 0])(w), (Subscript[α, 1])(w), (Subscript[β, 0])(w), (Subscript[β, 1])(w)}}))(({{r(u,j)}, {r(u,j)}, {(Subscript[r, v])(u,j)}, {(Subscript[r, v])(u,j+1)}}))-
(({{(Subscript[α, 0])(t), (Subscript[α, 1])(t), (Subscript[β, 0])(t), (Subscript[β, 1])(t)}}))(({{r(i,j), r(i,j+1), (Subscript[r, v])(i,j), (Subscript[r, v])(i,j+1)}, {r(i+1,j), r(i+1,j+1), (Subscript[r, v])(i+1,j), (Subscript[r, v])(i+1,j+1)}, {(Subscript[r, u])(i,j), (Subscript[r, u])(i,j+1), (Subscript[r, u v])(i,j), (Subscript[r, u v])(i,j+1)}, {(Subscript[r, u])(i+1,j), (Subscript[r, u])(i+1,j+1), (Subscript[r, u v])(i+1,j), (Subscript[r, u v])(i+1,j+1)}}))(({{(Subscript[α, 0])(w)}, {(Subscript[α, 1])(w)}, {(Subscript[β, 0])(w)}, {(Subscript[β, 1])(w)}})),

t=u-i, w=v-j, 0≤u≤m,0 ≤ v ≤ n. Здесь i=[u] – целая часть u, j=[v] – целая часть v, t=u-i={u} – дробная часть u, w=v-j={v} – дробная часть v, t,w∈[0,1) – местные параметры в данной ячейке сетки.  Области изменения параметров u и v таковы: 0 ≤ u ≤ m, 0 ≤ v ≤ n.

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

(Subscript[r, u])(i,v)=((r(i+1,v)-r(i-1,v))/2),
(Subscript[r, v])(u,j)=((r(u,j+1)-r(u,j-1))/2)

при 1 ≤ u ≤ m-1, 1 ≤ v ≤ n-1 (т.е. для внутренних линий сетки) и по формулам

(Subscript[r, u])(0,v)=2r(1,v)-2r(0,v)-(Subscript[r, u])(1,v),
(Subscript[r, u])(m,v)=2r(m,v)-2r(m-1,v)-(Subscript[r, u])(m-1,v),
(Subscript[r, v])(u,0)=2r(u,1)-2r(u,0)-(Subscript[r, v])(u,1),
(Subscript[r, v])(u,n)=2r(u,n)-2r(u,n-1)-(Subscript[r, v])(u,n-1)
для крайних линий сетки.

Геометрически это означает, что производная (Subscript[r, u])(i,v) в произвольной точке r(i,v) внутренней линии сетки полагается пропорциональной вектору, соединяющему точки r(i-1,v) и r(i+1,v) на соседних с ней линиях сетки, отвечающие тому же значению параметра v, что и данная точка. Аналогичное утверждение справедливо и для производной (Subscript[r, v])(u,j).

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

2.2.4.6. ПОВЕРХНОСТИ ТЕНЗОРНОГО ПРОИЗВЕДЕНИЯ

Рассмотрим поверхность Гордона, построенную по сетке кривых, которые являются кубическими сплайнами Эрмита. Ясно, что такая поверхность целиком определяется узлами сетки и первыми производными в узлах, заданными в обоих параметрических направлениях (поскольку этими величинами целиком определяются кубические сплайны Эрмита, являющиеся кривыми сетки). Таким образом, в данном случае поверхность строится фактически не по сетке кривых, а по набору точек (узлов) с заданными в них производными. Такие поверхности называются ПОВЕРХНОСТЯМИ ТЕНЗОРНОГО ПРОИЗВЕДЕНИЯ

Радиус-вектор поверхности тензорного произведения (в матричных обозначениях) имеет следующий вид:
 r(u,v)=(({{(Subscript[α, 0])(t), (Subscript[α, 1])(t), (Subscript[β, 0])(t), (Subscript[β, 1])(t)}}))(({{r(i,j), r(i,j+1), (Subscript[r, v])(i,j), (Subscript[r, v])(i,j+1)}, {r(i+1,j), r(i+1,j+1), (Subscript[r, v])(i+1,j), (Subscript[r, v])(i+1,j+1)}, {(Subscript[r, u])(i,j), (Subscript[r, u])(i,j+1), (Subscript[r, u v])(i,j), (Subscript[r, u v])(i,j+1)}, {(Subscript[r, u])(i+1,j), (Subscript[r, u])(i+1,j+1), (Subscript[r, u v])(i+1,j), (Subscript[r, u v])(i+1,j+1)}}))(({{(Subscript[α, 0])(w)}, {(Subscript[α, 1])(w)}, {(Subscript[β, 0])(w)}, {(Subscript[β, 1])(w)}})),

t=u-i, w=v-j, 0≤u≤m,0 ≤ v ≤ n. Здесь i=[u] – целая часть u, j=[v] – целая часть v, t=u-i={u} – дробная часть u, w=v-j={v} – дробная часть v, t,w∈[0,1) – местные параметры в данной ячейке сетки.  Области изменения параметров u и v таковы: 0 ≤ u ≤ m, 0 ≤ v ≤ n.

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

2.2.5. ПОВЕРХНОСТИ С ТРЕУГОЛЬНОЙ ПАРАМЕТРИЧЕСКОЙ ОБЛАСТЬЮ

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

2.2.5.1. БИЛИНЕЙНАЯ ТРЕУГОЛЬНАЯ ПОВЕРХНОСТЬ

Пусть даны три точки (Subscript[P, 1]), (Subscript[P, 2]), (Subscript[P, 3]) в (R^3). В качестве области параметров возьмем треугольник A B C в (R^2). Точка параметрического пространства описывается координатами 0≤a,b,c≤1, a+b+c=1. Радиус-вектор треугольной поверхности с вершинами в точках (Subscript[P, 1]), (Subscript[P, 2]), (Subscript[P, 3]) задается формулой r(a,b,c)=a (Subscript[r, 1])+b (Subscript[r, 2])+c (Subscript[r, 3]), где (Subscript[r, i]) – радиус-вектор точки (Subscript[P, i]).

2.2.5.2. ТРЕУГОЛЬНАЯ ПОВЕРХНОСТЬ НА ТРЕХ КРИВЫХ

Пусть даны три попарно пересекающиеся кривые r=(Subscript[r, 1])(u), r=(Subscript[r, 2])(v), r=(Subscript[r, 3])(w), 0≤ u,v,w≤1, пересекающиеся в точках (Subscript[P, 1])=(Subscript[r, 2])(0)=(Subscript[r, 3])(1), (Subscript[P, 2])=(Subscript[r, 1])(1)=(Subscript[r, 3])(0), (Subscript[P, 3])=(Subscript[r, 1])(0)=(Subscript[r, 2])(1). Определим треугольную поверхность формулой:
r(a,b,c)=a((Subscript[r, 3])(1-b)+(Subscript[r, 2])(c)-(Subscript[P, 1]))+b((Subscript[r, 1])(1-c)+(Subscript[r, 3])(a)-(Subscript[P, 2]))+c((Subscript[r, 2])(1-a)+(Subscript[r, 1])(b)-(Subscript[P, 3])), где 0≤a,b,c≤1, a+b+c=1. Легко проверить, что при c=0 формула описывает кривую (Subscript[r, 3])(a) (0 ≤ a ≤ 1), при a=0 – кривую (Subscript[r, 1])(b) (0 ≤ b ≤ 1), при b=0 – кривую (Subscript[r, 2])(c) (0 ≤ c ≤ 1).

ПРОГРАММА 8. Поверхность, соединяющая три поверхности по кривым на их границах. При изменении поверхностей, имеющих вид ``треугольника'', построенная поверхность не изменяется.

DynamicModule[{pe200,pe002,pe101,te,tee,pe020,pe110,pe011,r2,rf2,rf1,u,v,w,s,t,pb020,pb002,pb011,tb,tbb,pb200,pb101,pb110,r1,pj200,pj020,pj110,tj,tjj,pj002,pj101,pj011,r3,rf3,f1perem,f2perem,f3perem,rfel},
pe200={0,10,0};pe002={10,0,0};pe101={10,10,0};
pb020={0,0,20};pb002={10,0,0};pb011={10,0,20};
pj200={0,10,0};pj020={0,0,20};pj110={0,10,20};
pe020[te_]:=If[te≤0.5,{8,8,10}+2 tee ({20,20,0}-{8,8,10})/.tee->te,{20,20,0}+2 (tee-0.5) ({8,8,-10}-{20,20,0})/.tee->te];
pe110[te_]:=If[te≤0.5,{4,9,5}+2 tee ({10,15,0}-{4,9,5})/.tee->te,{10,15,0}+2 (tee-0.5) ({4,9,-5}-{10,15,0})/.tee->te];
pe011[te_]:=If[te≤0.5,{9,4,5}+2 tee ({15,10,0}-{9,4,5})/.tee->te,{15,10,0}+2 (tee-0.5) ({9,4,-5}-{15,10,0})/.tee->te];
r2[te_]:=pe200 u u+pe020[te] v v+pe002 w w+2 pe101 u w+2 pe110[te] u v+2 pe011[te] v w;
rf2[te_,s_,t_]:=r2[te]/.w->1-u-v/.{u->t-s t,v->s t};
f2perem=ParametricPlot3D[Evaluate[r2[te]/.w->1-u-v/.v->0],{u,0,1},Mesh->18,PlotStyle->{Red,AbsoluteThickness[2]}];
pb020={0,0,20};pb002={10,0,0};pb011={10,0,20};
pb200[tb_]:=If[tb≤0.5,{8,6,16}+2 tbb ({20,0,40}-{8,6,16})/.tbb->tb,{20,0,40}+2 (tbb-0.5) ({8,-6,16}-{20,0,40})/.tbb->tb];
pb101[tb_]:=If[tb≤0.5,{9,3,8}+2 tbb ({15,0,20}-{9,3,8})/.tbb->tb,{15,0,20}+2 (tbb-0.5) ({9,-3,8}-{15,0,20})/.tbb->tb];
pb110[tb_]:=If[tb≤0.5,{4,3,18}+2 tbb ({10,0,30}-{4,3,18})/.tbb->tb,{10,0,30}+2 (tbb-0.5) ({4,-3,18}-{10,0,30})/.tbb->tb];
r1[tb_]:=pb200[tb] u u+pb020 v v+pb002 w w+2 pb101[tb] u w+2 pb110[tb] u v+2 pb011 v w;
rf1[tb_,s_,t_]:=r1[tb]/.w->1-u-v/.{u->t-s t,v->s t};
f1perem=ParametricPlot3D[Evaluate[r1[tb]/.w->1-u-v/.u->0],{v,0,1},Mesh->18,PlotStyle->{Red,AbsoluteThickness[2]}];
pj200={0,10,0};pj020={0,0,20};pj110={0,10,20};
pj002[tj_]:=If[tj≤0.5,{6,8,16}+2 tjj ({0,20,40}-{6,8,16})/.tjj->tj,{0,20,40}+2 (tjj-0.5) ({-6,8,16}-{0,20,40})/.tjj->tj];
pj101[tj_]:=If[tj≤0.5,{5,9,8}+2 tjj ({0,15,20}-{5,9,8})/.tjj->tj,{0,15,20}+2 (tjj-0.5) ({-5,9,8}-{0,15,20})/.tjj->tj];
pj011[tj_]:=If[tj≤0.5,{5,4,18}+2 tjj ({0,10,30}-{5,4,18})/.tjj->tj,{0,10,30}+2 (tjj-0.5) ({-5,4,18}-{0,10,30})/.tjj->tj];
r3[tj_]:=pj200 u u+pj020 v v+pj002[tj] w w+2 pj101[tj] u w+2 pj110 u v+2 pj011[tj] v w;
rf3[tj_,s_,t_]:=r3[tj]/.w->1-u-v/.{u->t-s t,v->s t};
f3perem=ParametricPlot3D[Evaluate[r3[tj]/.w->0/.v->1-u],{u,0,1},Mesh->18,PlotStyle->{Red,AbsoluteThickness[2]}];

rfel[tb_,te_,tj_,u_,v_]:=u(rf3[tj,v,1]+rf2[te,0,u+v]-rf2[te,0,1])+v(rf1[tb,1,u+v]+rf3[tj,1-u,1]-rf1[tb,1,1])+(1-u-v)(rf2[te,0,u]+rf1[tb,1,v]-rf1[tb,1,0]); (*три кривые rf1[tb,1,t], rf2[te,0,1-t] и rf3[tj,1-t,1], t на [0,1]*)
Manipulate[Show[{
ParametricPlot3D[Evaluate[rf1[tb,s,t]],{s,0,1},{t,0,1},Mesh->None,PlotRange->All,PlotStyle->{Opacity[0.8],Green}],ParametricPlot3D[Evaluate[rf2[te,s,t]],{s,0,1},{t,0,1},Mesh->None,PlotRange->All,PlotStyle->{Opacity[0.8],Yellow}],ParametricPlot3D[Evaluate[rf3[tj,s,t]],{s,0,1},{t,0,1},Mesh->None,PlotRange->All,PlotStyle->{Opacity[0.8],Blue}],ParametricPlot3D[Evaluate[rfel[tb,te,tj,s,t]],{s,0,1},{t,0,1-s},Mesh->{18,18},PlotRange->All],
f1perem,
f2perem,
f3perem},
ViewPoint->{-0.7,0.4,1.3},ImageSize->{300,200},Boxed->False,Axes->None,PlotRange->{{-6.,20.},{-6.,20.},{-10.,40.}},SphericalRegion->True,ViewAngle->0.7
],
{{tb,0.4,"Зеленый ``треугольник''"},0.,1.,Appearance->"Labeled"},{{te,0.4,"Желтый ``треугольник''"},0.,1.,Appearance->"Labeled"},{{tj,0.4,"Синий ``треугольник''"},0.,1.,Appearance->"Labeled"},SaveDefinitions->True]
]