Опубликован: 27.12.2010 | Уровень: специалист | Доступ: платный
Лекция 6:

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

Алгоритм вычисления радиус-вектора 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 -кривой

r(t)=\frac{\sum_{i=1}^n N_{i,m}(t) \omega_i p_i}{\sum_{i=1}^n N_{i,m}(t) \omega_i}.
  1. Фиксируем 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.

  2. Для определенного выше индекса i_0, если 1 \le i_0 \le n + m, вычисляем единственное ненулевое значение ненормированного B -сплайна первого уровня (m=1):

    M_{i_0,1}(t)=\frac{1}{t_{i_0+1}-t_{i_0}}\\
(1 \le i_0 \le n+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} не вычисляется.

  3. С помощью соотношений Кокса - де Бура вычисляем все отличные от нуля в точке 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})}, ( 6.6)

    где 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}
  4. Вычисляем нормированные сплайны 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 ).

  5. Окончательно вычисляем 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}

Алгоритм де Бура вычисления радиус-вектора B-кривой

Теорема 6.4. Радиус-вектор r(t) B -кривой:

r(t)=\frac{\sum_{i=1}^nN_{i, m}(t) \omega_i p_i}{\sum_{i=1}^n N_{i, m}(t) \omega_i}=\frac{\omega r}{\omega}

может быть вычислен с помощью следующего алгоритма:

  1. Для заданного 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
  2. Вычисляем r(t) по формуле

    r(t)=\frac{r_{i_0}^{(m-1)}}{\omega_{i_0}^{(m-1)}}

Данный алгоритм иллюстрируется следующей диаграммой ( 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, образующих прямоугольную матрицу (сетку Безье):

P_{mn}=\begin{pmatrix}
p_{00} &\dots &p_{0n}\\
p_{10} &\dots &p_{1n}\\
\vdots &\ddots &\vdots \\
p_{m0} &\dots &p_{mn}
\end{pmatrix}

Поверхностью Безье порядка m * n, соответствующей сетке P_{mn}, называется поверхность

S_{00}(u, v; m, n)=(1-v+Fv)^n(1-u+Fu)^mp_{00}=\\
=(1-u+Fu)^m(1-v+Fv)^np_{00},\\
0 \le umv \le 1 ( 6.7)

где E - оператор сдвига вперед по первому индексу, F - оператор сдвига вперед по второму индексу:

Ep_{ij}=p_{(i+1)j},\\
Fp_{ij}=p_{i(j+1)}.

Операторы 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)

Формула (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])]

Геометрический смысл поверхности Безье

Поверхность Безье можно получить следующим образом:

  1. Строим (n + 1) кривую Безье по столбцам матрицы P_{mn}:

    r_j(u) = (1 - u + Eu)^mp_{0_j},        \\
j = 0, \dots ,n.
  2. Далее, начиная с каждой точки кривой Безье r_0(u) , строим кривые Безье, имеющие опорные точки на кривых r_j(u) (0 \le j \le n) , соответствующие одному и тому же значению параметра 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_{00}(u,v;0,0)=p_{00},\\
S_{00}(u, v;1,1) = (1-u + Eu)(1 -v + Fv)p_{00} =\\
= (1 - u)(1 - v)p_{00} + (1 - u)vp_{01} + u(1 - v)p_{10} + uvp_{11} ( 6.9)
S_{00}(u, v;2,2) = (1-u + Eu)^2(1 -v + Fv)^2p_{00} =\\
= (1 - u + Eu)(1 -v + Fv)S_{00}(u, v; 1,1) = \\
= (1 - u)(1 - v)S_{00}(u, v; 1,1) + (1 - u)vS_{01}(u, v; 1,1)+\\
+ u(1 - v)S_{10}(u, v; 1,1) + uvS_{11}(u, v; 1,1). ( 6.10)

Здесь через S_{ij} обозначены четырехугольные листы Безье, построенные из точек p ij, взятых в качестве угловых точек, аналогично формуле (6.9).

Аналогично формуле (6.10), любая точка на поверхности Безье степени m*n может быть представлена в виде точки четырехугольного листа Безье с углами в точках соответственно подобранных поверхностей Безье порядка (m-1)*(n-1):

S_{00}(u, v; m, n) = (1 - u)(1 - v)S_{00}(u, v;m-1,n- 1)+\\
+ (1 - u)vS_{01}(u, v;m-1,n-1)+u(1- v)S_{10}(u, v;m-1,n- 1)+\\
+ uvS_{11}(u,v;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) может быть получена с помощью следующего алгоритма:

  1. соединяем p_{00} и p_{01} рациональной кривой Безье S_{00}(0, v) с весами \omega_{00} и \omega_{01} ;
  2. соединяем p_{10} и p_{11} рациональной кривой Безье S_{00}(1, v) с весами \omega_{10}
  3. фиксируем v = v_0, соединяем точки S_{00}(0,v_0) и S_{00}(1,v_0) с весами (1 - v_0 + Fv_0) \omega_{00} и (1 - v_0 + Fv_0) \omega_{10} рациональной кривой Безье и берем на ней точку с параметром u = u_0 \in [0,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}.
Светлана Петрова
Светлана Петрова
Украина
Марина Семенова
Марина Семенова
Россия, г. Чебоксары