Основными недостатками кривых Безье являются следующие:
Поэтому при работе с кривыми Безье часто приходиться прибегать к делению кривой Безье, т. е. по сути дела - к замене одной кривой другими. Поэтому естественным образом возникает задача построения такого класса кривых, который бы позволял избежать недостатков 1), 2) и решать соответствующие задачи управления кривой путем изменения управляющих параметров кривой, а не путем ее замены на другие кривые.
Вспомним, что рациональная кривая Безье, построенная по точкам $$p_0, \dots, p_n$$ с весами $$\omega_0, \dots , \omega_n,$$ определялась формулой
$$r_0(t;n)=\frac{\sum_{i=0}^{n}C_n^it^i(1-t)^{n-i} \omega_ip_i}{\sum_{i=0}^n C_n^it^i(1-t)^{n-i} \omega_i}=\frac{\sum_{i=0}^n B_i^n(t) \omega_ip_i}{\sum_{i=0}^n B_i^n(t) \omega_i},$$где $$\{B_i^n(t) = C_n^ii t^i (1 -t)^{n-i}\}$$ - базис Бернштейна на [0,1].
Из формулы (6.1) видно, что глобальное влияние каждой опорной точки $$p_i$$ на всю кривую происходит из-за того, что каждая из функций $$B_i^n$$ не равна нулю на всем интервале (0,1). Если бы вместо функций Бернштейна в формуле (6.1) выступали функции с локальными носителями, существенно меньшими, чем область определения параметра кривой $$t,$$ то удалось бы как локализовать влияние отдельной точки на всю кривую, так и добиться (за счет выбора таких функций) прохождения кривой через некоторую заданную точку.
Поэтому поставим задачу следующим образом. Пусть даны опорные точки $$p_1,\dots., p_n$$ с весами $$\omega_1,\dots ,\omega_n.$$ Определим кривую $$r(t) $$ по формуле
$$r(t)=\frac{\sum_{i=1}^nN_i(t) \omega_ip_i}{\sum_{i=1}^nN_i(t) \omega_i}, t_{min} \le t \le t_{max},$$где $$N_i ( t) $$ - набор кусочно-полиномиальных функций таких, что
для любой гладкой функции $$f(t) $$ существует линейная комбинация $$\sum c_iN_i(t) ,$$ интерполирующая функцию $$f(t) $$ в заданных узлах $$\{t_i\}_{i=1}^n$$
$$f(t_i)=\sum_{j=1}^nc_jN_j(t_i), i=1, \dots, n;$$Последнее условие вводится для того, чтобы в случае совпадения всех опорных точек кривая превращалась бы в ту же точку.
Решение поставленной задачи основывается на понятии
Пусть имеется гладкая на каждом отрезке $$[t_i,t_{i+1}]$$ функция $$f(t) $$ и фиксированы некоторые значения $$t_1, \dots , t_{m+1}$$ переменной $$t$$ (не обязательно различные и не обязательно упорядоченные по возрастанию или убыванию $$t$$ ). Мы также будем рассматривать переменное значение $$t,$$ которое будем считать первым в последовательности. Поскольку оно предшествует значению $$t_1,$$ припишем ему индекс 0 (хотя в большинстве случаев будем обозначать его просто $$t$$ без индекса): $$t_0 = t, t_1, \dots, t_{m+1}.$$ Обозначим $$f_i = f(t_i), f_t = f(t) = f(t_0) = f_0.$$ Определим оператор правого сдвига $$E$$:
$$Et_i = t_{i+1},\\ Ef_i = f_{i+1}\\ (0 \le i \le m).$$Определение 6.1.1. Разделенной разностью функции $$f(t) $$ в узле $$t_i$$ порядка $$k \le m - i + 1$$ называется величина
$$R^{(k)}f_i=\frac{R^{(k-1)}(E-1)f_i}{(E^k-1)t_i}=\frac{R^{(k-1)}f_{i+1}-R^{(k-1)}f_i}{t_{i+k}-t_i},\\ R^{(0)}f_i=f_i.$$При этом если $$t_i = t_{i+1} = \dots = t_{i+k},$$ то в правой части (6.3) появляется отношение вида $$\frac 00,$$ которое по определению полагается равным
$$R^{(k)}f_i=lim_{\delta_1, \dots, \delta_k \to 0, \delta_i > 0} R_{\delta_1, \dots, \delta_k}^{(k)}f_i$$где $$R_{\delta_1, \dots, \delta_k}^{(k)}f_i $$ определяется по возмущенным узлам $$\tilde t_{i+j}=t_i+\delta_1 + \dots, + \delta_j,$$ где $$1 \le j \le k, \delta_j > 0 .$$
Теорема 6.1. Если $$t_i = t_{i+1} = \dots = t_{i+k},$$ то
Пример 6.1.1.
In[1]:=
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 [l/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[[.#l]] , {-20, 20, .01}, Appearance -> "UpArrow" ,
ImageSize -> {200, 12}], Dynamic@ss [[#1]] } /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols -> {ss} ,
SaveDef initions -> True ] ]
Задача 6.1.1. Обратите внимание, что
Задача 6.1.2. Написать программу, которая бы учитывала совпадение узлов.
Пусть $$m \ge 0$$ - целое число. Рассмотрим функцию $$\sigma_m(z,t) = ( max(0, z -t) )^m.$$ При фиксированном $$t$$ функция $$\sigma_m(z,t),$$ рассматриваемая как функция $$z,$$ называется усеченной степенной функцией.
Лемма 6.1. Для производных усеченной степенной функции справедливы формулы
$$\frac{d^k \sigma_m}{dz^k}=\frac{m!}{(m-k)!}\sigma_{m-k}, \frac{d^k \sigma_m}{dt^k}=\frac{(-1)^km!}{(m-k)!}\sigma_{m-k}, k=1, \dots, m$$Обозначим через $$f[t_1,\dots, t_k] $$ разделенную разность порядка $$(k-1) $$ функции $$f(x) ,$$ построенную по узлам $$t_1,\dots, t_k.$$ Фиксируем набор узлов $$t_0 < t_1 \lt; \dots < t_m$$ и рассмотрим функцию $$M(t) =\sigma_{m-1}[t_0, \dots, t_m](t) ,$$ где
Лемма 6.2. Пусть $$m > 0.$$ Тогда $$M(t) = 0,$$ если $$t\le t_0$$ или $$t \ge t_m.$$
Замечание 6.1.1. Из леммы 6.2 вытекает, что $$\sigma_{m-1}[t_0, \dots ,t_m](t) $$ является непрерывной при $$m \ge 2$$ функцией параметра $$t,$$ отличной от нуля лишь на $$(t_0,t_m) .$$
Пример 6.1.2.
In[2]:=f[t_, z_]:=If[t < z, (z-t)^3, 0];
Plot[f[-6, z], {z, -20, 20}]
In[4] :=DynamicModule[ {h, z, t, st, M, G, i, x, y} ,
h[x_, y_] : = If [x < y, (y - x)^3, О] ;
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- l]] = 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[l, Length[st]]},
Grid[{Spacer[3], Subscript[
Style["t", Italic], #-1],
Slider [Dynamic@st[[#1]] , {-15, 13, .01},
Appearance -> "UpArrow" , ImageSize -> {200, 12} ] ,
Dynamic@st|[#l]l} /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols ->{st},
SaveDefinitions -> True] ]
Задача 6.1.3. Написать программу, которая бы учитывала совпадение узлов.
На основе
Определение 6.1.2. Нормированным B-сплайном m-го порядка $$N_{i,m}(t) $$ для неубывающей последовательности узлов $$t_i \le t_{i+1} \le \dots t_{i+m},$$ отсчитываемой от первого узла $$t_i,$$ называется функция
$$N_{i,m}(t)=(t_{i+m}-t_i) \sigma_{m-1}[t_i, \dots, t_{i+m}](t)$$Соответственно, ненормированным B-сплайном $$M_{i,m}(t) $$ называется функция
$$M_{i,m}(t)=\sigma_{m-1}[t_i,\dots, t_{i+m}](t).$$Определение 6.1.3. Нормированным B-сплайном m-го порядка $$N_{m,i}(t) $$ для неубывающей последовательности узлов $$t_{i-m} \le t_{i-m+1} \le \dots \le t_{i-1} \le t_i,$$ отсчитываемой от последнего узла $$t_i,$$ называется функция
$$N_{m,i}(t) = (t_i -t_{i-m})\sigma_{m-1}[t_{i-m}, \dots ,t_i](t).$$Соответственно определяется и ненормированный B-сплайн $$M_{m,i}(t) $$:
$$M_{m,i}(t) =\sigma_{m-1}[t_{i-m}, \dots ,t_i](t).$$Замечание 6.1.2. Если порядок B -сплайна определять не числом узлов минус единица, а степенью усеченной функции (как многочлена t ), то введенные выше B -сплайны будут иметь (по определению) порядок не m, а (m - 1). В этом случае для них применяются другие обозначения: $$N_{-i}^{m-1}, N_{i-}^{m-1}, M_{i-}^{m-1}, M_{-i}^{m-1}.$$ При этом
И те и другие обозначения присутствуют в литературе по B -кривым. Мы будем пользоваться первыми. В Mathematica используются вторые.
Рассмотрим (для простоты) бесконечную в обе стороны последовательность неубывающих узлов $$\{t_i\},$$ где $$t_i \in R, i \in Z, t_i \le t_{i+1}$$ для любого $$i.$$ Фиксируем $$m > 0.$$ Тогда определены множества B -сплайнов: сплайны $$N_{i, m}(t) $$ привязаны к узлу $$i$$ как к первому узлу, сплайны $$N_{m,i}(t) $$ привязаны к узлу $$i$$ как к последнему узлу $$(i \in Z) .$$
Лемма 6.3. Если $$t \in [t_j,t_{j+1}],$$ то среди функций $$N_{i,m}(t) $$ и $$N_{m,i}(t) $$ только m функций в каждом из этих двух классов могут быть отличны от нуля в точке $$t.$$ А именно, отличны от нуля в точке $$t$$ могут быть лишь следующие функции:
$$N_{j-m+1,m}( t ), \dots ,N_{j,m}(t),$$и
$$N_{m,j+1}(t), \dots , N_{m,j+m}(t).$$Лемма 6.4. Для любого $$t \in R$$ имеет место
$$\sum_i N_{i, m}(t)=\sum_iN_{m,i}(t)\equiv 1$$Теорема 6.2 (формула Кокса - де Бура). Справедливы следующие формулы для ненормированных B-сплайнов при любых $$t, m \ge 2, \alpha \in Z$$ и $$\beta = \alpha + m$$:
$$M_{\alpha, m}(t)\frac{(t_{beta}-t)M_{\alpha+1, m-1}(t)+(t-t_{\alpha})M_{\alpha, m-1}(t)}{t_{\beta}-t_{\alpha}},\\ M_{m, \beta}(t)=\frac{(t_{\beta}-t)M_{m-1, \beta}(t)+(t-t_{\alpha})M_{m-1, \beta-1}(t)}{t_{\beta}-t_{\alpha}}\\ M_{\alpha-}^{m-1}(t)=\frac{(t_{\beta}-t)M_{\alpha+1-}^{m-2}(t)+(t-t_{\alpha})M_{\alpha-}^{m-2}(t)}{t_{\beta}-t_{\alpha}}\\ M_{-\beta}^{m-1}(t)=\frac{(t_{\beta}-t)M_{-\beta}^{m-2}(t)+(t-t_{\alpha}M_{-\beta-1}^{m-2}(t)}{t_{\beta}-t_{\alpha}}$$Здесь ненормированные B-сплайны первого порядка равны по определению разделенным разностям функции $$\sigma_0(z) = (z - t)_+^0$$:
$$M_{i,1}(t)=M_{i-}^0=\begin{cases} \frac{1}{t_{i+1}-t_i}, \qquad \mbox{если}t \in [t_i, t_{i+1}),\\ 0, \qquad \mbox{иначе}; \end{cases}\\ M_{1,i}(t)=M_{-i}^0=\begin{cases} \frac{1}{t_i-t_{i-1}}, \qquad \mbox{если } t \in [t_{i-1}, t_i), \\ 0, \qquad \mbox{иначе}. \end{cases}$$Теорема Кокса - де Бура позволяет вычислять B -сплайны рекуррентно с использованием формул (6.4) и (6.5).
Следствие 6.1. Пусть имеется произвольная последовательность (не обязательно неубывающих) узлов $$T = \{_t0, \dots ,t_m\}.$$ Пусть среди них есть 2 неравных друг другу узла $$t_a$$ и $$t_b: t_a \ne t_b,$$ где $$0 \le a,b \le m.$$ Обозначим
$$M^T(t) =\sigma_{m-1}[T] = \sigma_{m-1}[t_0,\dots ,_tm],\\ M^{T/a}(t) = \sigma_{m-2}[T\ a] = \sigma_{m-2}[t_0, \dots ,\hat t_a, \dots ,t_m],\\ M^{T/b}(t) =\sigma_{m-2}[T\b] =\sigma_{m-2}[t_0, \dots ,\hat t_b, \dots ,t_m],\\ p_{ab}(t) =\frac{t_b-t}{t_b-t_a}, q_{ab}(t)=\frac{t-t_a}{t_b-t_a}, (p_{ab}+q_{ab} \evqiv 1) $$где "домиком" отмечены отброшенные узлы. Тогда
$$M^T(t) = p_{ab}(t)M^{T\a}(t) + q_{ab}(t)M^{T\b}(t),\\ M ^{\{t_a, t_b\}} (t)=\begin{cases} \frac{1}{|t_b-t_a|}, \qquad \mbox{если } t \in [t_a \wedge t_b, t_a \vee t_b),\\ 0, \qquad \mbox{иначе}, \end{cases}$$где $$\alpha \vee \beta = max\{\alpha, \beta \}, \alpha \wedge \beta = min\{\alpha, \beta \}.$$
Следствие 6.2. B - сплайны - неотрицательные функции:
Следствие 6.3. Для нормированных B-сплайнов формулы Кокса - де Бура принимают вид $$(\beta = \alpha + m) $$:
$$N_{\alpha,m}(t)=\frac{t_{\beta}-t}{t_{\beta}-t_{\alpha+1}}N_{\alpha+1, m-1}(t)+\frac{t-t_{\alpha}}{t_{\beta-1}-t_{\alpha}}N_{\alpha, m-1}(t),\\ N_{m, \beta}(t)=\frac{t_{\beta}-t}}{t_{\beta}-t_{\alpha+1}}N_{m-1, \beta}(t)+\frac{t-t_{\alpha}}{t_{\beta - 1}-t_{\alpha}}N_{m-1, \beta-1}(t),\\ N_{i-}^0(t)=N_{i,1}(t)=\begin{cases} 1, \qquad \mbox{если } t \in [t_i, t_{i+1})\\ 0, \qquad \mbox{иначе;} \end{cases}\\ n_{-i}^0(t)=N_{1,i}=\begin{cases} 1, \qquad \mbox{если } t \in [t_i, t_{i+1}),\\ 0, \qquad \mbox{иначе}. \end{cases} $$Для сокращения объема вычислений следует сначала вычислить ненормированные B -сплайны нужного порядка по формулам Кокса - де Бура, а затем нормировать их. Если требуется вычислить все ненулевые B -сплайны порядка m для $$t \in [t_i,t_{i+1}]$$ в неубывающей последовательности узлов $$\{t_i\},$$ то вычисления производятся по схеме
Аналогичная схема вычислений применяется для B -сплайнов остальных трех типов.
Следствие 6.4. Для любых $$m \ge 2$$ и $$\alpha \in Z$$ функции $$N_{\alpha,m}(t) $$ и $$N_{m, \beta}(t) ,$$ где $$\beta = \alpha + m,$$ построенные по неубывающей последовательности узлов $$\{t_i\}_{i \in Z},$$ являются непрерывными $$(m - 2) $$ раза дифференцируемыми функциями, обращающимися в ноль вне отрезка $$[t_{\alpha},t_{\beta}].$$ На отрезке $$[t_i,t_{i+1}]$$ они имеют производные до порядка $$(m - 1) $$ включительно. Функции $$N_{\alpha, m}(t) $$ и $$N_{m, \beta}(t) $$ зависят (как от параметров) лишь от значений узлов $$t_{\alpha}, \dots , t_{\beta}.$$
Задача 6.1.4. С помощью встроенной функции в пакет Mathematica $$BSplineBasis[\{m,{t_0, \dots ,t_n\}\}, i, t] ,$$ вычисляющей тот из $$n - m$$ нормированныx B -сплайнов $$N_{i, m+1}(t) $$ порядка $$m + 1$$ (или m при втором способе определения порядка, см. выше), который привязан справа к i -му ( $$0 \le i \le n - (m + 1) $$ ) узлу из расширенного множества неубывающих узлов $$t_0 \le t_1 \le \dots \le t_n$$ (см. определение ниже), наглядно изучить ненормированные и нормированные B -сплайны, привязанные к узлам справа. Обратите внимание, что в Mathematica узлы в расширенном множестве нумеруются начиная с нуля. Соответственно, в Mathematica B -сплайны в базисе нумеруются с нулевого до n - (m + 1) -го.
Замечание 6.1.3. При работе с нижеследующей программой обратите внимание, что при совпадении первых m + 1 узлов значение нулевого (самого первого) нормированного B -сплайна из базиса в нулевом (самом первом) узле становится равно 1, тогда как все остальные B -сплайны базиса обращаются в этом узле в ноль. Это значит, что когда такой базис из нормированных B -сплайнов используется для вычисления B -кривой (аналогично базису Бернштейна для кривых Безье), получившаяся кривая будет начинаться в своей крайней левой опорной точке. Аналогично, при совпадении последних m + 1 узлов расширенного множества, последний сплайн базиса в последнем узле примет значение 1, а остальные будут равны в нем нулю, и, следовательно, правый конец соответствующей кривой будет находиться в ее последней опорной точке. См. ниже теорему 6.3, в формулировке которой порядок B -сплайна m соответствует значению m
+ 1 в обозначениях Mathematica (как уже говорилось, Mathematica использует второй способ определения порядка B -сплайна, а значит, уменьшает его на единицу по сравнению с тем способом, который используем мы).
Задача 6.1.5. Написать программу на пакете Mathematica, вычисляющую B -сплайны $$N_{m,i}(t) ,$$ привязанные к узлам слева.
Пример 6.1.3. Следующие две программы строят графики B -сплайнов, привязанных к узлам справа. В данном примере n = 9, то есть расширенное множество узлов состоит из 10 точек, начиная с нулевой, а m = 3, то есть порядок B -сплайна при принятом нами определении равен 4.
In[5]:=
DynamicModule [ {m = 3, i, tt, knots} ,
For [i = 1, i Sm + 1, i++, tt [i] = 0] ; (*m+l начальных узло равны 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 , 2 m + 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 [N0,m+i [t] , 14] ] ,
Do [If [ss[[i - 1]] > ss[[1]] , ss[[i - 1]] = ss[[i]] , {i, 2, Length [ss] }] ,
Plot[BSplineBasis[{m, knots}, 0, t],
{t, knots[[1]1 , knots[[Length[knots]]] +1} ,
PlotRange -> + {{0, 1.1}, {0, 1. 2} }], PlotLabel -> Style [N0, 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[[#l]] , {0, 1, .01}, Appearance -> " UpArrow" ,
ImageSize -> {200, 12}], Dynamic@ss[el]]} /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols -> {ss} , SaveDefinitions -> True 1 ]]
In|6] :=DynamicModule [{m = 3, n = 9, i, j , tt, knots, a) ,
For [i = 1, i <= m+1, i + +, tt[i] = 0] ;
(*m+l начальных уело равны 0*)
For[j = 1, j im + 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;
Show[If[OrderedQ[ss], Plot[Evaluate[
Table[BSplineBasis[{m, knots}, 1, t] ,
{1, 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 [Hi, 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}, 1, t] ,
{1, 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 [Ni, т+1 [t] , 14]] ] ,
Graphics[Arrow[{{2.5, 1.1},
, {i + 0.1, BSplineBasis[{m, knots}, i, i + 0.1]}}] ] ,
], {{i, 3, " Номер полинома 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[[#l]], {0, 2, .01},
Appearance -> "UpArrow" , ImageSize -> {200 , 12} ] ,
Dynamic@ss|[#l]]} /@r]] ,
ControlPlacement -> Bottom} ,
{{fill, True, "Заполнение"}, {True, False),
ControlPlacement -> Bottom} , SaveDefinitions -> True] ]
Определение 6.1.4. B - кривой порядка m, построенной по вершинам $$p_1, \dots , p_n, n \ge m,$$ с весами $$\omega_1, \dots , \omega_n,$$ называется кривая
определяемая опорными точками $$\{p_i\}_{i=1}^n,$$ весами $$\{\omega_i\}_{i=1}^n$$ и последовательностью узлов $$t_1, \dots ,t_{n+m}.$$
Определение 6.1.5. Множество узлов $$\{t_1, \dots , t_{n+m}\}$$ называется расширенным множеством узлов B -кривой.
Замечание 6.1.4.
Управляющими параметрами B -кривой являются:
m,неубывающая последовательность вещественных чисел $$\{t_1, \dots , _{tn+m}\}.$$
Таким образом, B -кривые обладают дополнительным множеством управляющих параметров по сравнению с рациональными кривыми Безье - расширенным множеством узлов. Именно эти дополнительные параметры и обеспечивают гораздо большую гибкость в управлении B -кривыми, чем кривыми Безье. В современной компьютерной геометрии в основном применяются именно B -кривые.
B -кривые, вообще говоря, не проходят через свои крайние опорные точки. Однако такого прохождения нетрудно добиться с помощью следующего утверждения.Теорема 6.3. Если первые m узлов расширенного множества совпадают, то B -кривая начинается в своей первой опорной точке $$p_1.$$ Аналогично, если последние m узлов расширенного множества совпадают, то B -кривая заканчивается в своей последней опорной точке $$p_n.$$
B -кривые в пакете Mathematica. Mathematica имеет встроенные функции BSplineFunction и BSplineCurve, которые позволяют строить B -кривые по заданному множеству опорных точек пространства $$R^n, n = 2,3,$$ автоматическим или полуавтоматическим подбором параметров. Расширенное множество узлов по умолчанию подбирается автоматически, однако его можно задать самому с помощью параметра SplineKnots. Например, $$SplineKnots \to Automatic$$ означает автоматический подбор расширенного множества узлов, а $$SplineKnots \to \{0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 2\}$$ задает расширенное множество узлов непосредственно.
Так же дело обстоит и с весами опорных точек. Веса можно задавать с помощью SplineWeights. Например, $$SplineWeights \to Automatic$$ означает автоматический выбор весов, а $$SplineWeights \to \{1, 1, 10, 1, 1, 1\}$$ задает конкретные веса для шести опорных точек.
Порядок сплайнов можно задавать с помощью SplineDegree, но при этом должно быть $$SplineKnots \to Automatic.$$
Вычисление порядка происходит так. Число сплайнов в базисе, равное n - m (см. выше), должно равняться числу опорных точек. Поэтому порядок сплайнов может быть вычислен как разность числа узлов и числа опорных точек. Однако для Mathematica эта разность должна быть подобрана так, чтобы она равнялась 4, иначе B -кривая может быть построена в Mathematica неправильно.
Чтобы воспользоваться возможностью полноценного управления B-кривой, в том числе и весов, нужно запрограммировать B -кривую самостоятельно с помощью встроенной функции BSplineBasis.
Задача 6.1.6. С помощью встроенных функций Mathematica SplineFunction и BSplineCurve построить B -кривую по заданным опорным точкам с автоматическим выбором узлов на плоскости и в трехмерном пространстве.
Пример 6.1.4. B -кривая на плоскости с автоматическим выбором узлов, построенная двумя способами: с помощью функций BSplineFunction и BSplineCurve соответственно. На втором графике для сравнения нарисована также простая кривая Безье, построенная по тем же опорным точкам.
In[7] : = DynamicModule [ {ptsO, f, g, n, ptsl, pts2 , t} ,
n = Length [ptsO] ;
Manipulate[f = BSplineFunction[ptsl] ;
Show[{
Graphics[{Green, Line[ptsl], Black,
Text["B-кривая, построенная через BSplineFunction", {12, 10-0}],
Text[ToString[# - 1] , ptsl [ [#] ] +{0.7, 0}] /@ Range [n] } ,
PlotRange -> {{-4, 20}, {-4, 10}}, Axes -> True] ,
ParametricPlot[f[t], {t, 0, 1}]}, ImageSize -> 300] ,
{{ptsl, 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}})]
In[8] :=
DynamicModule [ {ptsO , f, g, n, ptsl, pts2 , t} ,
n = Length [ptsO] ;
Manipulate[g = BezierFunction[pts2];
Show [ {
Graphics [ {Red, BSplineCurve [pts2 , SplineKnots -> Automatic] , Black,
Text["В-кривая, построенная через BSplineCurve (красная)", {10, 10.0}],
Техt["Криая Беэье (голубая)", {10, 9.0}],
Text [ToSt ring] # - 1] , pts2 [ [#] ] + {0.7, 0}] /@ Range [n] } ,
PlotRange -> {{-4, 20} , {-4, 10}} , Axes -> True] ,
ParametricPlot[g[t] , {t, 0, 1}, ElotStyle -> {Blue}]
}, ImageSize -> 300] , { {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}})]
Пример 6.1.5. B -кривая в пространстве с автоматическим выбором узлов в сравнении с кривой Безье, построенной по тем же опорным точкам. Обратите внимание, насколько более послушно ведет себя B -кривая при управлении опорными точками, чем кривая Безье. Кривая Безье по-другому - более глобально - усредняет влияние опорных точек на кривую, что делает ее локальное управление за счет смещения отдельных опорных точек гораздо более вязким, чем в случае B -кривой.
In[9]:=
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["t1", Italic] , #- l] ,
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 ["t2 и t3", Italic], #-l],
Slider2D [Dynamic® {pts[[1]] [ [2] ] , pts[[1]] [ [3] ] } ,
{{-15, -15}, {15, 15}}, Appearance -> "Left",
ImageSize -> {50, 30}],
Dynamic@ {pts [[#1]] [ [2] ] , pts [[#2]] [[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.1.6. B -кривая на плоскости с автоматическим выбором порядка, но с управляемыми мышью опорными точками и узлами (узлы должны быть неубывающими).
In[10]:=
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[l, m] } ,
Grid[
{Spacer[3] , Subscript [Style ["t" , Italic] , # - 1] ,
Slider [Dynamic@tt[[#l]] , {0, 15, .01}, Appearance -> "UpArrow" ,
ImageSize -> {200, 12}], Dynamic@tt[[i]]} /@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])]
Пример 6.1.7. B -кривая на плоскости с автоматическим выбором порядка и узлов, но с управляемыми мышью опорными точками и весами, и простая кривая Безье, управляемая теми же опорными точками:
In[11]: =
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 [To St ring [# - 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, "Степень В-кривой" } ,
Dynamic@RadioButtonBar [Dynamic [d] , Range[0, 7]] } ,
{{curve, 1, "Кривая"} , {1 -> "В-кривая" , 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[l, {i, l,n}])]
Задача 6.1.7. Самостоятельно запрограммировать B -кривую с помощью встроенной функции BSplineBasis, так, чтобы можно было воспользоваться возможностью полноценного управления B -кривой. Обеспечить возможность произвольного допустимого порядка B -сплайнов.
Пусть даны опорные точки $$p_1, \dots , p_n$$ и их веса $$\omega_1, \dots , \omega_n,$$ а также расширенное множество неубывающих узлов $$\{t_1, \dots ,t_{n+m} : t_1 \le \dots \le t_{n+m}\}.$$ Опишем алгоритм вычисления B -кривой
Фиксируем t. Положим $$t_0 = -\infty, t_{n+m+1} = + \infty.$$ Тогда существует единственное $$i_0$$ такое, что $$t \in [t_{i_0},t_{i_0+1}),$$ где $$0 \le i_0 \le n + m.$$
Для определенного выше индекса $$i_0,$$ если $$1 \le i_0 \le n + m,$$ вычисляем единственное ненулевое значение ненормированного B -сплайна первого уровня (m=1):
Напомним, что $$i_0$$ было выбрано на первом шаге, так чтобы $$t \in [t_{i_0}, t_{i_0+1}).$$ При этом, в силу условия $$1 \le i_0 \le n + m - 1,$$ знаменатель $$t_{i_0+1} - t_{i_0} < \infty,$$ а в силу условия $$t \in [t_{i_0}, t_{i_0+1})$$ имеем $$t_{i_0+1} - t_{i_0} > 0.$$ Следовательно, $$0 < M_{i_0, 1}(t) > \infty.$$ Для индекса $$i_0,$$ находящегося вне множества $$\{1, \dots , n+m-1\},$$ значение $$M_{i_0,1}$$ не вычисляется.
С помощью соотношений Кокса - де Бура вычисляем все отличные от нуля в точке $$t$$ ненормированные сплайны $$m$$ -го порядка $$M_{j,m}(t) $$ при $$1 \le j \le n$$:
$$M_{j, m}(t)= \frac{(t_{j+m}-t)M_{j+1, m-1}(t)+(t-t_j)M_{j, m-1}(t)}{t_{j+m}-t_j}$$
В частности,
$$M_{i_0, s}(t)=\frac{(t-t_{i_0})^{s-1}}{(t_{i_0+1}-t_{i_0}) \dots (t_{i_0+s}-t_{i_0})},\\ N_{i_o, s}(t)=\frac{(t-t_{i_0})^{s-1}}{(t_{i_0+1}-t_{i_0}) \dots (t_{i_0+s-1}-t_{i_0})},$$где $$1 \le s \le m, t \in [t_{i_0}, t_{i_0+1}).$$ Положив в (6.6) $$i_0 = 1,$$ получим
$$N_{1,s}(t)=\frac{(t-t_1)^{s-1}}{(t_2-t_1) \dots (t_s-t_1)},\\ 1 \le s \le m,\\ t \in [t_1, t_2).$$Лемма 6.5. Имеет место соотношение:
$$N_{1, m}(t)= \begin{cases} 1, \qquad \mbox{если } t_1= \dots = t_m,\\ 0, \qquad \mbox{иначе}. \end{cases}$$Вычисляем нормированные сплайны $$N_{j, m}(t) $$ для каждого $$1 \le j \le n$$ по формуле $$N_{j, m}(t) = (t_{j+m} - t_j)M_{j, m} (t) $$ (при этом $$1 \le j,j + m \le n + m$$ ).
Окончательно вычисляем $$r(t) $$ по формуле
$$r(t)=\frac{\sum_{i=1}^nN_{i, m}(t) \omega_i p_i}{\sum_{i=1}^n N_{i, m}(t) \omega_i}$$Теорема 6.4. r(t) B -кривой:
может быть вычислен с помощью следующего алгоритма:
Для заданного $$t \in [t_{i_0},t_{i_0+1}) (m \le i_0 \le n) $$ вычисляем величины $$r_i^{(k)}, \omega_i^{(k)} $$ при $$k = 1, \dots , m - 1, i = i_0 - m + k + 1, \dots , i_0$$ по рекуррентным формулам:
$$r_i^{(k)}=\frac{t_{i+m-k}-t}{t_{i+m-k}-t_i}r_{i-1}^{(k-1)}+\frac{t-t_i}{t_{i+m-k}-t_i}r_i^{(k-1)},\\ \omega_i^{(k)}=\frac{t_{i+m-k}-t}{t_{i+mk}-t_i}\omega_{i-1}^{(k-1)}+\frac{t-t_i}{t_{i+m-k}-t_i} \omega_i^{(k-1)},\\ r_o^{(0)}= \omega_i p_i, \omega_i^{(0)}=\omega_i, i=i_0-m+1, \dots, i_0 $$Вычисляем r(t) по формуле
Данный алгоритм иллюстрируется следующей диаграммой ( $$t \in [t_{i_0}, t_{i_0+1})$$ ):
$$\begin{matrix} r_{i_0-m+1}^{(0)} r_{i_0-m+2}^{(0)} \dots r_{i_0-1}^{(0)} r_{i_0}^{(0)}\\ r_{i_0-m+2}^{(1)} \dots r_{i_0-1}^{(1)}r_{i_0}^{(1)}\\ \ddots \vdots \vdots \\ r_{i_0-1}^{(m-2)} r_{i_0}^{(m-2)}\\ r_{i_0}^{(m-1)} \end{matrix}$$Определение 6.2.1. Пусть дано (m +1) * (n+1) точек в пространстве $$R^3,$$ образующих прямоугольную матрицу (сетку Безье):
Поверхностью Безье порядка m * n, соответствующей сетке $$P_{mn},$$ называется поверхность
где E - оператор сдвига вперед по первому индексу, F - оператор сдвига вперед по второму индексу:
Операторы E и F очевидно коммутируют друг с другом: $$EFp_{ij} = FEp_{ij} = p_{(i+1)(j+1)} .$$ Поэтому $$(1-v + Fv)^n(1 -u + Eu)^mp_{00} = (1 - u + Eu )^m(1 -v + Fv)^np_{00},$$ т. е. формула (6.7) непротиворечива.
Определение 6.2.2. Рациональная поверхность Безье, построенная по точкам $$p_{ij}$$ с весами $$\omega_{ij}$$ определяется следующим образом:
$$S_{00}(u, v; 1,1)=\frac{(1-u+Eu)(1-v+Fv) \omega_{00}p_{00}}{(1-u+Eu)(1-v+Fv) \omega_{00}}$$Формула (6.8) означает, что мы строим поверхность Безье в $$R_+^{k+1}$$ по точкам $$(\omega_{ij}; р_{ij} ,\omega_{ij}) \in R_+^{k+1},$$ а затем применяем преобразование $$\varphi R_+^{k+1} \to R^ k , \varphi (x, \omega) =\frac{x}{\omega}.$$
Задача 6.2.1. Наглядно изучить влияние опорных точек и их весов на рациональную поверхность Безье, меняя исходные данные pts (опорные точки) и $$\omega$$ (их веса) в нижеследующей программе. В Mathematica применять функцию BezierFunction к матрицам из скаляров.
Пример 6.2.1. Рациональная поверхность Безье с возможностью непосредственного управления весами с помощью движков.
In[12]: =
DynamicModule [ {pts, a, w0, pw, w, g, f, w6, w7, wl0, wll, 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}, {wl0, 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]]}],
ParametrioPlot3D[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},
{{wl0, 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])]
Поверхность Безье можно получить следующим образом:
Строим (n + 1) кривую Безье по столбцам матрицы $$P_{mn}$$:
u.При этом соответствующий оператор перехода от опорной точки, взятой на кривой $$r_j(u) ,$$ к следующей опорной точке, взятой на кривой $$r_{j+1}(u) ,$$ в терминах узлов сетки pij описывается оператором сдвига F по второму индексу. Следовательно,
$$S(u,v) = (1-v + Fv)^nr_0(u) = (1-v + Fv)^n (1-u + Eu )^m p_{00} = S_{00}(u,v;m,n),$$где m и n - количество шагов по первому и второму индексу соответственно, начиная от угловой точки сетки Безье, 00 - индекс этой угловой точки, с которой мы начинаем образовывать все остальные узлы сетки операторами сдвига E и F. Имеем
Здесь через $$S_{ij}$$ обозначены четырехугольные листы Безье, построенные из точек p ij, взятых в качестве угловых точек, аналогично формуле (6.9).
Аналогично формуле (6.10), любая точка на поверхности Безье степени m*n может быть представлена в виде точки четырехугольного листа Безье с углами в точках соответственно подобранных поверхностей Безье порядка (m-1)*(n-1):
Поверхность Безье можно представить в следующим виде:
$$S_{00}(u,v;m,n) = F(u,m)PF^{\top}(v,n),$$где $$P = (p_{ij})$$ и
$$F(u,m) = (f_0(u,m), \dots ,f_m(u,m) ) ,\\ F(v,n)= ( f_0(v,n), \dots ,f_n(v,n) ) ,\\ f_i ( \xi,k) = C_k^i(1- \xi)^{k-i} \xi^i.$$Аналогично обычной поверхности Безье, рациональная поверхность Безье $$S_{00}(u, u; 1,1) $$ может быть получена с помощью следующего алгоритма:
Полученная точка и будет точкой с параметрами $$(u_0,v_0) $$ на рациональной поверхности Безье, построенной по $$(p_{ij},\omega_{ij})_{i,j=0}^1.$$
Деление поверхностей Безье целиком основано на делении кривых Безье. Как мы уже знаем, каждая из кривых Безье (обычных или рациональных), построенных по столбцам сетки опорных точек $$P_{mn},$$ может быть разделена в параметрической точке $$u = u^* \in [0,1] $$ без изменения ее формы. При таком делении получаются две поверхности Безье (составляющие вкупе исходную поверхность) и соответствующие им две матрицы опорных точек Безье, каждая из которых имеет тот же порядок, что и исходная матрица: $$P_{тп}^a$$ и $$P_{тп}^b.$$ Затем, каждую из получившихся двух поверхностей Безье надо поделить в параметрической точке $$v = v^*$$ (с помощью деления кривых Безье, построенных по строкам соответствующих матриц). В итоге получится четыре поверхности Безье и соответствующие им четыре матрицы опорных точек: $$P_{mn}^{ac}, P_{mn}^{ad}, P_{mn}^{bc}$$ и $$P_{mn}^{bd}.$$ Эти четыре поверхности Безье будут лежать на исходной поверхности Безье, делить ее на четыр е части и в совокупности составлять исходную поверхность.
Новые сети опорных точек $$P_{mn}^{ac}, P_{mn}^{ad}, P_{m n}^{bc}$$ и $$P_{mn}^{bd}$$ получаются из исходной сети $$P_{mn}$$ по следующим формулам, аналогичным формулам (5.26):
$$p_{ij}^{ac} = (1-u^*+Eu^*)^i(1-v^*+Fv^*)^jp_{00}\\ p_{(m-i)j}^{bc} = ( (1 - u^*)E^{-1} + u^* )^i (1 -v^* + Fv^*)^jp_{m0},\\ p_{i(n-j)}^{ad} = (1-u^*+ Eu^*)^i ( (1 - v^*)F^{-1} + v^* )^j p_{0n},\\ p_{(m -i)(n-j)}^{bd} = ( (1 - u^*)E^{-1} + u^* )^i ( (1 - v^*)F^{-1} + v^* )^j p_{mn}.$$Добавление новых узлов в сетку Безье без изменения самой поверхности происходит по формулам, полностью аналогичным формулам добавления опорных точек в кривую Безье. Если дана сетка Безье $$P_{mn} = (p_{ij})_{i,j}^{m,n}=0,$$ то увеличенная на одну строку и на один столбец сетка $$Q_{(m+1)(n+1)} $$ для той же самой поверхности Безье строится по формулам
$$q_{ij}=\frac{1}{m+1}\frac{1}{n+1}((m+1-i)(1 - \delta_{(m+1)i})+i(1-\delta_{0i})E^{-1})\times\\ \times((n+1-j)(1-\delta_{(n+1)j})+j(1-\delta_{0j}F^{-1})p_{ij}.$$При этом углы новой сетки совпадут с углами старой: $$q_{00} = p_{00}, q_{(m+1)0} = p_{m0}, q_{0(n+1)} = p_{0n}, q_{(m+1)(n+1)} = p_{mn}.$$
Прежде всего заметим, что с помощью деления поверхности Безье любая наперед заданная точка этой поверхности может быть представлена как угловая точка новой поверхности Безье, являющейся частью исходной, а потому имеющей в каждой своей точке те же самые дифференциально-геометрические свойства, что и исходная поверхность Безье. Следовательно, для вычисления дифференциально-геометрических характеристик поверхности Безье в произвольной точке достаточно уметь вычислять эти характеристики лишь для угловой точки поверхности Безье с нулевыми значениями параметрических координат.
Таким образом, достаточны следующие формулы.
Касательные векторы к поверхности Безье в угловой точке $$p_{00}$$ могут быть вычислены по формулам $$T_u(0) = m(p_{10} - p_{00})$$ и $$T_v(0) = n(p_{01} - p_{00}).$$
Нормаль к поверхности в угловой точке дается формулой $$N =\frac{T_u(0) \times T_v(0)}{|T_u(0) \times T_v(0)|} $$
Вторые производные
Вычислив, согласно приведенным формулам, первые и вторые производные в угловой точке, а также нормаль к поверхности в этой точке, мы сможем вычислить гауссову и среднюю кривизны, главные направления и главные кривизны в угловой точке. Все эти величины будут выражаться через узлы сетки: $$p_{00}, p_{01}, p_{02}, p_{10}, p_{20}.$$
Для произвольной точки $$S_{00}(u^*,v^*;m,n) $$ поступаем так:
Пример 6.2.2. Чайник, нарисованный с помощью поверхности Безье. Чайник состоит из заплат:
In[13]:=
DynamicModule[{potPoints, рр2, potPatches, t, i, n, Bernstein,
BezierSurface}, potPoints = {{1.4, 0, 2.4},
(* Оставшиеся 304 точки *)
{1.425, -0.798, 0}}; (* это точки каркаса чайника, их 306 *)
рр2 = With[{avg = Mean [potPoints]}, Map[(j" - avg) , potPoints]];
(*усреднили все точки и составили из них список*)
potPatches = {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16},
(* оставшиеся тридцать заплат *)
{270, 270, 270, 270, 300, 305, 306, 279, 297, 303, 304, 275, 294,
301, 302, 271}}; ("чайник сосотоит из заплат, их 32 штуки*)
Bernstein [i_Tnteger, 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_, y_, 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, "Разбиение"}, 0, 1},
SaveDefinitions -> True]]
Пусть задана матрица опорных точек $$(p_{ij} ), 1 \le i \le n, 1 \le j \le m,$$ с весами ( $$\omega_{ij}$$ ). Фиксируем числа $$p$$ и $$q$$ - порядки B -сплайнов в двух координатных направлениях u и v соответственно. Тогда B -поверхность порядка (p,q), построенная по матрице опорных точек $$(p_{ij}) $$ с весами $$(\omega_{ij}) $$ и по двум расширенным множествам узлов в направлениях u и v соответственно, определяется формулой
Если ввести новые обозначения для входящих в эту формулу произведений B -сплайнов
то формула (6.11) перепишется в виде
$$r(u,v)=\frac{\sum_{i=1}^n \sum_{j=1}^m N_{ij,qp}(v,u) \omega_{ij}p_{ij}}{\sum_{i=1}^n \sum_{j=1}^m N_{ij,qp}(v,u) \omega_{ij}}$$Если расширенные множества узлов в u и v -направлениях имеют вид
то область изменения параметров (u,v) на B -поверхности (6.12) будет следующей:
Задача 6.2.2. Наглядно изучить влияние опорных точек и расширенного множества узлов на B-поверхность, порождаемую пакетом Mathematica с помощью следующей программы. Поменять в программе автоматический выбор расширенных множеств узлов на прямое задание этих узлов. Изучить влияние расширенных множеств узлов на B -поверхность.
Пример 6.2.3. B -поверхность. Для обеих поверхностей задан одинаковый остов. В первом случае поверхность "замыкается" только по оси Ox, во втором случаем - по обеим осям.
In[14] : = 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 = BSplineEunction [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}]]
Пример 6.2.4. Моделирование кузова автомобиля с помощью В -поверхности.
In[19]: =
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, 140}}], 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}}
(* Другие опорные точки *)
{{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}}
})]
Пример 6.2.5. Интерполяция заданного рисунка с помощью поверхности Безье и B -поверхности.
In[20]: =
DynamicModule {f, colfun, prevn =5, prevapprox =True, start, locators, grid,
order, approx, images},
Manipulate {* Создали список, который задает рисунок,
и для этого списка построили В-функцию, т.е. мы "синтерполировали" рисунок*)
colfun = BSplineFunction[ImageDatafimages[[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"1}] ,
f = Interpolation[
Flatten[Table[{{N[1 / (n-1)] i, N[1/ (n-1)] j } , pts [ [nit] + 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, {}]},
ImageSise -> 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 -> "В-поверхность" } ,
ControlType -> RadioButton} ,
{{order, 3, "Порядок В-поверхности"},
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, l/(n-l)}, {j, 0, 1, l/(n-l)}], 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 = {
};) ]]
В ряде приложений (например, в проектировании корпусов и крыльев самолетов и кораблей) требуется строить гладкие поверхности, затягивающие заданные остовы из кривых - криволинейные сетки, криволинейные направляющие (например, несущий остов корпуса корабля), различные граничные контуры. Во многих случаях важное значение имеет также и вычислительная сторона вопроса - сложность вычисления
В компьютерной геометрии существует несколько различных подходов к решению подобных задач. Ниже будет рассмотрены некоторые из этих подходов, ставшие классическими.
Линейчатая поверхность - геометрическое место точек отрезков прямых, соединяющих соответствующие друг другу точки двух заданных кривых. Пусть $$\rho_1(t), t_0 \le t \le t_1,$$ и $$\rho_2(w), w_0 \le w \le w_1,$$- две заданные кривые. Тогда определяемая ими линейчатая поверхность имеет вид
$$r(u,v) = (1-v)\rho_1(t) + v \rho_2(w),$$где $$t = t_0(1-u)+t_1u, w = w_0(1 - u) + w_1u, 0 \le u, v \le 1.$$
Если обе кривые $$\rho_1$$ и $$\rho_2$$ замкнуты, то линейчатая поверхность замкнута по u. По параметру v линейчатая поверхность всегда незамкнута.
Если обе кривые $$\rho_1$$ и $$\rho_2$$ - отрезки прямых, то линейчатая поверхность называется билинейной и определяется только своими четырьмя угловыми точками $$p_1, p_2, p_3, p_4$$:
$$r(u, v) = (1 - v) ( (1 - u)p_1 + up_2 ) + v ( (1 - u)p_3 + u p_4 ) =\\ = (1 - u)(1 - v)p_1 + u(1 - v)p_2 + (1 - u)vp_3 + uvp_4,\\ 0 \le u, v \le 1.$$Секториальная поверхность - частный случай линейчатой поверхности, когда одна из кривых $$\rho_1$$ и $$\rho_2$$ вырождается в точку:
$$r(u,v) = (1-v) \rho (u)+vp,$$где $$\rho ( u) $$ - кривая, p - точка, $$u_1 \le u \le u_2, 0 \le v \le 1.$$
Секториальная билинейная поверхность называется треугольной поверхностью.
Линейные поверхности Кунса. Вычтем и добавим к уравнению линейчатой поверхности уравнение билинейной поверхности, построенной по четырем угловым точкам (см. выше). Уравнение линейчатой поверхности примет вид
$$r(u,v) = (1-v)\rho_1(u)+v \rho_2(u) = (1-v) \rho_1(u)+v \rho_2(u)+(1-u) ( (1-v)p_1+vp_3 ) +\\ + u ( (1 - v)p_2 + vp_4 ) - (1 - u)(1 - v)p_1 - u(1 - v)p_2 - (1 - u)vp_3 - uvp_4 =\\ = (1 - v) \rho_1(u) + v \rho_2(u) + (1 - u) \rho_3(v) + u \rho_4(v) - (1 - u)(1 - v)p_1-\\ - u(1 - v)p_2 - (1 - u)vp_3 - uvp_4,\\ 0 \le u, v \le 1,$$где $$\rho_1 (u) $$ и $$\rho_2 (u) $$ - кривые, по которым строится данная линейчатая поверхность, причем считается, что на них обеих задан один и тот же параметр $$0\leu\le1,$$и
$$\rho_3(v) = (1-v)p_1+vp_3,\\ \rho_4(v) = (1 - v)p_2 + vp_4$$- отрезки прямых.
Подставим теперь в уравнение (6.13) вместо отрезков прямых (6.14) уравнения двух произвольных кривых $$\rho_3(v) $$ и $$\rho_4(v) ,$$ соединяющих соответственно пары точек $$(p_1, p_3) $$ и $$(p_2, p_4) .$$ Тогда новая поверхность, задаваемая уравнением (6.13), будет определяться четырьмя кривыми $$\rho_1, \rho_2, \rho_3, \rho_4,$$ образующими криволинейный четырехугольник, натянутый на угловые точки $$p_1, p_2, p_3, p_4.$$ При этом все четыре кривые $$\rho_1, \rho_2, \rho_3, \rho_4$$ войдут в уравнение (6.13) симметрично. Получится линейная поверхность Кунса.
Определение 6.3.1. Пусть заданы четыре кривые, образующие криволинейный четырехугольник с углами в точках $$p_1, p_2, p_3, p_4$$:
$$\begin{martix} \gamma_a=r(u,0),0 \le u \le 1,r(0,0)=p_1,r(1,0)=p_2,\\ \gamma_b=r(0,v) 0 \le v \le 1, r(0,0)=p_1 r(1,0)=p_3),\\ \gamma_c(1,v) 0 \le v 1 r(0,0)=p_1 r(0,1)=p_3,\\ \gamma_d=r(u,1),0 \le u \le 1, r(0,1)=p_3, r(1,1)=p_4 \end{matrix} $$Обозначим
$$\alpha_0(t) = 1-t,\\ \alpha_1(t)=t.$$Функции $$\alpha_0(t), \alpha_1(t) $$ будем называть функциями смещения линейной поверхности Кунса. Сама линейная поверхность Кунса определяется следующей формулой:
$$r(u, v) = \alpha_0(v)r(u, 0) + \alpha_0(u)r(0, v) + \alpha_1(v)r(u, 1)+\\ + \alpha_1(u)r(1,v) -\alpha_0(u) \alpha_0(v)r(0,0) -\alpha_1(u)\alpha_0(v)r(1,0)-\\ - \alpha_0(u)\alpha_1(v)r(0,1) -\alpha_1(u)\alpha_1(v)r(1, 1).$$Замечание 6.3.1. Линейная поверхность Кунса не является, вообще говоря, линейчатой поверхностью. Название "линейная" обусловлено тем, что в ее определении используются линейные функции смещения $$\alpha_0$$ и $$\alpha_1.$$
Замечание 6.3.2. Параметрическая область поверхности Кунса может иметь вид прямоугольника, не обязательно совпадающего с квадратом $$[0,1] \times [0,1] .$$
Матричный вид уравнения поверхности Кунса. Преобразуем формулу (6.15) к другому виду. Имеем в матричных обозначениях
$$r(u,v) = ( \alpha_0(u) \alpha_1(u) ) \begin{pmatrix} r(0,v)\\ r(1,v) \end{pmatrix} + ( \alpha_0(v) \alpha_1(v) ) \begin{pmatrix} r(u,0)\\ r(u,1) \end{pmatrix} -\\ -( \alpha_0(u) \alpha_1(u)) \begin{pmatrix} r(0,0) r(0,1) r(1,0) r(1,1) \end{pmatrix} \begin {pmatrix} \alpha_0(v)\\ \alpha_1(v) \end{pmatrix} =\\ = \sum_{i=0}^1 \alpha_i(u)r(i, v) + \sum_{j=0}^1 \alpha_j(v)r(u,j) - \sum_{i=0}^1 \sum_{j-0}^1 \alpha_i(u)\alpha_j(v)r(i, j).$$Пример 6.3.1. Зададим с помощью поверхности Кунса область D плоскости, ограниченную эллипсом $$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$ или $$\begin{cases}
x=a\ cos\ \varphi\\
y=b\ cos\ \varphi,
\end{cases}$$ где $$0 \le \varphi < 2 \pi,$$ с полуосями $$a > 0, b > 0.$$ Представим границу области D в виде четырех кривых - дуг эллипса, стыкующихся в точках $$p_1, p_2, p_3, p_4.$$ Соответствующая линейная поверхность Кунса имеет вид
Преимущество задания плоских областей в виде линейных поверхностей Кунса состоит в том, что у поверхности Кунса параметрическая область -всегда прямоугольник.
Пример 6.3.2. Поверхность Кунса, построенная по кривым, лежащим на заданных поверхностях. Кривые раскрашены красным, зеленым, желтым и оранжевым.
In[21]: =
DynamicModule[{rl = 30 . , r2 = 30. , r3 = 30 . , d = 50 . , u , v, ml, m2 , m3, m4 , fl,
f2, f3, f4, rfel, α, MV, MU, laml, lam2, rfk}, rl = 30. ; r2 = 30. ; r3 = 30. ;
d = 50. ; laml = 2. ; lam2 = 1. ;
f1 [ml_] : = {rl - rl *Cos[u * Pi] , laml * 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, m1* (Sin[v* 3* Pi] )^2} :
f4 [m4_] :={1.6*r1+0.4*rl*u+(v^4-2*v^3+v^2)*100*m4*u, d*v, 0} ;
Manipulate[
Show[{ParametricPlot3D[Evaluate@fl [ml] , {u, 0, 1}, {v, -1, 0}, Mesh -> {24, 24},
PlotRange -> All] ,
ParametricPlot3D[Evaluate@{fl [ml] /. 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 -> l}, {u, 0, 1}, Mesh -> {24, 24},
PlotRange -> All, PlotStyle -> {Green} ] ,
ParametricPlot3D[Evaiuate@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 -> l}, {v, 0, 1}, Mesh -> {24, 24},
PlotRange -> All, PlotStyle -> {Orange} ] ,
ParametricPlot3D [Evaluate@rfel [ml, 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}}
], {{ml, 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 : -> (
rl = 30 . ; r2 = 30. ; r3 = 30. ; d = 50. ; laml = 2. , lam2 = 1. ;
fl [ml_] : = {rl - rl *Cos [u * Pi] , laml *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*rl+0.4*rl*u+(v^4-2*v^3+v^2)*100*m4*u, d*v, 0} ;
\alpha [t_] := t; MV[v_] : = {1 - \alpha [v] , \alpha [v] } ; MU[u_]: = {1 - \alpha [u] , \alpha [u] } ;
(* наши функции а *)
Table [rfk [ml, m2, m3, m4] [i] , {i, 3}] ;
Do[rfk[ml_, m2_, тЗ_, т4_] [i] =
MU[u].{f3[m3][[i]] /. u ->0, f4[m4][[i]] /. u ->l} +
MV[v] . {fl [ml] [ [i] ] /. v->0, f2[m2][[i]] /. V->1} -
MU[u] . {{f3[m3] [[i]] /. {u ->0, v-> 0} , f3[m3] [[i]] /. {u -> 0, v->l}},
{f4 [m4] [ [i] ] /. {u -> 1 , v-> 0} , f4 [m4] [ [i] ] /. {u-> 1, v -> 1} } } .MV[v] //
Together;
, {i, 3}] ; rfel[ml_, m2_, тЗ_, m4_] =
{rfk [ml, m2, m3, m4] [1] , rfk [ml, m2, m3, m4] [2] , rfk [ml, m2, m3, m4] [3] };) ]]
Обобщенные поверхности Кунса. При стыковке линейных поверхностей Кунса по граничным кривым производные в направлении, ортогональном границе, претерпевают разрыв. Чтобы гладко склеивать поверхности Кунса, необходимо строить их таким образом, чтобы они имели на границе заданные (согласованные) производные до k -го порядка включительно. Производные в касательном направлении определяются уравнениями самих граничных кривых, а производные в трансверсальном к границе направлении необходимо задавать дополнительно.
Определение 6.3.2. Обобщенной поверхностью Кунса n-го порядка называется поверхность
$$r(u,v)=\sum_{i=0}^1 \sum_{k=0}^n \alpha_{k,i}(u)r_u^{(k)}(I,v)+\sum_{j=0}^1 \sum_{k=0}^n \alpha_{k,j}(v)r_v^{(k)}(u,j)-\\ -\sum_{i=0}^1 \sum_{j=0}^1 \sum_{k=0}^n \sum_{l=0}^n \alpha_{k,i}(u) \alpha_{l,j}(v)r_{uv}^{(k,l)}(I,j),\\ 0 \le u,v \le 1$$задаваемая четырьмя граничными кривыми $$r(u,0), r(0,v), r(1,v), r(u, 1) ,$$ стыкующимися в точках $$p_1 = r(0, 0), p_2 = r(1,0), p_3 = r(0,1) $$ и $$p_4 = r(1,1) ,$$ а также чистыми частными производными в трансверсальном направлении на граничных кривых до $$n$$ -го порядка включительно и смешанными частными производными в угловых точках $$p_1, p_2, p_3, p_4$$ до $$2n$$ -го порядка включительно:
$$\begin{matrix} r_u^{(k)}(I,v)=\frac{\partial^k r(u,v)}{ \partial u^k} \left |_{u=i}, i=0,1, k=1, \dots, n\\ r_v^{(k)}(u,j)=\frac{\partial^kr(u,v)}{ \partial v^k} \left |_{v=j}, j=0,1, k=1, \dots, n\\ r_{uv}^{(k,l)}=\frac{\partial^{k+l}(u,v)}{ \partial u^k \partial v^l} \left|_{u=i, v=j} i,j=0,1, k,l=1, \dots, n \end{matrix} $$а также обобщенными функциями смещения $$\alpha_{k ,i}(t) ,$$ где $$0 \le k \le n, i = 0, 1,$$ в качестве которых могут быть взяты произвольные гладкие функции, удовлетворяющие соотношениям
$$\frac{d^m \alpha_{k,i}(t)}{dt^m} \left|_{t=j}=\delta_i^j \delta_k^m=\begin{cases} 1, \mbox{если} i=j, k=m,\\ 0, \mbox{иначе}. \end{cases}$$Здесь $$0\le k, m \le n.$$
Задача 6.3.1. Проверить, что поверхность (6.16), определенная с помощью (6.17) и (6.18), имеет заданные в (6.17) производные на границе.
Замечание 6.3.3. В качестве обобщенных функций смещения всегда можно выбрать многочлены.
Пример 6.3.3. Определим функции смещения формулами
$$\alpha_0(t) = 1-3t^2 + 2t^3, \alpha_1(t) = 3t^2-2t^3,\\ \beta_0(t)=t-2t^2 +t^3 , \beta_1(t) = -t^2 +t^3$$(здесь для простоты использованы обозначения $$\alpha_{0,0}(t) = \alpha_0(t), \alpha_{0,1}(t) = \alpha_1(t), \alpha_{1,0}(t) = \beta_0(t), \alpha_{1, 1}(t) = \beta_ 1(t)) .$$ Пусть заданы , четыре кривые $$r(u, i), r(i,v) ,$$ где $$i = 0, 1,$$ стыкующиеся в точках $$r(0,0), r(0,1), r(1,0), r(1,1) .$$ Пусть заданы первые производные на границе в трансверсальном направлении: $$r_u(i, v), r_v(u, i), i = 0, 1,$$ и вторые смешанные производные в угловых точках: $$r_{uv}(i,j),i,j = 0,1.$$
Кубическая поверхность Кунса определяется формулой (в матричных обозначениях)
$$r(u,v)=(\alpha_0(u) \alpha_1(u) \beta_0(u) \baeta_1(u)) \left( \begin{matrix} r(0,v)\\ r(1,v)\\ r_u(0,v)\\ r_u(1,v) \end{matrix} \right)+\\ +(\alpha_0(v) \alpha_1(v) \beta_0(v) \beta_1(v) ) \left( \begin{matrix} r(u,0)\\ r(u,1)\\ r_v(u,0)\\ r_v(u,1) \end{matrix} \right)-\\ - \left( \begin{matrix} \alpha_0(u)\\ \alpha_1(u)\\ \beta_0(u)\\ \beta_1(u) \end{matrix} \right )^{\top} \begin{pmatrix} r(0,0) r(0,1) r_v(0,0) r_v(0,1)\\ r(1,0) r(1,1) r_v(1,0) r_v(1,1)\\ r_u(0,0) r_u(0,1) r_{uv}(0,0) r_{uv}(0,1)\\ r_u(1,0) r_u(1,1) r_{uv}(1,0) r_{uv}(1,1) \end{pmatrix} \begin{pmatrix} \alpha_0(v)\\ \alpha_1(v)\\ \beta_0(v)\\ \beta_1(v) \end{pmatrix} ,\\ 0 \le u,v \le 1$$Задача 6.3.2. Написать на пакете Mathematica программу построения кубической поверхности Кунса.
В проектировании поверхностей иногда применяется метод построения поверхности по дискретной совокупности ее сечений. Для затягивания остова, образованного сечениями, гладкой поверхностью применяются методы, аналогичные построению интерполяционных кривых. Поэтому получаемые таким способом поверхности имеют такие же названия, как и соответствующие интерполяционные сплайны.
Поверхности Эрмита. Напомним, что сплайном Эрмита называется кусочно-кубическая кривая, проходящая через заданные точки $$p_0, \dots , p n \in R^N$$ и имеющая в этих точках заданные производные $$q_0, \dots , q_n \in R^N.$$ Если производные $$q_i$$ не заданы, то в точках $$p_1, \dots , p_{ n-1}$$ их можно определить из некоторых "естественных" соображений. Одним из примеров такого рода является псевдоупругий сплайн Эрмита. Еще один способ - положить по определению
$$q_i=\frac{p_{i+1}-p_i}{t_{i+1}-t_i}$$(есть и другие способы).
Задание производных на концах $$p_0$$ и $$p_n$$ иногда заменяют на условие обращения в ноль третьей производной на концах кривой. Это условие равносильно следующим формулам:
$$\begin{cases} q_0=2(p_1-p_o)-q_1,\\ q_n=2(p_n-p_{n-1})-q_{n-1} \end{cases} $$Напомним, что вторые производные в опорных точках сплайна Эрмита (где происходит стык кубических кривых) претерпевают разрыв, т. е. сплайны Эрмита принадлежат классу $$C^1.$$ Аналогично сплайнам Эрмита строятся поверхности Эрмита.
Пусть имеется семейство кривых
$$\begin{cases} r=r_0(u)=r(u,0),\\ ..................\\ r=r_i(u)=r(u,i),\\ ..................\\ r=r_n(u)=r(u,n) \end{cases} $$Без ограничения общности будем считать, что параметр u для всех кривых пробегает один и тот же отрезок: $$u_0 \le u \le u_1.$$ Тогда поверхность Эрмита, затягивающая остов кривых $$r(u,i), 0 \le i \le n,$$ определяется формулой, аналогичной формуле для сплайна Эрмита:
если $$j \le v \le j + 1.$$ Здесь $$w = v - j$$ - местный параметр на данном участке поверхности ( $$0 \le w \le 1$$ ), $$r_v(u,j) $$ и $$r_v(u,j + 1) $$ - производные в трансверсальном направлении на заданном остове кривых (6.20), $$\alpha_0(w) = 1-3w^2 + 2w^3, \alpha_1(w) = 3w^2-2w^3, \beta_0( w ) = w-2w^2+w^3, \beta_1(w) = -w^2+w^3.$$ Параметрическая область поверхности Эрмита имеет вид
$$\begin{cases} u_0 \le u \le u_1,\\ 0 \le v \le n \end{cases}$$Если трансверсальные производные $$r_v ( u, i) $$ на линиях остова не заданы, то их, аналогично случаю сплайна Эрмита, можно определить, например, так:
$$\begin{cases} r_v(u,j)=\frac{r(u,j+1)-r(u,j-1)}{2}, 1 \le j \le n-1,\\ r_v(u,0)=2r(u,1)-2r(u,0)-r_v(u,1),\\ r_v(u,n)=2r(u,n)-2r(u,n-1)-r_v(u,n-1). \end{cases} $$При таком определении третьи производные на кривых $$r_0$$ и $$r_n$$ в трансверсальном направлении обратятся в ноль:
$$\frac{\partial^3r(u,v)}{ \partial v^3} \left|_{v=0}=\frac{\partial^3r(u,v)}{ \partial v^3} \left|_{v=n}=0$$Сами исходные кривые $$r_0, \dots ,r_n,$$ образующие остов поверхности Эрмита, могут быть построены как кривые Безье или B -кривые. При этом получаются поверхности класса $$C^1,$$ имеющие на двух участках границы управляемые производные по трансверсальному направлению (а именно, на участках $$r(u,0) $$ и $$r(u,n) $$ ). На участках границы $$r(u_0,v) $$ и $$r(u_1,v) $$ производные, вообще говоря, управляются лишь приблизительно - с помощью выбора кривых остова $$r_i ( u) .$$
Применение поверхностей Эрмита: поверхность перехода. Примером применения поверхностей Эрмита может служить задача объединения двух поверхностей в одну общую поверхность путем соединения их краев с помощью поверхности перехода.
Пусть даны две поверхности M и N, которые надо сопрячь, соединив поверхностью перехода край $$r = r_M(t,u_1^M) $$ первой поверхности и край $$r = r_N(t,u_0^N) $$ второй поверхности. Будем считать, что оба сопрягаемых края отнесены к одному и тому же параметру t, изменяющемуся в пределах $$t_0 \le t \le t_1.$$ Второй параметр на поверхностях M и N обозначим (в обоих случаях) через u. Пусть параметрические пространства поверхностей имеют вид
для поверхности $$M : r = r_M ( t, u ) $$ и
$$\begin{cases} t_0 \le t \le t_1,\\ u_0^N \le u \le u_1^N \end{cases}$$для поверхности $$N : r = r_N(t, u) .$$ Без ограничения общности будем считать, что $$u_1^M \le u_0^N.$$ Тогда поверхность перехода от M к N может быть описана формулой поверхности Эрмита:
где $$w =\frac{u-u_1^M}{u_0^N-u_1^M}$$ - местный параметр по направлению $$u, t_0 \le t \le t_1, 0 \le w \le 1.$$ Данная поверхность Эрмита гладко сопрягает поверхности M и N по заданным краям. Класс гладкости сопряжения $$- C^1.$$
Пример 6.3.4. Приклейка ручки с помощью поверхности перехода:
In[22] := Style [Module [{v, rl, r2, s1, s2 , drl, dr2 , dsl, ds2 , r, s},
r1[u_, v_] : = {-3 + (7 + 1.5 Cos [2 π u]) Sin [2 π v + π] , 1.5 Sin [2 π u] ,
{7 + 1.5 Cos [2 π u]) Cos [2 π v + π] } ; (*тор, u и v в пределах [0,1]*)
r2 [u_, y_] := {5 Cos [2 π u] * Sin [π (v - 1) ] , 5 Sin [2 π] * Sin [π (v - 1) ] ,
5 Cos [π {v - 1) ] } ; {*сфера, u в пределах [0,1], v в пределах [1,2]*)
drl[t_] :=D[rl[t, v], v] * (7 / 6 - (1/2-1/10)) /.v -> l/2-l/10;
dr2[t_] := D[r2[t, v] , v] * (7 / 6 - (1/2-1/10)) / . v -> 7 / 6;
r[t_, w_] :=rl[fc, 1/2-1/10] (l- 3 w2 + 2 w3) + r2[t, 7/6] (з w2 - 2 w3) +
drl[t] (w- 2 w2 + w3) + dr2[t] (-w2 + w3 ) ;
sl[u_, v_] :=r2[u, v + 2] ; (* u в пределах [0,1], v в пределах [-1,0]*)
s2[u_, v_] :=rl[u, v] ;
dsl[t_] := D[sl[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_] :=sl[t, -1/6] (l- 3 w2 + 2 w3) + s2[t, 1/10] (з w2 - 2 w3) +
dsl[t] (w- 2 w2 + w3) + ds2[fc] (-w2 + w3 ) ;
Manipulate[Show[{
PararaetricPlot3D[rl[u, v] , {u, 0, 1} , {v, 1/10, 1 /2 - 1 / 10} ,
PlotStyle -> FaceForm[Orange, Blue]],
PararaetricPlot3D[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, wl},
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
] , {{wl, 1} , 0.001, 1}, {{w2, 0.001}, 0.001, 1}, TrackedSymbols -> {wl, w2},
SaveDefinitions -> True, ControlPlacement -> Lef t] 1 , Magnification -> 2]
Поверхности Лагранжа. Аналогично интерполяционной кривой Лагранжа, можно построить поверхность Лагранжа, затягивающую заданный остов кривых (6.20). Эта поверхность описывается
где
$$L_i(v)=\frac{\Pi_{j=0, j \ne i}^{n}(v-j)}{\Pi_{j=0, j \ne i}^{n}(i-j)}$$Последняя формула получается из формулы для коэффициентов интерполяционного полинома Лагранжа подстановкой $$t_i = i (0 \le i \le n) .$$
Поверхность Лагранжа имеет те же недостатки, что и v поверхность будет бесконечно гладкой и представляет собой многочлен. Гладкость по параметру u зависит от гладкости кривых остова.
Поверхности Гордона. Пусть задан остов поверхности, состоящий из двух семейств кривых в обоих параметрических направлениях, которые образуют сетку - т. е. каждая кривая одного семейства пересекает в одной и только одной точке каждую кривую другого семейства. Различные кривые одного и того же семейства не пересекаются. Обозначим заданную таким образом сетку кривых через
$$r = r(u,0), \dots ,r = r(u,n) \mbox{(первое семейство)},\\ r = r(0,v), \dots ,r = r(m,v) \mbox{ (второе семейство)}.$$Здесь $$0 \le u \le m, 0 \le v \le n.$$ j -й кривой первого семейства с i -й кривой второго семейства равен $$r(i,j), 0 \le i \le m,$$
Определение 6.3.3. Поверхностью Гордона, затягивающей сетку кривых (6.21), называется поверхность
$$r(u,v)=\sum_{j=0}^nL_j^n(v)r(u,j)+ \sum_{i=0}^m L_i^m(u)r(I,v)-\\ - \sum_{i=0}^m \sum_{j=0}^nL_i^m(u)L_j^n(v)r(i,j),\\ 0 \le u \le m,\\ 0 \le v \le n,$$где
$$L_i^m(u)=\frac{\Pi_{r=0, r \ne 0}^m(u-r)}{\Pi_{r=0, r \ne 0}^m(i-r)},\\ L_j^n(v)=\frac{\Pi_{s=0, s \ne j}^n(v-s)}{\Pi_{s=0, s \ne 0}^n(j-s)}$$- коэффициенты Лагранжа.
Обобщенная поверхность Гордона задается формулой (6.22), в которой вместо коэффициентов Лагранжа (6.23) используются произвольные $$C^1$$ -функции $$\alpha_k(w) ,$$ которые в целых точках соответствующего промежутка удовлетворяют условию
$$\begin{cases} 1, \mbox{если } w=k,\\ 0, \mbox{если } w \ne k \end{cases}$$Например, в качестве $$\alpha_k(w) $$ можно взять функции смещения (6.19). Гладкость обобщенной поверхности Гордона зависит от гладкости используемых функций $$\alpha_k(w) .$$ Поверхности Гордона, построенные по коэффициентам Лагранжа, имеют гладкость $$C^{\infty},$$ но могут образовывать необоснованные изгибы.
Поверхности, затягивающие сетку кривых заплатами Кунса. Остов, состоящий из сетки кривых (6.21), можно затянуть поверхностью, состоящей из заплат, затягивающих каждую отдельную ячейку сетки и гладко сопрягающихся между собой. В качестве таких заплат могут выступать поверхности Кунса. Получается следующая поверхность:
$$r(u,v)=(\alpha_0(t) \alpha_1(t) \beta_0(t) \beta_1(t)) \left (\begin{matrix} r(i,v)\\ r(i+1,v)\\ r_u(i, v)\\ r_u(i_1,v) \end{matrix} \right )+\\ +(\alpha_0(w) \alpha_1(w) \beta_0(w) \beta_1(w)) \left (\begin{matrix} r(u,j)\\ r(u, j+1)\\ r_v(u,j)\\ r_v(u,j+1) \end{matrix} \right )-\\ (\alpha_0(t) \alpha_1(t) \beta_0(t) \beta_1(t)) \times\\ \times \begin{pmatrix} r(i,j) r(i,j+1) r_v(i,j) r_v(i,j+1)\\ r(i+1,j) r(i+1, j+1) r_v(i+1,j) r_v(i+1, j+1)\\ r_u(I,j) r_u(i,j+1) r_{uv}(i,j) r_{uv}(i,j+1)\\ r_u(i+1,i) r_u(i+1,j+1) r_{uv}(i+1,j) r_{uv}(i+1,j+1) \end{pmatrix} \times\\ \times \begin{pmatrix} \alpha_0(w)\\ \alpha_1(w)\\ \beta_0(w)\\ \beta_1(w) \end{pmatrix} $$$$t = u ? i, w = v ? j, 0 \le u \le m, 0 \le v \le n.$$ Здесь $$i = [u] $$ - целая часть $$u, j = [v] $$ - целая часть $$v,t = u-i =\{u\}$$ - дробная часть $$u, w = v - j = \{v\}$$ - дробная часть $$v, t,w \in [0,1) $$ - местные параметры в данной ячейке сетки (6.21). Согласно (6.21), области изменения параметров u и v таковы: $$0 \le u \le m, 0 \le v \le n.$$
Для вычисления
при $$1 \le u \le m - 1, 1 \le v \le n - 1$$ (т. е. для внутренних линий сетки) и по формулам
$$r_u(0, v) = 2r(1, v) - 2r(0, v) - r_u(1, v),\\ r_u(m, v) = 2r(m, v) - 2r(m - 1, v) - r_u(m - 1, v),\\ r_v(u,0) = 2r(u, 1) - 2r(u, 0) - r_v(u, 1),\\ r_v (u, n) = 2r(u, n) - 2r(u, n - 1) - r_v (u, n - 1)$$для крайних линий сетки.
Геометрически это означает, что производная $$r_u(i,v) $$ в произвольной точке $$r(i,v) $$ внутренней линии сетки полагается пропорциональной вектору, соединяющему точки $$r(i - 1,v) $$ и $$r(i + 1,v) $$ на соседних с ней линиях сетки, отвечающие тому же значению параметра v, что и данная точка. Аналогичное утверждение справедливо и для производной $$r_v(u,j) .$$
Так же как и в случае кубических сплайнов Эрмита, формулы (6.25) обеспечивают равенство нулю третьих производных поверхности на границе в трансверсальном к граничной кривой направлении.
Поверхности тензорного произведения. Рассмотрим обобщенную поверхность Гордона (6.22), построенную по сетке кривых, которые являются
Для того чтобы получить формулу поверхности тензорного произведения, поставим в формулу (6.22) выражения, определяющие u и v:
где, напомним:
$$\alpha_0(x) = 1 - 3x^2 + 2x^3,\\ \alpha_1(x) = 3x^2-2x^3,\\ \beta_0(x) = x-2x^2+x^3,\\ \beta^1(x) = -x^2 + x^3,$$где $$x \in [0,1] $$ - местный параметр. В формулах (6.26) через $$i $$ и $$j$$ обозначены целые части параметров u и v соответственно. Местные параметры $$t$$ и $$w$$ определяются как $$t = u - i, w = v - j.$$ Вторые две формулы в (6.26) получаются дифференцированием формулы (6.22). После подстановки формул (6.26) в (6.22) и несложных алгебраических преобразований получим следующее выражение для
Из формулы (6.27) хорошо видно, что поверхность тензорного произведения определяется величинами, заданными лишь в узлах сетки.
При триангуляции поверхностей иногда бывает удобно рассматривать криволинейные поверхности треугольной формы, имеющие в качестве параметрической области треугольник на плоскости. При этом удобно пользоваться не декартовыми, а барицентрическими параметрическими координатами.
Лемма 6.6. Пусть даны три точки $$A, B, C \in R^2,$$ не лежащие на одной прямой. Тогда для любой точки $$P \in R^2$$ существует единственная тройка чисел $$\{a, b, c\}$$ таких, что
$$\begin{cases} r_P= ar_A + br_B + cr_C, \\ a + b + c=1. \end{cases}$$Здесь через $$r_P, r_A, r_B, r_C$$ обозначены радиус-векторы точек $$P, A, B, C$$ соответственно.
Числа $$a, b, c$$ называются барицентрическими координатами точки P относительно треугольника $$\Delta ABC.$$
Следствие 6.5. Имеют место формулы
$$|a|=\frac{S_{PBC}}{S_{ABC}},\\ |b|=\frac{S_{APC}}{S_{ABC}},\\ |c|=\frac{S_{ABP}}{S_{ABC}} $$т. е. коэффициенты $$a, b, c$$ по модулю равны отношению площади треугольника $$\Delta ABC,$$ в котором соответствующая вершина заменена на $$P,$$ к площади исходного треугольника $$\Delta ABC.$$
Следствие 6.6. Все три барицентрические координаты точки положительны тогда и только тогда, когда точка принадлежит внутренности треугольника $$\Delta ABC.$$ Все три барицентрические координаты точки неотрицательны тогда и только тогда, когда точка принадлежит треугольнику $$\Delta ABC.$$
Билинейная треугольная поверхность. Пусть даны три точки $$P_1, P_2, P_3$$ в $$R^3.$$ В качестве области параметров возьмем треугольник $$\Delta ABC \in R^2.$$ Точка параметрического пространства описывается координатами $$0 \le a, b, c \le 1, a+b+c = 1.$$
Треугольная поверхность на трех кривых. Пусть даны три попарно пересекающиеся кривые $$r = r_1(u), r = r_2(v), r = r_3(w ), 0 \le u,v,w \le 1,$$ пересекающиеся в точках $$P_1 = r_2(0) = r_3(1), P_2 = r_1(1) = r_3(0), P_3 = r_1(0) = r_2(1) .$$ Определим треугольную поверхность формулой
$$r(a, b, c) = a (r_3(1 - b) + r_2(c) - Р_1) + b (r_1(1 - c) + r_3(a) - P_2 ) + c (r_2(1-a) + r_1(b)-P_3 ),$$где $$0 \le a, b, c \le 1, a + b + c = 1.$$ Легко проверить, что при $$c = 0$$ формула (6.28) описывает кривую $$r_3(a) (0 \le a \le 1) ,$$ при $$a = 0$$ - кривую $$r_1(b) (0 \le b \le 1) ,$$ при $$b = 0$$ - кривую $$r_2(c) (0 \le c \le 1) .$$ Например, если $$a = 0,$$ то $$c = 1 - b$$ и, согласно (6.28), имеем
$$r(0, b,1 - b) = b (r_1(b) + r_3(0) - P_2 ) + (1 - b) (r_2(1) + r_1(b) - P_3 ) =\\ = b (r_1(b) + P_2 -P_2 ) + (1 - b) ( P_3 + r_1(b) -P_3) = r_1(b).$$Пример 6.3.5. Поверхность, соединяющая три поверхности по кривым на их границах. При изменении поверхностей, имеющих вид "треугольника", построенная поверхность не изменяется:
In[23]: =
DynamicModule[{pe200, pe002, pel0l, te, tee, pe020, pell0, pe0ll, r2 , rf2 ,
rf1, u, v, w, s, t, рb020, рb002, pb0ll, tb, tbb, pb200, pbl0l, pbll0, rl,
pj200, pj020, pjll0, tj , tjj, pj002, pjl0l, pj0ll, гЗ, rf3, flperem, f2perem,
f3perem, rfel},
pe200 = {0, 10, 0} ; pe002 = {10, 0, 0} ; pel01 = {10, 10, 0}; pb020 = {0, 0, 20};
pb002 = {10, 0, 0} ; pb0ll = {10, 0, 20} ; pj200 = {0, 10, 0}; pj020 = {0, 0, 20};
pjll0 = {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] ;
pell0[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 pel0l u w + 2 pel 10 [ te] u v +
2 pe011[te] vw;
rf2[ te_, s_, t_] := r2 [te] /.w ->l-u-v/. {u -> t- s t, v-> s t} ;
f2perem = ParametricPlot3D[Evaluate[r2 [te] /. w -> l-u-v/. v->0], {u, 0, 1},
Mesh-> 18, PlotStyle -> {Red, AbsoluteThickness [2] } ] ;
pb020 = {0, 0, 20} ; pb002 = {10, 0, 0} ; pb0ll = {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] ;
pbl01[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] ;
pbll0[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] ;
rl [ tb_] : = pb200 [ tb] u u + pb020 v v + pb002 w w + 2 pbl0l [ tb] uw + 2 pbll0 [ tb] u v +
2 pb0ll vw;
rfl [ tb_, s_, t_] := rl [tb] /.w->l-u-v/. {u -> t- s t, v-> s t} ;
flperem = ParametricPlot3D [Evaluate [rl [tb] /. w -> l - u - v / . u -> 0] , {v, 0, 1},
Mesh -> 18, PlotStyle -> {Red, AbsoluteThickness[2]}];
pj200 = {0, 10, 0} ; pj020 = {0, 0, 20} ; pjll0 = {0 , 10, 20};
pj002[tj_] := If [tj<= 0.5, {6, 8, 16} + 2 tj j ({0, 20, 40} - {6, 8, 16}) /. tj j -> tj,
{0, 20, 40} + 2 (tj j - 0.5) ({-6, 8, 16} - {0, 20, 40}) /. tj j -> tj] ;
pjl01[tj_] := If {tj <= 0,5, {5, 9, 3} + 2 tjj ({0, 15, 20}-{5, 9, 8}) /. tj j -> tj,
{0, 15, 20} + 2 (tjj - 0.5) ({-5, 9, 8} - {0, 15, 20}) / . tj j -> tj] ;
pj0ll [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}) /. tj j -> tj] ;
r3[ tj_] : = pj200uu+pj020vv + pj002[fcj] v w + 2 pjl0l[tj] u w + 2 pjll0 u v +
2 pj0ll [tj] vw;
rf3[tj_, s_, t_] := r3[tj] /.w->l-u-v/. {u -> t- s t, v-> s t} ;
f3perem = ParametricPlot3D [Evaluate [r3 [tj ] /. w -> 0 / . v->l-u], {u, 0, 1},
Mesh ->18, PlotStyle -> {Red, AbsoluteThickness[2]}];
rfel[tb_, te_, tj_, u_, v_] : =
и (rf3[tj, v, 1] + rf2[te, 0, u + v] - rf2[te, 0, 1] ) +
v (rfl [tb, 1, u+ v] + rf3[tj, 1 - u, 1] - rfl[ tb, 1, 1] ) +
(1-u- v) (rf2[te, 0, u] + rfl[ tb, 1, v] - rf 1 [ tb, 1, 0] ) ; (*три кривые rfl[tb,l,t] ,
rf2[te,0,l-t] и rЈ3[tj,l-t,l] , t Ha [0,1]*)
Manipulate[
Show[{ParametricPlot3D[Evaluate[rfl[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 . 9] , 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, t j , s, t] ] , {s, 0, 1}, {t, 0, 1 - s} ,
Mesh -> {13, 18}, PlotRange-> All] ,
f1perem, f2perem, fЗрегеm} , Viewpoint-> {-0 . 7 , 0.4, 1-3}, ImageSize -> {300 , 200},
Boxed -> False, Axes -> None, PlotRange -> {{-6. , 20.}, {-6., 20.}, {-10., 40.}},
Spherical Region -> True, ViewAngle-> 0 . 7 ],
{{tb, 0.4, "Зеленый ""треугольник '"}, 0., 1., Appearance->'Labeled" } ,
{ {te, 0.4, " Желтый ' ' треугольник '" } , 0 . , 1. , Appearance -> " Labeled" } ,
{ {t j , 0.4, " Синий " " треугольник' ' " } , 0 . , 1. , Appearance -> " Labeled" } ,
SaveDefinitions -> True] ]
Основными недостатками кривых Безье являются следующие:
Поэтому при работе с кривыми Безье часто приходиться прибегать к делению кривой Безье, т. е. по сути дела - к замене одной кривой другими. Поэтому естественным образом возникает задача построения такого класса кривых, который бы позволял избежать недостатков 1), 2) и решать соответствующие задачи управления кривой путем изменения управляющих параметров кривой, а не путем ее замены на другие кривые.
Вспомним, что рациональная кривая Безье, построенная по точкам $$p_0, \dots, p_n$$ с весами $$\omega_0, \dots , \omega_n,$$ определялась формулой
$$r_0(t;n)=\frac{\sum_{i=0}^{n}C_n^it^i(1-t)^{n-i} \omega_ip_i}{\sum_{i=0}^n C_n^it^i(1-t)^{n-i} \omega_i}=\frac{\sum_{i=0}^n B_i^n(t) \omega_ip_i}{\sum_{i=0}^n B_i^n(t) \omega_i},$$где $$\{B_i^n(t) = C_n^ii t^i (1 -t)^{n-i}\}$$ - базис Бернштейна на [0,1].
Из формулы (6.1) видно, что глобальное влияние каждой опорной точки $$p_i$$ на всю кривую происходит из-за того, что каждая из функций $$B_i^n$$ не равна нулю на всем интервале (0,1). Если бы вместо функций Бернштейна в формуле (6.1) выступали функции с локальными носителями, существенно меньшими, чем область определения параметра кривой $$t,$$ то удалось бы как локализовать влияние отдельной точки на всю кривую, так и добиться (за счет выбора таких функций) прохождения кривой через некоторую заданную точку.
Поэтому поставим задачу следующим образом. Пусть даны опорные точки $$p_1,\dots., p_n$$ с весами $$\omega_1,\dots ,\omega_n.$$ Определим кривую $$r(t) $$ по формуле
$$r(t)=\frac{\sum_{i=1}^nN_i(t) \omega_ip_i}{\sum_{i=1}^nN_i(t) \omega_i}, t_{min} \le t \le t_{max},$$где $$N_i ( t) $$ - набор кусочно-полиномиальных функций таких, что
для любой гладкой функции $$f(t) $$ существует линейная комбинация $$\sum c_iN_i(t) ,$$ интерполирующая функцию $$f(t) $$ в заданных узлах $$\{t_i\}_{i=1}^n$$
$$f(t_i)=\sum_{j=1}^nc_jN_j(t_i), i=1, \dots, n;$$Последнее условие вводится для того, чтобы в случае совпадения всех опорных точек кривая превращалась бы в ту же точку.
Решение поставленной задачи основывается на понятии
Пусть имеется гладкая на каждом отрезке $$[t_i,t_{i+1}]$$ функция $$f(t) $$ и фиксированы некоторые значения $$t_1, \dots , t_{m+1}$$ переменной $$t$$ (не обязательно различные и не обязательно упорядоченные по возрастанию или убыванию $$t$$ ). Мы также будем рассматривать переменное значение $$t,$$ которое будем считать первым в последовательности. Поскольку оно предшествует значению $$t_1,$$ припишем ему индекс 0 (хотя в большинстве случаев будем обозначать его просто $$t$$ без индекса): $$t_0 = t, t_1, \dots, t_{m+1}.$$ Обозначим $$f_i = f(t_i), f_t = f(t) = f(t_0) = f_0.$$ Определим оператор правого сдвига $$E$$:
$$Et_i = t_{i+1},\\ Ef_i = f_{i+1}\\ (0 \le i \le m).$$Определение 6.1.1. Разделенной разностью функции $$f(t) $$ в узле $$t_i$$ порядка $$k \le m - i + 1$$ называется величина
$$R^{(k)}f_i=\frac{R^{(k-1)}(E-1)f_i}{(E^k-1)t_i}=\frac{R^{(k-1)}f_{i+1}-R^{(k-1)}f_i}{t_{i+k}-t_i},\\ R^{(0)}f_i=f_i.$$При этом если $$t_i = t_{i+1} = \dots = t_{i+k},$$ то в правой части (6.3) появляется отношение вида $$\frac 00,$$ которое по определению полагается равным
$$R^{(k)}f_i=lim_{\delta_1, \dots, \delta_k \to 0, \delta_i > 0} R_{\delta_1, \dots, \delta_k}^{(k)}f_i$$где $$R_{\delta_1, \dots, \delta_k}^{(k)}f_i $$ определяется по возмущенным узлам $$\tilde t_{i+j}=t_i+\delta_1 + \dots, + \delta_j,$$ где $$1 \le j \le k, \delta_j > 0 .$$
Теорема 6.1. Если $$t_i = t_{i+1} = \dots = t_{i+k},$$ то
Пример 6.1.1.
In[1]:=
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 [l/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[[.#l]] , {-20, 20, .01}, Appearance -> "UpArrow" ,
ImageSize -> {200, 12}], Dynamic@ss [[#1]] } /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols -> {ss} ,
SaveDef initions -> True ] ]
Задача 6.1.1. Обратите внимание, что
Задача 6.1.2. Написать программу, которая бы учитывала совпадение узлов.
Пусть $$m \ge 0$$ - целое число. Рассмотрим функцию $$\sigma_m(z,t) = ( max(0, z -t) )^m.$$ При фиксированном $$t$$ функция $$\sigma_m(z,t),$$ рассматриваемая как функция $$z,$$ называется усеченной степенной функцией.
Лемма 6.1. Для производных усеченной степенной функции справедливы формулы
$$\frac{d^k \sigma_m}{dz^k}=\frac{m!}{(m-k)!}\sigma_{m-k}, \frac{d^k \sigma_m}{dt^k}=\frac{(-1)^km!}{(m-k)!}\sigma_{m-k}, k=1, \dots, m$$Обозначим через $$f[t_1,\dots, t_k] $$ разделенную разность порядка $$(k-1) $$ функции $$f(x) ,$$ построенную по узлам $$t_1,\dots, t_k.$$ Фиксируем набор узлов $$t_0 < t_1 \lt; \dots < t_m$$ и рассмотрим функцию $$M(t) =\sigma_{m-1}[t_0, \dots, t_m](t) ,$$ где
Лемма 6.2. Пусть $$m > 0.$$ Тогда $$M(t) = 0,$$ если $$t\le t_0$$ или $$t \ge t_m.$$
Замечание 6.1.1. Из леммы 6.2 вытекает, что $$\sigma_{m-1}[t_0, \dots ,t_m](t) $$ является непрерывной при $$m \ge 2$$ функцией параметра $$t,$$ отличной от нуля лишь на $$(t_0,t_m) .$$
Пример 6.1.2.
In[2]:=f[t_, z_]:=If[t < z, (z-t)^3, 0];
Plot[f[-6, z], {z, -20, 20}]
In[4] :=DynamicModule[ {h, z, t, st, M, G, i, x, y} ,
h[x_, y_] : = If [x < y, (y - x)^3, О] ;
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- l]] = 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[l, Length[st]]},
Grid[{Spacer[3], Subscript[
Style["t", Italic], #-1],
Slider [Dynamic@st[[#1]] , {-15, 13, .01},
Appearance -> "UpArrow" , ImageSize -> {200, 12} ] ,
Dynamic@st|[#l]l} /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols ->{st},
SaveDefinitions -> True] ]
Задача 6.1.3. Написать программу, которая бы учитывала совпадение узлов.
На основе
Определение 6.1.2. Нормированным B-сплайном m-го порядка $$N_{i,m}(t) $$ для неубывающей последовательности узлов $$t_i \le t_{i+1} \le \dots t_{i+m},$$ отсчитываемой от первого узла $$t_i,$$ называется функция
$$N_{i,m}(t)=(t_{i+m}-t_i) \sigma_{m-1}[t_i, \dots, t_{i+m}](t)$$Соответственно, ненормированным B-сплайном $$M_{i,m}(t) $$ называется функция
$$M_{i,m}(t)=\sigma_{m-1}[t_i,\dots, t_{i+m}](t).$$Определение 6.1.3. Нормированным B-сплайном m-го порядка $$N_{m,i}(t) $$ для неубывающей последовательности узлов $$t_{i-m} \le t_{i-m+1} \le \dots \le t_{i-1} \le t_i,$$ отсчитываемой от последнего узла $$t_i,$$ называется функция
$$N_{m,i}(t) = (t_i -t_{i-m})\sigma_{m-1}[t_{i-m}, \dots ,t_i](t).$$Соответственно определяется и ненормированный B-сплайн $$M_{m,i}(t) $$:
$$M_{m,i}(t) =\sigma_{m-1}[t_{i-m}, \dots ,t_i](t).$$Замечание 6.1.2. Если порядок B -сплайна определять не числом узлов минус единица, а степенью усеченной функции (как многочлена t ), то введенные выше B -сплайны будут иметь (по определению) порядок не m, а (m - 1). В этом случае для них применяются другие обозначения: $$N_{-i}^{m-1}, N_{i-}^{m-1}, M_{i-}^{m-1}, M_{-i}^{m-1}.$$ При этом
И те и другие обозначения присутствуют в литературе по B -кривым. Мы будем пользоваться первыми. В Mathematica используются вторые.
Рассмотрим (для простоты) бесконечную в обе стороны последовательность неубывающих узлов $$\{t_i\},$$ где $$t_i \in R, i \in Z, t_i \le t_{i+1}$$ для любого $$i.$$ Фиксируем $$m > 0.$$ Тогда определены множества B -сплайнов: сплайны $$N_{i, m}(t) $$ привязаны к узлу $$i$$ как к первому узлу, сплайны $$N_{m,i}(t) $$ привязаны к узлу $$i$$ как к последнему узлу $$(i \in Z) .$$
Лемма 6.3. Если $$t \in [t_j,t_{j+1}],$$ то среди функций $$N_{i,m}(t) $$ и $$N_{m,i}(t) $$ только m функций в каждом из этих двух классов могут быть отличны от нуля в точке $$t.$$ А именно, отличны от нуля в точке $$t$$ могут быть лишь следующие функции:
$$N_{j-m+1,m}( t ), \dots ,N_{j,m}(t),$$и
$$N_{m,j+1}(t), \dots , N_{m,j+m}(t).$$Лемма 6.4. Для любого $$t \in R$$ имеет место
$$\sum_i N_{i, m}(t)=\sum_iN_{m,i}(t)\equiv 1$$Теорема 6.2 (формула Кокса - де Бура). Справедливы следующие формулы для ненормированных B-сплайнов при любых $$t, m \ge 2, \alpha \in Z$$ и $$\beta = \alpha + m$$:
$$M_{\alpha, m}(t)\frac{(t_{beta}-t)M_{\alpha+1, m-1}(t)+(t-t_{\alpha})M_{\alpha, m-1}(t)}{t_{\beta}-t_{\alpha}},\\ M_{m, \beta}(t)=\frac{(t_{\beta}-t)M_{m-1, \beta}(t)+(t-t_{\alpha})M_{m-1, \beta-1}(t)}{t_{\beta}-t_{\alpha}}\\ M_{\alpha-}^{m-1}(t)=\frac{(t_{\beta}-t)M_{\alpha+1-}^{m-2}(t)+(t-t_{\alpha})M_{\alpha-}^{m-2}(t)}{t_{\beta}-t_{\alpha}}\\ M_{-\beta}^{m-1}(t)=\frac{(t_{\beta}-t)M_{-\beta}^{m-2}(t)+(t-t_{\alpha}M_{-\beta-1}^{m-2}(t)}{t_{\beta}-t_{\alpha}}$$Здесь ненормированные B-сплайны первого порядка равны по определению разделенным разностям функции $$\sigma_0(z) = (z - t)_+^0$$:
$$M_{i,1}(t)=M_{i-}^0=\begin{cases} \frac{1}{t_{i+1}-t_i}, \qquad \mbox{если}t \in [t_i, t_{i+1}),\\ 0, \qquad \mbox{иначе}; \end{cases}\\ M_{1,i}(t)=M_{-i}^0=\begin{cases} \frac{1}{t_i-t_{i-1}}, \qquad \mbox{если } t \in [t_{i-1}, t_i), \\ 0, \qquad \mbox{иначе}. \end{cases}$$Теорема Кокса - де Бура позволяет вычислять B -сплайны рекуррентно с использованием формул (6.4) и (6.5).
Следствие 6.1. Пусть имеется произвольная последовательность (не обязательно неубывающих) узлов $$T = \{_t0, \dots ,t_m\}.$$ Пусть среди них есть 2 неравных друг другу узла $$t_a$$ и $$t_b: t_a \ne t_b,$$ где $$0 \le a,b \le m.$$ Обозначим
$$M^T(t) =\sigma_{m-1}[T] = \sigma_{m-1}[t_0,\dots ,_tm],\\ M^{T/a}(t) = \sigma_{m-2}[T\ a] = \sigma_{m-2}[t_0, \dots ,\hat t_a, \dots ,t_m],\\ M^{T/b}(t) =\sigma_{m-2}[T\b] =\sigma_{m-2}[t_0, \dots ,\hat t_b, \dots ,t_m],\\ p_{ab}(t) =\frac{t_b-t}{t_b-t_a}, q_{ab}(t)=\frac{t-t_a}{t_b-t_a}, (p_{ab}+q_{ab} \evqiv 1) $$где "домиком" отмечены отброшенные узлы. Тогда
$$M^T(t) = p_{ab}(t)M^{T\a}(t) + q_{ab}(t)M^{T\b}(t),\\ M ^{\{t_a, t_b\}} (t)=\begin{cases} \frac{1}{|t_b-t_a|}, \qquad \mbox{если } t \in [t_a \wedge t_b, t_a \vee t_b),\\ 0, \qquad \mbox{иначе}, \end{cases}$$где $$\alpha \vee \beta = max\{\alpha, \beta \}, \alpha \wedge \beta = min\{\alpha, \beta \}.$$
Следствие 6.2. B - сплайны - неотрицательные функции:
Следствие 6.3. Для нормированных B-сплайнов формулы Кокса - де Бура принимают вид $$(\beta = \alpha + m) $$:
$$N_{\alpha,m}(t)=\frac{t_{\beta}-t}{t_{\beta}-t_{\alpha+1}}N_{\alpha+1, m-1}(t)+\frac{t-t_{\alpha}}{t_{\beta-1}-t_{\alpha}}N_{\alpha, m-1}(t),\\ N_{m, \beta}(t)=\frac{t_{\beta}-t}}{t_{\beta}-t_{\alpha+1}}N_{m-1, \beta}(t)+\frac{t-t_{\alpha}}{t_{\beta - 1}-t_{\alpha}}N_{m-1, \beta-1}(t),\\ N_{i-}^0(t)=N_{i,1}(t)=\begin{cases} 1, \qquad \mbox{если } t \in [t_i, t_{i+1})\\ 0, \qquad \mbox{иначе;} \end{cases}\\ n_{-i}^0(t)=N_{1,i}=\begin{cases} 1, \qquad \mbox{если } t \in [t_i, t_{i+1}),\\ 0, \qquad \mbox{иначе}. \end{cases} $$Для сокращения объема вычислений следует сначала вычислить ненормированные B -сплайны нужного порядка по формулам Кокса - де Бура, а затем нормировать их. Если требуется вычислить все ненулевые B -сплайны порядка m для $$t \in [t_i,t_{i+1}]$$ в неубывающей последовательности узлов $$\{t_i\},$$ то вычисления производятся по схеме
Аналогичная схема вычислений применяется для B -сплайнов остальных трех типов.
Следствие 6.4. Для любых $$m \ge 2$$ и $$\alpha \in Z$$ функции $$N_{\alpha,m}(t) $$ и $$N_{m, \beta}(t) ,$$ где $$\beta = \alpha + m,$$ построенные по неубывающей последовательности узлов $$\{t_i\}_{i \in Z},$$ являются непрерывными $$(m - 2) $$ раза дифференцируемыми функциями, обращающимися в ноль вне отрезка $$[t_{\alpha},t_{\beta}].$$ На отрезке $$[t_i,t_{i+1}]$$ они имеют производные до порядка $$(m - 1) $$ включительно. Функции $$N_{\alpha, m}(t) $$ и $$N_{m, \beta}(t) $$ зависят (как от параметров) лишь от значений узлов $$t_{\alpha}, \dots , t_{\beta}.$$
Задача 6.1.4. С помощью встроенной функции в пакет Mathematica $$BSplineBasis[\{m,{t_0, \dots ,t_n\}\}, i, t] ,$$ вычисляющей тот из $$n - m$$ нормированныx B -сплайнов $$N_{i, m+1}(t) $$ порядка $$m + 1$$ (или m при втором способе определения порядка, см. выше), который привязан справа к i -му ( $$0 \le i \le n - (m + 1) $$ ) узлу из расширенного множества неубывающих узлов $$t_0 \le t_1 \le \dots \le t_n$$ (см. определение ниже), наглядно изучить ненормированные и нормированные B -сплайны, привязанные к узлам справа. Обратите внимание, что в Mathematica узлы в расширенном множестве нумеруются начиная с нуля. Соответственно, в Mathematica B -сплайны в базисе нумеруются с нулевого до n - (m + 1) -го.
Замечание 6.1.3. При работе с нижеследующей программой обратите внимание, что при совпадении первых m + 1 узлов значение нулевого (самого первого) нормированного B -сплайна из базиса в нулевом (самом первом) узле становится равно 1, тогда как все остальные B -сплайны базиса обращаются в этом узле в ноль. Это значит, что когда такой базис из нормированных B -сплайнов используется для вычисления B -кривой (аналогично базису Бернштейна для кривых Безье), получившаяся кривая будет начинаться в своей крайней левой опорной точке. Аналогично, при совпадении последних m + 1 узлов расширенного множества, последний сплайн базиса в последнем узле примет значение 1, а остальные будут равны в нем нулю, и, следовательно, правый конец соответствующей кривой будет находиться в ее последней опорной точке. См. ниже теорему 6.3, в формулировке которой порядок B -сплайна m соответствует значению m
+ 1 в обозначениях Mathematica (как уже говорилось, Mathematica использует второй способ определения порядка B -сплайна, а значит, уменьшает его на единицу по сравнению с тем способом, который используем мы).
Задача 6.1.5. Написать программу на пакете Mathematica, вычисляющую B -сплайны $$N_{m,i}(t) ,$$ привязанные к узлам слева.
Пример 6.1.3. Следующие две программы строят графики B -сплайнов, привязанных к узлам справа. В данном примере n = 9, то есть расширенное множество узлов состоит из 10 точек, начиная с нулевой, а m = 3, то есть порядок B -сплайна при принятом нами определении равен 4.
In[5]:=
DynamicModule [ {m = 3, i, tt, knots} ,
For [i = 1, i Sm + 1, i++, tt [i] = 0] ; (*m+l начальных узло равны 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 , 2 m + 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 [N0,m+i [t] , 14] ] ,
Do [If [ss[[i - 1]] > ss[[1]] , ss[[i - 1]] = ss[[i]] , {i, 2, Length [ss] }] ,
Plot[BSplineBasis[{m, knots}, 0, t],
{t, knots[[1]1 , knots[[Length[knots]]] +1} ,
PlotRange -> + {{0, 1.1}, {0, 1. 2} }], PlotLabel -> Style [N0, 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[[#l]] , {0, 1, .01}, Appearance -> " UpArrow" ,
ImageSize -> {200, 12}], Dynamic@ss[el]]} /@ r] ] ,
ControlPlacement -> Bottom} , TrackedSymbols -> {ss} , SaveDefinitions -> True 1 ]]
In|6] :=DynamicModule [{m = 3, n = 9, i, j , tt, knots, a) ,
For [i = 1, i <= m+1, i + +, tt[i] = 0] ;
(*m+l начальных уело равны 0*)
For[j = 1, j im + 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;
Show[If[OrderedQ[ss], Plot[Evaluate[
Table[BSplineBasis[{m, knots}, 1, t] ,
{1, 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 [Hi, 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}, 1, t] ,
{1, 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 [Ni, т+1 [t] , 14]] ] ,
Graphics[Arrow[{{2.5, 1.1},
, {i + 0.1, BSplineBasis[{m, knots}, i, i + 0.1]}}] ] ,
], {{i, 3, " Номер полинома 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[[#l]], {0, 2, .01},
Appearance -> "UpArrow" , ImageSize -> {200 , 12} ] ,
Dynamic@ss|[#l]]} /@r]] ,
ControlPlacement -> Bottom} ,
{{fill, True, "Заполнение"}, {True, False),
ControlPlacement -> Bottom} , SaveDefinitions -> True] ]
Определение 6.1.4. B - кривой порядка m, построенной по вершинам $$p_1, \dots , p_n, n \ge m,$$ с весами $$\omega_1, \dots , \omega_n,$$ называется кривая
определяемая опорными точками $$\{p_i\}_{i=1}^n,$$ весами $$\{\omega_i\}_{i=1}^n$$ и последовательностью узлов $$t_1, \dots ,t_{n+m}.$$
Определение 6.1.5. Множество узлов $$\{t_1, \dots , t_{n+m}\}$$ называется расширенным множеством узлов B -кривой.
Замечание 6.1.4.
Управляющими параметрами B -кривой являются:
m,неубывающая последовательность вещественных чисел $$\{t_1, \dots , _{tn+m}\}.$$
Таким образом, B -кривые обладают дополнительным множеством управляющих параметров по сравнению с рациональными кривыми Безье - расширенным множеством узлов. Именно эти дополнительные параметры и обеспечивают гораздо большую гибкость в управлении B -кривыми, чем кривыми Безье. В современной компьютерной геометрии в основном применяются именно B -кривые.
B -кривые, вообще говоря, не проходят через свои крайние опорные точки. Однако такого прохождения нетрудно добиться с помощью следующего утверждения.Теорема 6.3. Если первые m узлов расширенного множества совпадают, то B -кривая начинается в своей первой опорной точке $$p_1.$$ Аналогично, если последние m узлов расширенного множества совпадают, то B -кривая заканчивается в своей последней опорной точке $$p_n.$$
B -кривые в пакете Mathematica. Mathematica имеет встроенные функции BSplineFunction и BSplineCurve, которые позволяют строить B -кривые по заданному множеству опорных точек пространства $$R^n, n = 2,3,$$ автоматическим или полуавтоматическим подбором параметров. Расширенное множество узлов по умолчанию подбирается автоматически, однако его можно задать самому с помощью параметра SplineKnots. Например, $$SplineKnots \to Automatic$$ означает автоматический подбор расширенного множества узлов, а $$SplineKnots \to \{0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 2\}$$ задает расширенное множество узлов непосредственно.
Так же дело обстоит и с весами опорных точек. Веса можно задавать с помощью SplineWeights. Например, $$SplineWeights \to Automatic$$ означает автоматический выбор весов, а $$SplineWeights \to \{1, 1, 10, 1, 1, 1\}$$ задает конкретные веса для шести опорных точек.
Порядок сплайнов можно задавать с помощью SplineDegree, но при этом должно быть $$SplineKnots \to Automatic.$$
Вычисление порядка происходит так. Число сплайнов в базисе, равное n - m (см. выше), должно равняться числу опорных точек. Поэтому порядок сплайнов может быть вычислен как разность числа узлов и числа опорных точек. Однако для Mathematica эта разность должна быть подобрана так, чтобы она равнялась 4, иначе B -кривая может быть построена в Mathematica неправильно.
Чтобы воспользоваться возможностью полноценного управления B-кривой, в том числе и весов, нужно запрограммировать B -кривую самостоятельно с помощью встроенной функции BSplineBasis.
Задача 6.1.6. С помощью встроенных функций Mathematica SplineFunction и BSplineCurve построить B -кривую по заданным опорным точкам с автоматическим выбором узлов на плоскости и в трехмерном пространстве.
Пример 6.1.4. B -кривая на плоскости с автоматическим выбором узлов, построенная двумя способами: с помощью функций BSplineFunction и BSplineCurve соответственно. На втором графике для сравнения нарисована также простая кривая Безье, построенная по тем же опорным точкам.
In[7] : = DynamicModule [ {ptsO, f, g, n, ptsl, pts2 , t} ,
n = Length [ptsO] ;
Manipulate[f = BSplineFunction[ptsl] ;
Show[{
Graphics[{Green, Line[ptsl], Black,
Text["B-кривая, построенная через BSplineFunction", {12, 10-0}],
Text[ToString[# - 1] , ptsl [ [#] ] +{0.7, 0}] /@ Range [n] } ,
PlotRange -> {{-4, 20}, {-4, 10}}, Axes -> True] ,
ParametricPlot[f[t], {t, 0, 1}]}, ImageSize -> 300] ,
{{ptsl, 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}})]
In[8] :=
DynamicModule [ {ptsO , f, g, n, ptsl, pts2 , t} ,
n = Length [ptsO] ;
Manipulate[g = BezierFunction[pts2];
Show [ {
Graphics [ {Red, BSplineCurve [pts2 , SplineKnots -> Automatic] , Black,
Text["В-кривая, построенная через BSplineCurve (красная)", {10, 10.0}],
Техt["Криая Беэье (голубая)", {10, 9.0}],
Text [ToSt ring] # - 1] , pts2 [ [#] ] + {0.7, 0}] /@ Range [n] } ,
PlotRange -> {{-4, 20} , {-4, 10}} , Axes -> True] ,
ParametricPlot[g[t] , {t, 0, 1}, ElotStyle -> {Blue}]
}, ImageSize -> 300] , { {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}})]
Пример 6.1.5. B -кривая в пространстве с автоматическим выбором узлов в сравнении с кривой Безье, построенной по тем же опорным точкам. Обратите внимание, насколько более послушно ведет себя B -кривая при управлении опорными точками, чем кривая Безье. Кривая Безье по-другому - более глобально - усредняет влияние опорных точек на кривую, что делает ее локальное управление за счет смещения отдельных опорных точек гораздо более вязким, чем в случае B -кривой.
In[9]:=
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["t1", Italic] , #- l] ,
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 ["t2 и t3", Italic], #-l],
Slider2D [Dynamic® {pts[[1]] [ [2] ] , pts[[1]] [ [3] ] } ,
{{-15, -15}, {15, 15}}, Appearance -> "Left",
ImageSize -> {50, 30}],
Dynamic@ {pts [[#1]] [ [2] ] , pts [[#2]] [[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.1.6. B -кривая на плоскости с автоматическим выбором порядка, но с управляемыми мышью опорными точками и узлами (узлы должны быть неубывающими).
In[10]:=
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[l, m] } ,
Grid[
{Spacer[3] , Subscript [Style ["t" , Italic] , # - 1] ,
Slider [Dynamic@tt[[#l]] , {0, 15, .01}, Appearance -> "UpArrow" ,
ImageSize -> {200, 12}], Dynamic@tt[[i]]} /@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])]
Пример 6.1.7. B -кривая на плоскости с автоматическим выбором порядка и узлов, но с управляемыми мышью опорными точками и весами, и простая кривая Безье, управляемая теми же опорными точками:
In[11]: =
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 [To St ring [# - 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, "Степень В-кривой" } ,
Dynamic@RadioButtonBar [Dynamic [d] , Range[0, 7]] } ,
{{curve, 1, "Кривая"} , {1 -> "В-кривая" , 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[l, {i, l,n}])]
Задача 6.1.7. Самостоятельно запрограммировать B -кривую с помощью встроенной функции BSplineBasis, так, чтобы можно было воспользоваться возможностью полноценного управления B -кривой. Обеспечить возможность произвольного допустимого порядка B -сплайнов.
Пусть даны опорные точки $$p_1, \dots , p_n$$ и их веса $$\omega_1, \dots , \omega_n,$$ а также расширенное множество неубывающих узлов $$\{t_1, \dots ,t_{n+m} : t_1 \le \dots \le t_{n+m}\}.$$ Опишем алгоритм вычисления B -кривой
Фиксируем t. Положим $$t_0 = -\infty, t_{n+m+1} = + \infty.$$ Тогда существует единственное $$i_0$$ такое, что $$t \in [t_{i_0},t_{i_0+1}),$$ где $$0 \le i_0 \le n + m.$$
Для определенного выше индекса $$i_0,$$ если $$1 \le i_0 \le n + m,$$ вычисляем единственное ненулевое значение ненормированного B -сплайна первого уровня (m=1):
Напомним, что $$i_0$$ было выбрано на первом шаге, так чтобы $$t \in [t_{i_0}, t_{i_0+1}).$$ При этом, в силу условия $$1 \le i_0 \le n + m - 1,$$ знаменатель $$t_{i_0+1} - t_{i_0} < \infty,$$ а в силу условия $$t \in [t_{i_0}, t_{i_0+1})$$ имеем $$t_{i_0+1} - t_{i_0} > 0.$$ Следовательно, $$0 < M_{i_0, 1}(t) > \infty.$$ Для индекса $$i_0,$$ находящегося вне множества $$\{1, \dots , n+m-1\},$$ значение $$M_{i_0,1}$$ не вычисляется.
С помощью соотношений Кокса - де Бура вычисляем все отличные от нуля в точке $$t$$ ненормированные сплайны $$m$$ -го порядка $$M_{j,m}(t) $$ при $$1 \le j \le n$$:
$$M_{j, m}(t)= \frac{(t_{j+m}-t)M_{j+1, m-1}(t)+(t-t_j)M_{j, m-1}(t)}{t_{j+m}-t_j}$$В частности,
$$M_{i_0, s}(t)=\frac{(t-t_{i_0})^{s-1}}{(t_{i_0+1}-t_{i_0}) \dots (t_{i_0+s}-t_{i_0})},\\ N_{i_o, s}(t)=\frac{(t-t_{i_0})^{s-1}}{(t_{i_0+1}-t_{i_0}) \dots (t_{i_0+s-1}-t_{i_0})},$$где $$1 \le s \le m, t \in [t_{i_0}, t_{i_0+1}).$$ Положив в (6.6) $$i_0 = 1,$$ получим
$$N_{1,s}(t)=\frac{(t-t_1)^{s-1}}{(t_2-t_1) \dots (t_s-t_1)},\\ 1 \le s \le m,\\ t \in [t_1, t_2).$$Лемма 6.5. Имеет место соотношение:
$$N_{1, m}(t)= \begin{cases} 1, \qquad \mbox{если } t_1= \dots = t_m,\\ 0, \qquad \mbox{иначе}. \end{cases}$$Вычисляем нормированные сплайны $$N_{j, m}(t) $$ для каждого $$1 \le j \le n$$ по формуле $$N_{j, m}(t) = (t_{j+m} - t_j)M_{j, m} (t) $$ (при этом $$1 \le j,j + m \le n + m$$ ).
Окончательно вычисляем $$r(t) $$ по формуле
$$r(t)=\frac{\sum_{i=1}^nN_{i, m}(t) \omega_i p_i}{\sum_{i=1}^n N_{i, m}(t) \omega_i}$$Теорема 6.4. r(t) B -кривой:
может быть вычислен с помощью следующего алгоритма:
Для заданного $$t \in [t_{i_0},t_{i_0+1}) (m \le i_0 \le n) $$ вычисляем величины $$r_i^{(k)}, \omega_i^{(k)} $$ при $$k = 1, \dots , m - 1, i = i_0 - m + k + 1, \dots , i_0$$ по рекуррентным формулам:
$$r_i^{(k)}=\frac{t_{i+m-k}-t}{t_{i+m-k}-t_i}r_{i-1}^{(k-1)}+\frac{t-t_i}{t_{i+m-k}-t_i}r_i^{(k-1)},\\ \omega_i^{(k)}=\frac{t_{i+m-k}-t}{t_{i+mk}-t_i}\omega_{i-1}^{(k-1)}+\frac{t-t_i}{t_{i+m-k}-t_i} \omega_i^{(k-1)},\\ r_o^{(0)}= \omega_i p_i, \omega_i^{(0)}=\omega_i, i=i_0-m+1, \dots, i_0 $$Вычисляем r(t) по формуле
Данный алгоритм иллюстрируется следующей диаграммой ( $$t \in [t_{i_0}, t_{i_0+1})$$ ):
$$\begin{matrix} r_{i_0-m+1}^{(0)} r_{i_0-m+2}^{(0)} \dots r_{i_0-1}^{(0)} r_{i_0}^{(0)}\\ r_{i_0-m+2}^{(1)} \dots r_{i_0-1}^{(1)}r_{i_0}^{(1)}\\ \ddots \vdots \vdots \\ r_{i_0-1}^{(m-2)} r_{i_0}^{(m-2)}\\ r_{i_0}^{(m-1)} \end{matrix}$$Определение 6.2.1. Пусть дано (m +1) * (n+1) точек в пространстве $$R^3,$$ образующих прямоугольную матрицу (сетку Безье):
Поверхностью Безье порядка m * n, соответствующей сетке $$P_{mn},$$ называется поверхность
где E - оператор сдвига вперед по первому индексу, F - оператор сдвига вперед по второму индексу:
Операторы E и F очевидно коммутируют друг с другом: $$EFp_{ij} = FEp_{ij} = p_{(i+1)(j+1)} .$$ Поэтому $$(1-v + Fv)^n(1 -u + Eu)^mp_{00} = (1 - u + Eu )^m(1 -v + Fv)^np_{00},$$ т. е. формула (6.7) непротиворечива.
Определение 6.2.2. Рациональная поверхность Безье, построенная по точкам $$p_{ij}$$ с весами $$\omega_{ij}$$ определяется следующим образом:
$$S_{00}(u, v; 1,1)=\frac{(1-u+Eu)(1-v+Fv) \omega_{00}p_{00}}{(1-u+Eu)(1-v+Fv) \omega_{00}}$$Формула (6.8) означает, что мы строим поверхность Безье в $$R_+^{k+1}$$ по точкам $$(\omega_{ij}; р_{ij} ,\omega_{ij}) \in R_+^{k+1},$$ а затем применяем преобразование $$\varphi R_+^{k+1} \to R^ k , \varphi (x, \omega) =\frac{x}{\omega}.$$
Задача 6.2.1. Наглядно изучить влияние опорных точек и их весов на рациональную поверхность Безье, меняя исходные данные pts (опорные точки) и $$\omega$$ (их веса) в нижеследующей программе. В Mathematica применять функцию BezierFunction к матрицам из скаляров.
Пример 6.2.1. Рациональная поверхность Безье с возможностью непосредственного управления весами с помощью движков.
In[12]: =
DynamicModule [ {pts, a, w0, pw, w, g, f, w6, w7, wl0, wll, 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}, {wl0, 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]]}],
ParametrioPlot3D[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},
{{wl0, 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])]
Поверхность Безье можно получить следующим образом:
Строим (n + 1) кривую Безье по столбцам матрицы $$P_{mn}$$:
u.При этом соответствующий оператор перехода от опорной точки, взятой на кривой $$r_j(u) ,$$ к следующей опорной точке, взятой на кривой $$r_{j+1}(u) ,$$ в терминах узлов сетки pij описывается оператором сдвига F по второму индексу. Следовательно,
$$S(u,v) = (1-v + Fv)^nr_0(u) = (1-v + Fv)^n (1-u + Eu )^m p_{00} = S_{00}(u,v;m,n),$$где m и n - количество шагов по первому и второму индексу соответственно, начиная от угловой точки сетки Безье, 00 - индекс этой угловой точки, с которой мы начинаем образовывать все остальные узлы сетки операторами сдвига E и F. Имеем
Здесь через $$S_{ij}$$ обозначены четырехугольные листы Безье, построенные из точек p ij, взятых в качестве угловых точек, аналогично формуле (6.9).
Аналогично формуле (6.10), любая точка на поверхности Безье степени m*n может быть представлена в виде точки четырехугольного листа Безье с углами в точках соответственно подобранных поверхностей Безье порядка (m-1)*(n-1):
Поверхность Безье можно представить в следующим виде:
$$S_{00}(u,v;m,n) = F(u,m)PF^{\top}(v,n),$$где $$P = (p_{ij})$$ и
$$F(u,m) = (f_0(u,m), \dots ,f_m(u,m) ) ,\\ F(v,n)= ( f_0(v,n), \dots ,f_n(v,n) ) ,\\ f_i ( \xi,k) = C_k^i(1- \xi)^{k-i} \xi^i.$$Аналогично обычной поверхности Безье, рациональная поверхность Безье $$S_{00}(u, u; 1,1) $$ может быть получена с помощью следующего алгоритма:
Полученная точка и будет точкой с параметрами $$(u_0,v_0) $$ на рациональной поверхности Безье, построенной по $$(p_{ij},\omega_{ij})_{i,j=0}^1.$$
Деление поверхностей Безье целиком основано на делении кривых Безье. Как мы уже знаем, каждая из кривых Безье (обычных или рациональных), построенных по столбцам сетки опорных точек $$P_{mn},$$ может быть разделена в параметрической точке $$u = u^* \in [0,1] $$ без изменения ее формы. При таком делении получаются две поверхности Безье (составляющие вкупе исходную поверхность) и соответствующие им две матрицы опорных точек Безье, каждая из которых имеет тот же порядок, что и исходная матрица: $$P_{тп}^a$$ и $$P_{тп}^b.$$ Затем, каждую из получившихся двух поверхностей Безье надо поделить в параметрической точке $$v = v^*$$ (с помощью деления кривых Безье, построенных по строкам соответствующих матриц). В итоге получится четыре поверхности Безье и соответствующие им четыре матрицы опорных точек: $$P_{mn}^{ac}, P_{mn}^{ad}, P_{mn}^{bc}$$ и $$P_{mn}^{bd}.$$ Эти четыре поверхности Безье будут лежать на исходной поверхности Безье, делить ее на четыр е части и в совокупности составлять исходную поверхность.
Новые сети опорных точек $$P_{mn}^{ac}, P_{mn}^{ad}, P_{m n}^{bc}$$ и $$P_{mn}^{bd}$$ получаются из исходной сети $$P_{mn}$$ по следующим формулам, аналогичным формулам (5.26):
$$p_{ij}^{ac} = (1-u^*+Eu^*)^i(1-v^*+Fv^*)^jp_{00}\\ p_{(m-i)j}^{bc} = ( (1 - u^*)E^{-1} + u^* )^i (1 -v^* + Fv^*)^jp_{m0},\\ p_{i(n-j)}^{ad} = (1-u^*+ Eu^*)^i ( (1 - v^*)F^{-1} + v^* )^j p_{0n},\\ p_{(m -i)(n-j)}^{bd} = ( (1 - u^*)E^{-1} + u^* )^i ( (1 - v^*)F^{-1} + v^* )^j p_{mn}.$$Добавление новых узлов в сетку Безье без изменения самой поверхности происходит по формулам, полностью аналогичным формулам добавления опорных точек в кривую Безье. Если дана сетка Безье $$P_{mn} = (p_{ij})_{i,j}^{m,n}=0,$$ то увеличенная на одну строку и на один столбец сетка $$Q_{(m+1)(n+1)} $$ для той же самой поверхности Безье строится по формулам
$$q_{ij}=\frac{1}{m+1}\frac{1}{n+1}((m+1-i)(1 - \delta_{(m+1)i})+i(1-\delta_{0i})E^{-1})\times\\ \times((n+1-j)(1-\delta_{(n+1)j})+j(1-\delta_{0j}F^{-1})p_{ij}.$$При этом углы новой сетки совпадут с углами старой: $$q_{00} = p_{00}, q_{(m+1)0} = p_{m0}, q_{0(n+1)} = p_{0n}, q_{(m+1)(n+1)} = p_{mn}.$$
Прежде всего заметим, что с помощью деления поверхности Безье любая наперед заданная точка этой поверхности может быть представлена как угловая точка новой поверхности Безье, являющейся частью исходной, а потому имеющей в каждой своей точке те же самые дифференциально-геометрические свойства, что и исходная поверхность Безье. Следовательно, для вычисления дифференциально-геометрических характеристик поверхности Безье в произвольной точке достаточно уметь вычислять эти характеристики лишь для угловой точки поверхности Безье с нулевыми значениями параметрических координат.
Таким образом, достаточны следующие формулы.
Касательные векторы к поверхности Безье в угловой точке $$p_{00}$$ могут быть вычислены по формулам $$T_u(0) = m(p_{10} - p_{00})$$ и $$T_v(0) = n(p_{01} - p_{00}).$$
Нормаль к поверхности в угловой точке дается формулой $$N =\frac{T_u(0) \times T_v(0)}{|T_u(0) \times T_v(0)|} $$
Вторые производные
Вычислив, согласно приведенным формулам, первые и вторые производные в угловой точке, а также нормаль к поверхности в этой точке, мы сможем вычислить гауссову и среднюю кривизны, главные направления и главные кривизны в угловой точке. Все эти величины будут выражаться через узлы сетки: $$p_{00}, p_{01}, p_{02}, p_{10}, p_{20}.$$
Для произвольной точки $$S_{00}(u^*,v^*;m,n) $$ поступаем так:
Пример 6.2.2. Чайник, нарисованный с помощью поверхности Безье. Чайник состоит из заплат:
In[13]:=
DynamicModule[{potPoints, рр2, potPatches, t, i, n, Bernstein,
BezierSurface}, potPoints = {{1.4, 0, 2.4},
(* Оставшиеся 304 точки *)
{1.425, -0.798, 0}}; (* это точки каркаса чайника, их 306 *)
рр2 = With[{avg = Mean [potPoints]}, Map[(j" - avg) , potPoints]];
(*усреднили все точки и составили из них список*)
potPatches = {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16},
(* оставшиеся тридцать заплат *)
{270, 270, 270, 270, 300, 305, 306, 279, 297, 303, 304, 275, 294,
301, 302, 271}}; ("чайник сосотоит из заплат, их 32 штуки*)
Bernstein [i_Tnteger, 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_, y_, 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, "Разбиение"}, 0, 1},
SaveDefinitions -> True]]
Пусть задана матрица опорных точек $$(p_{ij} ), 1 \le i \le n, 1 \le j \le m,$$ с весами ( $$\omega_{ij}$$ ). Фиксируем числа $$p$$ и $$q$$ - порядки B -сплайнов в двух координатных направлениях u и v соответственно. Тогда B -поверхность порядка (p,q), построенная по матрице опорных точек $$(p_{ij}) $$ с весами $$(\omega_{ij}) $$ и по двум расширенным множествам узлов в направлениях u и v соответственно, определяется формулой
Если ввести новые обозначения для входящих в эту формулу произведений B -сплайнов
то формула (6.11) перепишется в виде
$$r(u,v)=\frac{\sum_{i=1}^n \sum_{j=1}^m N_{ij,qp}(v,u) \omega_{ij}p_{ij}}{\sum_{i=1}^n \sum_{j=1}^m N_{ij,qp}(v,u) \omega_{ij}}$$Если расширенные множества узлов в u и v -направлениях имеют вид
то область изменения параметров (u,v) на B -поверхности (6.12) будет следующей:
Задача 6.2.2. Наглядно изучить влияние опорных точек и расширенного множества узлов на B-поверхность, порождаемую пакетом Mathematica с помощью следующей программы. Поменять в программе автоматический выбор расширенных множеств узлов на прямое задание этих узлов. Изучить влияние расширенных множеств узлов на B -поверхность.
Пример 6.2.3. B -поверхность. Для обеих поверхностей задан одинаковый остов. В первом случае поверхность "замыкается" только по оси Ox, во втором случаем - по обеим осям.
In[14] : = 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 = BSplineEunction [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}]]
Пример 6.2.4. Моделирование кузова автомобиля с помощью В -поверхности.
In[19]: =
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, 140}}], 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}}
(* Другие опорные точки *)
{{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}}
})]
Пример 6.2.5. Интерполяция заданного рисунка с помощью поверхности Безье и B -поверхности.
In[20]: =
DynamicModule {f, colfun, prevn =5, prevapprox =True, start, locators, grid,
order, approx, images},
Manipulate {* Создали список, который задает рисунок,
и для этого списка построили В-функцию, т.е. мы "синтерполировали" рисунок*)
colfun = BSplineFunction[ImageDatafimages[[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"1}] ,
f = Interpolation[
Flatten[Table[{{N[1 / (n-1)] i, N[1/ (n-1)] j } , pts [ [nit] + 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, {}]},
ImageSise -> 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 -> "В-поверхность" } ,
ControlType -> RadioButton} ,
{{order, 3, "Порядок В-поверхности"},
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, l/(n-l)}, {j, 0, 1, l/(n-l)}], 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 = {
};) ]]
В ряде приложений (например, в проектировании корпусов и крыльев самолетов и кораблей) требуется строить гладкие поверхности, затягивающие заданные остовы из кривых - криволинейные сетки, криволинейные направляющие (например, несущий остов корпуса корабля), различные граничные контуры. Во многих случаях важное значение имеет также и вычислительная сторона вопроса - сложность вычисления
В компьютерной геометрии существует несколько различных подходов к решению подобных задач. Ниже будет рассмотрены некоторые из этих подходов, ставшие классическими.
Линейчатая поверхность - геометрическое место точек отрезков прямых, соединяющих соответствующие друг другу точки двух заданных кривых. Пусть $$\rho_1(t), t_0 \le t \le t_1,$$ и $$\rho_2(w), w_0 \le w \le w_1,$$- две заданные кривые. Тогда определяемая ими линейчатая поверхность имеет вид
$$r(u,v) = (1-v)\rho_1(t) + v \rho_2(w),$$где $$t = t_0(1-u)+t_1u, w = w_0(1 - u) + w_1u, 0 \le u, v \le 1.$$
Если обе кривые $$\rho_1$$ и $$\rho_2$$ замкнуты, то линейчатая поверхность замкнута по u. По параметру v линейчатая поверхность всегда незамкнута.
Если обе кривые $$\rho_1$$ и $$\rho_2$$ - отрезки прямых, то линейчатая поверхность называется билинейной и определяется только своими четырьмя угловыми точками $$p_1, p_2, p_3, p_4$$:
$$r(u, v) = (1 - v) ( (1 - u)p_1 + up_2 ) + v ( (1 - u)p_3 + u p_4 ) =\\ = (1 - u)(1 - v)p_1 + u(1 - v)p_2 + (1 - u)vp_3 + uvp_4,\\ 0 \le u, v \le 1.$$Секториальная поверхность - частный случай линейчатой поверхности, когда одна из кривых $$\rho_1$$ и $$\rho_2$$ вырождается в точку:
$$r(u,v) = (1-v) \rho (u)+vp,$$где $$\rho ( u) $$ - кривая, p - точка, $$u_1 \le u \le u_2, 0 \le v \le 1.$$
Секториальная билинейная поверхность называется треугольной поверхностью.
Линейные поверхности Кунса. Вычтем и добавим к уравнению линейчатой поверхности уравнение билинейной поверхности, построенной по четырем угловым точкам (см. выше). Уравнение линейчатой поверхности примет вид
$$r(u,v) = (1-v)\rho_1(u)+v \rho_2(u) = (1-v) \rho_1(u)+v \rho_2(u)+(1-u) ( (1-v)p_1+vp_3 ) +\\ + u ( (1 - v)p_2 + vp_4 ) - (1 - u)(1 - v)p_1 - u(1 - v)p_2 - (1 - u)vp_3 - uvp_4 =\\ = (1 - v) \rho_1(u) + v \rho_2(u) + (1 - u) \rho_3(v) + u \rho_4(v) - (1 - u)(1 - v)p_1-\\ - u(1 - v)p_2 - (1 - u)vp_3 - uvp_4,\\ 0 \le u, v \le 1,$$где $$\rho_1 (u) $$ и $$\rho_2 (u) $$ - кривые, по которым строится данная линейчатая поверхность, причем считается, что на них обеих задан один и тот же параметр $$0\leu\le1,$$и
$$\rho_3(v) = (1-v)p_1+vp_3,\\ \rho_4(v) = (1 - v)p_2 + vp_4$$- отрезки прямых.
Подставим теперь в уравнение (6.13) вместо отрезков прямых (6.14) уравнения двух произвольных кривых $$\rho_3(v) $$ и $$\rho_4(v) ,$$ соединяющих соответственно пары точек $$(p_1, p_3) $$ и $$(p_2, p_4) .$$ Тогда новая поверхность, задаваемая уравнением (6.13), будет определяться четырьмя кривыми $$\rho_1, \rho_2, \rho_3, \rho_4,$$ образующими криволинейный четырехугольник, натянутый на угловые точки $$p_1, p_2, p_3, p_4.$$ При этом все четыре кривые $$\rho_1, \rho_2, \rho_3, \rho_4$$ войдут в уравнение (6.13) симметрично. Получится линейная поверхность Кунса.
Определение 6.3.1. Пусть заданы четыре кривые, образующие криволинейный четырехугольник с углами в точках $$p_1, p_2, p_3, p_4$$:
$$\begin{martix} \gamma_a=r(u,0),0 \le u \le 1,r(0,0)=p_1,r(1,0)=p_2,\\ \gamma_b=r(0,v) 0 \le v \le 1, r(0,0)=p_1 r(1,0)=p_3),\\ \gamma_c(1,v) 0 \le v 1 r(0,0)=p_1 r(0,1)=p_3,\\ \gamma_d=r(u,1),0 \le u \le 1, r(0,1)=p_3, r(1,1)=p_4 \end{matrix} $$Обозначим
$$\alpha_0(t) = 1-t,\\ \alpha_1(t)=t.$$Функции $$\alpha_0(t), \alpha_1(t) $$ будем называть функциями смещения линейной поверхности Кунса. Сама линейная поверхность Кунса определяется следующей формулой:
$$r(u, v) = \alpha_0(v)r(u, 0) + \alpha_0(u)r(0, v) + \alpha_1(v)r(u, 1)+\\ + \alpha_1(u)r(1,v) -\alpha_0(u) \alpha_0(v)r(0,0) -\alpha_1(u)\alpha_0(v)r(1,0)-\\ - \alpha_0(u)\alpha_1(v)r(0,1) -\alpha_1(u)\alpha_1(v)r(1, 1).$$Замечание 6.3.1. Линейная поверхность Кунса не является, вообще говоря, линейчатой поверхностью. Название "линейная" обусловлено тем, что в ее определении используются линейные функции смещения $$\alpha_0$$ и $$\alpha_1.$$
Замечание 6.3.2. Параметрическая область поверхности Кунса может иметь вид прямоугольника, не обязательно совпадающего с квадратом $$[0,1] \times [0,1] .$$
Матричный вид уравнения поверхности Кунса. Преобразуем формулу (6.15) к другому виду. Имеем в матричных обозначениях
$$r(u,v) = ( \alpha_0(u) \alpha_1(u) ) \begin{pmatrix} r(0,v)\\ r(1,v) \end{pmatrix} + ( \alpha_0(v) \alpha_1(v) ) \begin{pmatrix} r(u,0)\\ r(u,1) \end{pmatrix} -\\ -( \alpha_0(u) \alpha_1(u)) \begin{pmatrix} r(0,0) r(0,1) r(1,0) r(1,1) \end{pmatrix} \begin {pmatrix} \alpha_0(v)\\ \alpha_1(v) \end{pmatrix} =\\ = \sum_{i=0}^1 \alpha_i(u)r(i, v) + \sum_{j=0}^1 \alpha_j(v)r(u,j) - \sum_{i=0}^1 \sum_{j-0}^1 \alpha_i(u)\alpha_j(v)r(i, j).$$Пример 6.3.1. Зададим с помощью поверхности Кунса область D плоскости, ограниченную эллипсом $$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$$ или $$\begin{cases}
x=a\ cos\ \varphi\\
y=b\ cos\ \varphi,
\end{cases}$$ где $$0 \le \varphi < 2 \pi,$$ с полуосями $$a > 0, b > 0.$$ Представим границу области D в виде четырех кривых - дуг эллипса, стыкующихся в точках $$p_1, p_2, p_3, p_4.$$ Соответствующая линейная поверхность Кунса имеет вид
Преимущество задания плоских областей в виде линейных поверхностей Кунса состоит в том, что у поверхности Кунса параметрическая область -всегда прямоугольник.
Пример 6.3.2. Поверхность Кунса, построенная по кривым, лежащим на заданных поверхностях. Кривые раскрашены красным, зеленым, желтым и оранжевым.
In[21]: =
DynamicModule[{rl = 30 . , r2 = 30. , r3 = 30 . , d = 50 . , u , v, ml, m2 , m3, m4 , fl,
f2, f3, f4, rfel, α, MV, MU, laml, lam2, rfk}, rl = 30. ; r2 = 30. ; r3 = 30. ;
d = 50. ; laml = 2. ; lam2 = 1. ;
f1 [ml_] : = {rl - rl *Cos[u * Pi] , laml * 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, m1* (Sin[v* 3* Pi] )^2} :
f4 [m4_] :={1.6*r1+0.4*rl*u+(v^4-2*v^3+v^2)*100*m4*u, d*v, 0} ;
Manipulate[
Show[{ParametricPlot3D[Evaluate@fl [ml] , {u, 0, 1}, {v, -1, 0}, Mesh -> {24, 24},
PlotRange -> All] ,
ParametricPlot3D[Evaluate@{fl [ml] /. 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 -> l}, {u, 0, 1}, Mesh -> {24, 24},
PlotRange -> All, PlotStyle -> {Green} ] ,
ParametricPlot3D[Evaiuate@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 -> l}, {v, 0, 1}, Mesh -> {24, 24},
PlotRange -> All, PlotStyle -> {Orange} ] ,
ParametricPlot3D [Evaluate@rfel [ml, 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}}
], {{ml, 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 : -> (
rl = 30 . ; r2 = 30. ; r3 = 30. ; d = 50. ; laml = 2. , lam2 = 1. ;
fl [ml_] : = {rl - rl *Cos [u * Pi] , laml *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*rl+0.4*rl*u+(v^4-2*v^3+v^2)*100*m4*u, d*v, 0} ;
\alpha [t_] := t; MV[v_] : = {1 - \alpha [v] , \alpha [v] } ; MU[u_]: = {1 - \alpha [u] , \alpha [u] } ;
(* наши функции а *)
Table [rfk [ml, m2, m3, m4] [i] , {i, 3}] ;
Do[rfk[ml_, m2_, тЗ_, т4_] [i] =
MU[u].{f3[m3][[i]] /. u ->0, f4[m4][[i]] /. u ->l} +
MV[v] . {fl [ml] [ [i] ] /. v->0, f2[m2][[i]] /. V->1} -
MU[u] . {{f3[m3] [[i]] /. {u ->0, v-> 0} , f3[m3] [[i]] /. {u -> 0, v->l}},
{f4 [m4] [ [i] ] /. {u -> 1 , v-> 0} , f4 [m4] [ [i] ] /. {u-> 1, v -> 1} } } .MV[v] //
Together;
, {i, 3}] ; rfel[ml_, m2_, тЗ_, m4_] =
{rfk [ml, m2, m3, m4] [1] , rfk [ml, m2, m3, m4] [2] , rfk [ml, m2, m3, m4] [3] };) ]]
Обобщенные поверхности Кунса. При стыковке линейных поверхностей Кунса по граничным кривым производные в направлении, ортогональном границе, претерпевают разрыв. Чтобы гладко склеивать поверхности Кунса, необходимо строить их таким образом, чтобы они имели на границе заданные (согласованные) производные до k -го порядка включительно. Производные в касательном направлении определяются уравнениями самих граничных кривых, а производные в трансверсальном к границе направлении необходимо задавать дополнительно.
Определение 6.3.2. Обобщенной поверхностью Кунса n-го порядка называется поверхность
$$r(u,v)=\sum_{i=0}^1 \sum_{k=0}^n \alpha_{k,i}(u)r_u^{(k)}(I,v)+\sum_{j=0}^1 \sum_{k=0}^n \alpha_{k,j}(v)r_v^{(k)}(u,j)-\\ -\sum_{i=0}^1 \sum_{j=0}^1 \sum_{k=0}^n \sum_{l=0}^n \alpha_{k,i}(u) \alpha_{l,j}(v)r_{uv}^{(k,l)}(I,j),\\ 0 \le u,v \le 1$$задаваемая четырьмя граничными кривыми $$r(u,0), r(0,v), r(1,v), r(u, 1) ,$$ стыкующимися в точках $$p_1 = r(0, 0), p_2 = r(1,0), p_3 = r(0,1) $$ и $$p_4 = r(1,1) ,$$ а также чистыми частными производными в трансверсальном направлении на граничных кривых до $$n$$ -го порядка включительно и смешанными частными производными в угловых точках $$p_1, p_2, p_3, p_4$$ до $$2n$$ -го порядка включительно:
$$\begin{matrix} r_u^{(k)}(I,v)=\frac{\partial^k r(u,v)}{ \partial u^k} \left |_{u=i}, i=0,1, k=1, \dots, n\\ r_v^{(k)}(u,j)=\frac{\partial^kr(u,v)}{ \partial v^k} \left |_{v=j}, j=0,1, k=1, \dots, n\\ r_{uv}^{(k,l)}=\frac{\partial^{k+l}(u,v)}{ \partial u^k \partial v^l} \left|_{u=i, v=j} i,j=0,1, k,l=1, \dots, n \end{matrix} $$а также обобщенными функциями смещения $$\alpha_{k ,i}(t) ,$$ где $$0 \le k \le n, i = 0, 1,$$ в качестве которых могут быть взяты произвольные гладкие функции, удовлетворяющие соотношениям
$$\frac{d^m \alpha_{k,i}(t)}{dt^m} \left|_{t=j}=\delta_i^j \delta_k^m=\begin{cases} 1, \mbox{если} i=j, k=m,\\ 0, \mbox{иначе}. \end{cases}$$Здесь $$0\le k, m \le n.$$
Задача 6.3.1. Проверить, что поверхность (6.16), определенная с помощью (6.17) и (6.18), имеет заданные в (6.17) производные на границе.
Замечание 6.3.3. В качестве обобщенных функций смещения всегда можно выбрать многочлены.
Пример 6.3.3. Определим функции смещения формулами
$$\alpha_0(t) = 1-3t^2 + 2t^3, \alpha_1(t) = 3t^2-2t^3,\\ \beta_0(t)=t-2t^2 +t^3 , \beta_1(t) = -t^2 +t^3$$(здесь для простоты использованы обозначения $$\alpha_{0,0}(t) = \alpha_0(t), \alpha_{0,1}(t) = \alpha_1(t), \alpha_{1,0}(t) = \beta_0(t), \alpha_{1, 1}(t) = \beta_ 1(t)) .$$ Пусть заданы , четыре кривые $$r(u, i), r(i,v) ,$$ где $$i = 0, 1,$$ стыкующиеся в точках $$r(0,0), r(0,1), r(1,0), r(1,1) .$$ Пусть заданы первые производные на границе в трансверсальном направлении: $$r_u(i, v), r_v(u, i), i = 0, 1,$$ и вторые смешанные производные в угловых точках: $$r_{uv}(i,j),i,j = 0,1.$$
Кубическая поверхность Кунса определяется формулой (в матричных обозначениях)
$$r(u,v)=(\alpha_0(u) \alpha_1(u) \beta_0(u) \baeta_1(u)) \left( \begin{matrix} r(0,v)\\ r(1,v)\\ r_u(0,v)\\ r_u(1,v) \end{matrix} \right)+\\ +(\alpha_0(v) \alpha_1(v) \beta_0(v) \beta_1(v) ) \left( \begin{matrix} r(u,0)\\ r(u,1)\\ r_v(u,0)\\ r_v(u,1) \end{matrix} \right)-\\ - \left( \begin{matrix} \alpha_0(u)\\ \alpha_1(u)\\ \beta_0(u)\\ \beta_1(u) \end{matrix} \right )^{\top} \begin{pmatrix} r(0,0) r(0,1) r_v(0,0) r_v(0,1)\\ r(1,0) r(1,1) r_v(1,0) r_v(1,1)\\ r_u(0,0) r_u(0,1) r_{uv}(0,0) r_{uv}(0,1)\\ r_u(1,0) r_u(1,1) r_{uv}(1,0) r_{uv}(1,1) \end{pmatrix} \begin{pmatrix} \alpha_0(v)\\ \alpha_1(v)\\ \beta_0(v)\\ \beta_1(v) \end{pmatrix} ,\\ 0 \le u,v \le 1$$Задача 6.3.2. Написать на пакете Mathematica программу построения кубической поверхности Кунса.
В проектировании поверхностей иногда применяется метод построения поверхности по дискретной совокупности ее сечений. Для затягивания остова, образованного сечениями, гладкой поверхностью применяются методы, аналогичные построению интерполяционных кривых. Поэтому получаемые таким способом поверхности имеют такие же названия, как и соответствующие интерполяционные сплайны.
Поверхности Эрмита. Напомним, что сплайном Эрмита называется кусочно-кубическая кривая, проходящая через заданные точки $$p_0, \dots , p n \in R^N$$ и имеющая в этих точках заданные производные $$q_0, \dots , q_n \in R^N.$$ Если производные $$q_i$$ не заданы, то в точках $$p_1, \dots , p_{ n-1}$$ их можно определить из некоторых "естественных" соображений. Одним из примеров такого рода является псевдоупругий сплайн Эрмита. Еще один способ - положить по определению
$$q_i=\frac{p_{i+1}-p_i}{t_{i+1}-t_i}$$(есть и другие способы).
Задание производных на концах $$p_0$$ и $$p_n$$ иногда заменяют на условие обращения в ноль третьей производной на концах кривой. Это условие равносильно следующим формулам:
$$\begin{cases} q_0=2(p_1-p_o)-q_1,\\ q_n=2(p_n-p_{n-1})-q_{n-1} \end{cases} $$Напомним, что вторые производные в опорных точках сплайна Эрмита (где происходит стык кубических кривых) претерпевают разрыв, т. е. сплайны Эрмита принадлежат классу $$C^1.$$ Аналогично сплайнам Эрмита строятся поверхности Эрмита.
Пусть имеется семейство кривых
$$\begin{cases} r=r_0(u)=r(u,0),\\ ..................\\ r=r_i(u)=r(u,i),\\ ..................\\ r=r_n(u)=r(u,n) \end{cases} $$Без ограничения общности будем считать, что параметр u для всех кривых пробегает один и тот же отрезок: $$u_0 \le u \le u_1.$$ Тогда поверхность Эрмита, затягивающая остов кривых $$r(u,i), 0 \le i \le n,$$ определяется формулой, аналогичной формуле для сплайна Эрмита:
если $$j \le v \le j + 1.$$ Здесь $$w = v - j$$ - местный параметр на данном участке поверхности ( $$0 \le w \le 1$$ ), $$r_v(u,j) $$ и $$r_v(u,j + 1) $$ - производные в трансверсальном направлении на заданном остове кривых (6.20), $$\alpha_0(w) = 1-3w^2 + 2w^3, \alpha_1(w) = 3w^2-2w^3, \beta_0( w ) = w-2w^2+w^3, \beta_1(w) = -w^2+w^3.$$ Параметрическая область поверхности Эрмита имеет вид
$$\begin{cases} u_0 \le u \le u_1,\\ 0 \le v \le n \end{cases}$$Если трансверсальные производные $$r_v ( u, i) $$ на линиях остова не заданы, то их, аналогично случаю сплайна Эрмита, можно определить, например, так:
$$\begin{cases} r_v(u,j)=\frac{r(u,j+1)-r(u,j-1)}{2}, 1 \le j \le n-1,\\ r_v(u,0)=2r(u,1)-2r(u,0)-r_v(u,1),\\ r_v(u,n)=2r(u,n)-2r(u,n-1)-r_v(u,n-1). \end{cases} $$При таком определении третьи производные на кривых $$r_0$$ и $$r_n$$ в трансверсальном направлении обратятся в ноль:
$$\frac{\partial^3r(u,v)}{ \partial v^3} \left|_{v=0}=\frac{\partial^3r(u,v)}{ \partial v^3} \left|_{v=n}=0$$Сами исходные кривые $$r_0, \dots ,r_n,$$ образующие остов поверхности Эрмита, могут быть построены как кривые Безье или B -кривые. При этом получаются поверхности класса $$C^1,$$ имеющие на двух участках границы управляемые производные по трансверсальному направлению (а именно, на участках $$r(u,0) $$ и $$r(u,n) $$ ). На участках границы $$r(u_0,v) $$ и $$r(u_1,v) $$ производные, вообще говоря, управляются лишь приблизительно - с помощью выбора кривых остова $$r_i ( u) .$$
Применение поверхностей Эрмита: поверхность перехода. Примером применения поверхностей Эрмита может служить задача объединения двух поверхностей в одну общую поверхность путем соединения их краев с помощью поверхности перехода.
Пусть даны две поверхности M и N, которые надо сопрячь, соединив поверхностью перехода край $$r = r_M(t,u_1^M) $$ первой поверхности и край $$r = r_N(t,u_0^N) $$ второй поверхности. Будем считать, что оба сопрягаемых края отнесены к одному и тому же параметру t, изменяющемуся в пределах $$t_0 \le t \le t_1.$$ Второй параметр на поверхностях M и N обозначим (в обоих случаях) через u. Пусть параметрические пространства поверхностей имеют вид
для поверхности $$M : r = r_M ( t, u ) $$ и
$$\begin{cases} t_0 \le t \le t_1,\\ u_0^N \le u \le u_1^N \end{cases}$$для поверхности $$N : r = r_N(t, u) .$$ Без ограничения общности будем считать, что $$u_1^M \le u_0^N.$$ Тогда поверхность перехода от M к N может быть описана формулой поверхности Эрмита:
где $$w =\frac{u-u_1^M}{u_0^N-u_1^M}$$ - местный параметр по направлению $$u, t_0 \le t \le t_1, 0 \le w \le 1.$$ Данная поверхность Эрмита гладко сопрягает поверхности M и N по заданным краям. Класс гладкости сопряжения $$- C^1.$$
Пример 6.3.4. Приклейка ручки с помощью поверхности перехода:
In[22] := Style [Module [{v, rl, r2, s1, s2 , drl, dr2 , dsl, ds2 , r, s},
r1[u_, v_] : = {-3 + (7 + 1.5 Cos [2 π u]) Sin [2 π v + π] , 1.5 Sin [2 π u] ,
{7 + 1.5 Cos [2 π u]) Cos [2 π v + π] } ; (*тор, u и v в пределах [0,1]*)
r2 [u_, y_] := {5 Cos [2 π u] * Sin [π (v - 1) ] , 5 Sin [2 π] * Sin [π (v - 1) ] ,
5 Cos [π {v - 1) ] } ; {*сфера, u в пределах [0,1], v в пределах [1,2]*)
drl[t_] :=D[rl[t, v], v] * (7 / 6 - (1/2-1/10)) /.v -> l/2-l/10;
dr2[t_] := D[r2[t, v] , v] * (7 / 6 - (1/2-1/10)) / . v -> 7 / 6;
r[t_, w_] :=rl[fc, 1/2-1/10] (l- 3 w2 + 2 w3) + r2[t, 7/6] (з w2 - 2 w3) +
drl[t] (w- 2 w2 + w3) + dr2[t] (-w2 + w3 ) ;
sl[u_, v_] :=r2[u, v + 2] ; (* u в пределах [0,1], v в пределах [-1,0]*)
s2[u_, v_] :=rl[u, v] ;
dsl[t_] := D[sl[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_] :=sl[t, -1/6] (l- 3 w2 + 2 w3) + s2[t, 1/10] (з w2 - 2 w3) +
dsl[t] (w- 2 w2 + w3) + ds2[fc] (-w2 + w3 ) ;
Manipulate[Show[{
PararaetricPlot3D[rl[u, v] , {u, 0, 1} , {v, 1/10, 1 /2 - 1 / 10} ,
PlotStyle -> FaceForm[Orange, Blue]],
PararaetricPlot3D[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, wl},
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
] , {{wl, 1} , 0.001, 1}, {{w2, 0.001}, 0.001, 1}, TrackedSymbols -> {wl, w2},
SaveDefinitions -> True, ControlPlacement -> Lef t] 1 , Magnification -> 2]
Поверхности Лагранжа. Аналогично интерполяционной кривой Лагранжа, можно построить поверхность Лагранжа, затягивающую заданный остов кривых (6.20). Эта поверхность описывается
где
$$L_i(v)=\frac{\Pi_{j=0, j \ne i}^{n}(v-j)}{\Pi_{j=0, j \ne i}^{n}(i-j)}$$Последняя формула получается из формулы для коэффициентов интерполяционного полинома Лагранжа подстановкой $$t_i = i (0 \le i \le n) .$$
Поверхность Лагранжа имеет те же недостатки, что и v поверхность будет бесконечно гладкой и представляет собой многочлен. Гладкость по параметру u зависит от гладкости кривых остова.
Поверхности Гордона. Пусть задан остов поверхности, состоящий из двух семейств кривых в обоих параметрических направлениях, которые образуют сетку - т. е. каждая кривая одного семейства пересекает в одной и только одной точке каждую кривую другого семейства. Различные кривые одного и того же семейства не пересекаются. Обозначим заданную таким образом сетку кривых через
$$r = r(u,0), \dots ,r = r(u,n) \mbox{(первое семейство)},\\ r = r(0,v), \dots ,r = r(m,v) \mbox{ (второе семейство)}.$$Здесь $$0 \le u \le m, 0 \le v \le n.$$ j -й кривой первого семейства с i -й кривой второго семейства равен $$r(i,j), 0 \le i \le m,$$
Определение 6.3.3. Поверхностью Гордона, затягивающей сетку кривых (6.21), называется поверхность
$$r(u,v)=\sum_{j=0}^nL_j^n(v)r(u,j)+ \sum_{i=0}^m L_i^m(u)r(I,v)-\\ - \sum_{i=0}^m \sum_{j=0}^nL_i^m(u)L_j^n(v)r(i,j),\\ 0 \le u \le m,\\ 0 \le v \le n,$$где
$$L_i^m(u)=\frac{\Pi_{r=0, r \ne 0}^m(u-r)}{\Pi_{r=0, r \ne 0}^m(i-r)},\\ L_j^n(v)=\frac{\Pi_{s=0, s \ne j}^n(v-s)}{\Pi_{s=0, s \ne 0}^n(j-s)}$$- коэффициенты Лагранжа.
Обобщенная поверхность Гордона задается формулой (6.22), в которой вместо коэффициентов Лагранжа (6.23) используются произвольные $$C^1$$ -функции $$\alpha_k(w) ,$$ которые в целых точках соответствующего промежутка удовлетворяют условию
$$\begin{cases} 1, \mbox{если } w=k,\\ 0, \mbox{если } w \ne k \end{cases}$$Например, в качестве $$\alpha_k(w) $$ можно взять функции смещения (6.19). Гладкость обобщенной поверхности Гордона зависит от гладкости используемых функций $$\alpha_k(w) .$$ Поверхности Гордона, построенные по коэффициентам Лагранжа, имеют гладкость $$C^{\infty},$$ но могут образовывать необоснованные изгибы.
Поверхности, затягивающие сетку кривых заплатами Кунса. Остов, состоящий из сетки кривых (6.21), можно затянуть поверхностью, состоящей из заплат, затягивающих каждую отдельную ячейку сетки и гладко сопрягающихся между собой. В качестве таких заплат могут выступать поверхности Кунса. Получается следующая поверхность:
$$r(u,v)=(\alpha_0(t) \alpha_1(t) \beta_0(t) \beta_1(t)) \left (\begin{matrix} r(i,v)\\ r(i+1,v)\\ r_u(i, v)\\ r_u(i_1,v) \end{matrix} \right )+\\ +(\alpha_0(w) \alpha_1(w) \beta_0(w) \beta_1(w)) \left (\begin{matrix} r(u,j)\\ r(u, j+1)\\ r_v(u,j)\\ r_v(u,j+1) \end{matrix} \right )-\\ (\alpha_0(t) \alpha_1(t) \beta_0(t) \beta_1(t)) \times\\ \times \begin{pmatrix} r(i,j) r(i,j+1) r_v(i,j) r_v(i,j+1)\\ r(i+1,j) r(i+1, j+1) r_v(i+1,j) r_v(i+1, j+1)\\ r_u(I,j) r_u(i,j+1) r_{uv}(i,j) r_{uv}(i,j+1)\\ r_u(i+1,i) r_u(i+1,j+1) r_{uv}(i+1,j) r_{uv}(i+1,j+1) \end{pmatrix} \times\\ \times \begin{pmatrix} \alpha_0(w)\\ \alpha_1(w)\\ \beta_0(w)\\ \beta_1(w) \end{pmatrix} $$$$t = u ? i, w = v ? j, 0 \le u \le m, 0 \le v \le n.$$ Здесь $$i = [u] $$ - целая часть $$u, j = [v] $$ - целая часть $$v,t = u-i =\{u\}$$ - дробная часть $$u, w = v - j = \{v\}$$ - дробная часть $$v, t,w \in [0,1) $$ - местные параметры в данной ячейке сетки (6.21). Согласно (6.21), области изменения параметров u и v таковы: $$0 \le u \le m, 0 \le v \le n.$$
Для вычисления
при $$1 \le u \le m - 1, 1 \le v \le n - 1$$ (т. е. для внутренних линий сетки) и по формулам
$$r_u(0, v) = 2r(1, v) - 2r(0, v) - r_u(1, v),\\ r_u(m, v) = 2r(m, v) - 2r(m - 1, v) - r_u(m - 1, v),\\ r_v(u,0) = 2r(u, 1) - 2r(u, 0) - r_v(u, 1),\\ r_v (u, n) = 2r(u, n) - 2r(u, n - 1) - r_v (u, n - 1)$$для крайних линий сетки.
Геометрически это означает, что производная $$r_u(i,v) $$ в произвольной точке $$r(i,v) $$ внутренней линии сетки полагается пропорциональной вектору, соединяющему точки $$r(i - 1,v) $$ и $$r(i + 1,v) $$ на соседних с ней линиях сетки, отвечающие тому же значению параметра v, что и данная точка. Аналогичное утверждение справедливо и для производной $$r_v(u,j) .$$
Так же как и в случае кубических сплайнов Эрмита, формулы (6.25) обеспечивают равенство нулю третьих производных поверхности на границе в трансверсальном к граничной кривой направлении.
Поверхности тензорного произведения. Рассмотрим обобщенную поверхность Гордона (6.22), построенную по сетке кривых, которые являются
Для того чтобы получить формулу поверхности тензорного произведения, поставим в формулу (6.22) выражения, определяющие u и v:
где, напомним:
$$\alpha_0(x) = 1 - 3x^2 + 2x^3,\\ \alpha_1(x) = 3x^2-2x^3,\\ \beta_0(x) = x-2x^2+x^3,\\ \beta^1(x) = -x^2 + x^3,$$где $$x \in [0,1] $$ - местный параметр. В формулах (6.26) через $$i $$ и $$j$$ обозначены целые части параметров u и v соответственно. Местные параметры $$t$$ и $$w$$ определяются как $$t = u - i, w = v - j.$$ Вторые две формулы в (6.26) получаются дифференцированием формулы (6.22). После подстановки формул (6.26) в (6.22) и несложных алгебраических преобразований получим следующее выражение для
Из формулы (6.27) хорошо видно, что поверхность тензорного произведения определяется величинами, заданными лишь в узлах сетки.
При триангуляции поверхностей иногда бывает удобно рассматривать криволинейные поверхности треугольной формы, имеющие в качестве параметрической области треугольник на плоскости. При этом удобно пользоваться не декартовыми, а барицентрическими параметрическими координатами.
Лемма 6.6. Пусть даны три точки $$A, B, C \in R^2,$$ не лежащие на одной прямой. Тогда для любой точки $$P \in R^2$$ существует единственная тройка чисел $$\{a, b, c\}$$ таких, что
$$\begin{cases} r_P= ar_A + br_B + cr_C, \\ a + b + c=1. \end{cases}$$Здесь через $$r_P, r_A, r_B, r_C$$ обозначены радиус-векторы точек $$P, A, B, C$$ соответственно.
Числа $$a, b, c$$ называются барицентрическими координатами точки P относительно треугольника $$\Delta ABC.$$
Следствие 6.5. Имеют место формулы
$$|a|=\frac{S_{PBC}}{S_{ABC}},\\ |b|=\frac{S_{APC}}{S_{ABC}},\\ |c|=\frac{S_{ABP}}{S_{ABC}} $$т. е. коэффициенты $$a, b, c$$ по модулю равны отношению площади треугольника $$\Delta ABC,$$ в котором соответствующая вершина заменена на $$P,$$ к площади исходного треугольника $$\Delta ABC.$$
Следствие 6.6. Все три барицентрические координаты точки положительны тогда и только тогда, когда точка принадлежит внутренности треугольника $$\Delta ABC.$$ Все три барицентрические координаты точки неотрицательны тогда и только тогда, когда точка принадлежит треугольнику $$\Delta ABC.$$
Билинейная треугольная поверхность. Пусть даны три точки $$P_1, P_2, P_3$$ в $$R^3.$$ В качестве области параметров возьмем треугольник $$\Delta ABC \in R^2.$$ Точка параметрического пространства описывается координатами $$0 \le a, b, c \le 1, a+b+c = 1.$$
Треугольная поверхность на трех кривых. Пусть даны три попарно пересекающиеся кривые $$r = r_1(u), r = r_2(v), r = r_3(w ), 0 \le u,v,w \le 1,$$ пересекающиеся в точках $$P_1 = r_2(0) = r_3(1), P_2 = r_1(1) = r_3(0), P_3 = r_1(0) = r_2(1) .$$ Определим треугольную поверхность формулой
$$r(a, b, c) = a (r_3(1 - b) + r_2(c) - Р_1) + b (r_1(1 - c) + r_3(a) - P_2 ) + c (r_2(1-a) + r_1(b)-P_3 ),$$где $$0 \le a, b, c \le 1, a + b + c = 1.$$ Легко проверить, что при $$c = 0$$ формула (6.28) описывает кривую $$r_3(a) (0 \le a \le 1) ,$$ при $$a = 0$$ - кривую $$r_1(b) (0 \le b \le 1) ,$$ при $$b = 0$$ - кривую $$r_2(c) (0 \le c \le 1) .$$ Например, если $$a = 0,$$ то $$c = 1 - b$$ и, согласно (6.28), имеем
$$r(0, b,1 - b) = b (r_1(b) + r_3(0) - P_2 ) + (1 - b) (r_2(1) + r_1(b) - P_3 ) =\\ = b (r_1(b) + P_2 -P_2 ) + (1 - b) ( P_3 + r_1(b) -P_3) = r_1(b).$$Пример 6.3.5. Поверхность, соединяющая три поверхности по кривым на их границах. При изменении поверхностей, имеющих вид "треугольника", построенная поверхность не изменяется:
In[23]: =
DynamicModule[{pe200, pe002, pel0l, te, tee, pe020, pell0, pe0ll, r2 , rf2 ,
rf1, u, v, w, s, t, рb020, рb002, pb0ll, tb, tbb, pb200, pbl0l, pbll0, rl,
pj200, pj020, pjll0, tj , tjj, pj002, pjl0l, pj0ll, гЗ, rf3, flperem, f2perem,
f3perem, rfel},
pe200 = {0, 10, 0} ; pe002 = {10, 0, 0} ; pel01 = {10, 10, 0}; pb020 = {0, 0, 20};
pb002 = {10, 0, 0} ; pb0ll = {10, 0, 20} ; pj200 = {0, 10, 0}; pj020 = {0, 0, 20};
pjll0 = {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] ;
pell0[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 pel0l u w + 2 pel 10 [ te] u v +
2 pe011[te] vw;
rf2[ te_, s_, t_] := r2 [te] /.w ->l-u-v/. {u -> t- s t, v-> s t} ;
f2perem = ParametricPlot3D[Evaluate[r2 [te] /. w -> l-u-v/. v->0], {u, 0, 1},
Mesh-> 18, PlotStyle -> {Red, AbsoluteThickness [2] } ] ;
pb020 = {0, 0, 20} ; pb002 = {10, 0, 0} ; pb0ll = {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] ;
pbl01[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] ;
pbll0[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] ;
rl [ tb_] : = pb200 [ tb] u u + pb020 v v + pb002 w w + 2 pbl0l [ tb] uw + 2 pbll0 [ tb] u v +
2 pb0ll vw;
rfl [ tb_, s_, t_] := rl [tb] /.w->l-u-v/. {u -> t- s t, v-> s t} ;
flperem = ParametricPlot3D [Evaluate [rl [tb] /. w -> l - u - v / . u -> 0] , {v, 0, 1},
Mesh -> 18, PlotStyle -> {Red, AbsoluteThickness[2]}];
pj200 = {0, 10, 0} ; pj020 = {0, 0, 20} ; pjll0 = {0 , 10, 20};
pj002[tj_] := If [tj<= 0.5, {6, 8, 16} + 2 tj j ({0, 20, 40} - {6, 8, 16}) /. tj j -> tj,
{0, 20, 40} + 2 (tj j - 0.5) ({-6, 8, 16} - {0, 20, 40}) /. tj j -> tj] ;
pjl01[tj_] := If {tj <= 0,5, {5, 9, 3} + 2 tjj ({0, 15, 20}-{5, 9, 8}) /. tj j -> tj,
{0, 15, 20} + 2 (tjj - 0.5) ({-5, 9, 8} - {0, 15, 20}) / . tj j -> tj] ;
pj0ll [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}) /. tj j -> tj] ;
r3[ tj_] : = pj200uu+pj020vv + pj002[fcj] v w + 2 pjl0l[tj] u w + 2 pjll0 u v +
2 pj0ll [tj] vw;
rf3[tj_, s_, t_] := r3[tj] /.w->l-u-v/. {u -> t- s t, v-> s t} ;
f3perem = ParametricPlot3D [Evaluate [r3 [tj ] /. w -> 0 / . v->l-u], {u, 0, 1},
Mesh ->18, PlotStyle -> {Red, AbsoluteThickness[2]}];
rfel[tb_, te_, tj_, u_, v_] : =
и (rf3[tj, v, 1] + rf2[te, 0, u + v] - rf2[te, 0, 1] ) +
v (rfl [tb, 1, u+ v] + rf3[tj, 1 - u, 1] - rfl[ tb, 1, 1] ) +
(1-u- v) (rf2[te, 0, u] + rfl[ tb, 1, v] - rf 1 [ tb, 1, 0] ) ; (*три кривые rfl[tb,l,t] ,
rf2[te,0,l-t] и rЈ3[tj,l-t,l] , t Ha [0,1]*)
Manipulate[
Show[{ParametricPlot3D[Evaluate[rfl[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 . 9] , 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, t j , s, t] ] , {s, 0, 1}, {t, 0, 1 - s} ,
Mesh -> {13, 18}, PlotRange-> All] ,
f1perem, f2perem, fЗрегеm} , Viewpoint-> {-0 . 7 , 0.4, 1-3}, ImageSize -> {300 , 200},
Boxed -> False, Axes -> None, PlotRange -> {{-6. , 20.}, {-6., 20.}, {-10., 40.}},
Spherical Region -> True, ViewAngle-> 0 . 7 ],
{{tb, 0.4, "Зеленый ""треугольник '"}, 0., 1., Appearance->'Labeled" } ,
{ {te, 0.4, " Желтый ' ' треугольник '" } , 0 . , 1. , Appearance -> " Labeled" } ,
{ {t j , 0.4, " Синий " " треугольник' ' " } , 0 . , 1. , Appearance -> " Labeled" } ,
SaveDefinitions -> True] ]
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.