Появление данной статьи обусловленно тем, что по специфике работы одно время занимался программированием компьютерной графикой с использованием библиотеки OpenGL.
В рамках этой работы возникла необходимость построения достаточно сложной поверхности. В учебнике [1] (Урок 28. Фрагменты поверхностей Безье) приведен пример поверхности Безье основанной на полигональной сетке размерностью (4*4). У меня возник закономерный вопрос, каким образом построить поверхность Безье базирующейся на произвольной полигональной сетке, размерностью n m?
Изучая данный вопрос в Интернете, я нашел книгу [2], в разделе «6-11 Поверхность Безье» которой, приводится алгоритм построения поверхности Безье базирующейся на произвольной полигональной сетке размерностью n m.
В моей статье речь пойдет о реализации поверхности Безье базирующейся на произвольной полигональной сетке размерностью n m с примером реализации на языке C++. При этом многие теоретические определения взяты из 1 и 2.
Кривая Безье задается многоугольником. Так как базис Безье называется бернштейновским, сразу же известны некоторые свойства кривых Безье. Например:
Функции базиса вещественны;
Степень многочлена, определяющего порядок кривой на единицу меньше количества точек соответствующего многочлена;
Основа формы кривой повторяет очертания многоугольника;
Первая и последняя точка кривой совпадает с соответствующими точками многоугольника;
Векторы касательных в концах кривой совпадают с первой и последней сторонами многоугольника;
Кривая лежит внутри выпуклой оболочки многоугольника, т. е. внутри самого большого многоугольника построенного по заданным точкам;
Кривая обладает свойством уменьшения вариации. Это означает, что кривая пересекает любую прямую не чаще, чем определяющий многоугольник;
Кривая инвариантна относительно аффинных преобразований.
Кривые Безье базируются на простой функции, из которой получаются более сложные версии. Эта функция выглядит следующим образом:
t + (1 - t) = 1 (1)
Поскольку эта простая функция истинна для всех значений t, мы можем взять, и возвести в любую степень это выражение с обеих сторон, например, в кубическую:
(2)
Когда мы хотим вычислить значение точки находящейся на кривой, мы умножаем каждую компоненту уравнения на свою контрольную точку (как вектор) и находим сумму.
,
где
.
Тогда при
и при
При обобщении формул (1), (2), (3) получим математическое параметрическое представление кривой Безье в общем виде для произвольной степени, которая имеет следующий вид:
где базис Безье или Бернштейна или функция аппроксимации
где
которая определяет общее число сочетаний из n элементов в группах по i.
Формула (5) является весовой функцией Безье/Бернштейна.
Формула (6) была мной запрограммирована следующим образом:
float fact(int n)
{
float nn=1;
if (n==0||n==1) return 1.0;
for(int i=1;i<=n;i++) nn*=(float)i;
return nn;
}
float cs(int n,int i)
{
float s,d;
int k=n-i;
s= fact(n);
d= fact(i);
return (s/(d*fact(k)));
}
Два замечания к функциям:
Для функции float fact(int n) использовался алгоритм прямого вычисления факториала, который работает быстрее рекурсивного вызова;
Тип float использован для совместимости с функциями OpenGL принимающих в качестве параметров числа с плавающей точкой.
Формулу (4) можно представить в следующем виде
где b – называется базисной матрицей Безье. Формирование данной матрицы будет рассмотрено на примере поверхностей Безье
Теперь перейдем к формированию поверхности Безье.
Поверхность Безье базируется на наборе кривых Безье объединенных в полигональную сетку.
Тогда декартово или тензорное произведение Безье задается в виде
где
и
базисные функции Бернштейна в
параметрических направлениях u
и w
При этом
и
- комбинаторные сочетания вычисляются
аналогично формулы (6).
Как и для кривых Безье, для смешивающих функций используется базис Бернштейна, при этом многие свойства поверхности известны.
В матричном виде декартово произведение поверхности Безье задается выражением
Матрицы [N] и [M] являются базисными матрицами Безье, которые задаются уравнением типа
Запрограммируем функцию, по которой вычисляют матрицы [N] и [M]. Данные матрицы вычисляются один раз.
Знак для каждого элемента матрицы определяется в следующей функции
int orts(int n)
{
if (n==0 || (n%2)==0) return 1;
return -1;
}
Функция для вычисления матриц выглядит следующим образом
void matr(float *x,int n)
{
float x1,x2,x3,x0;
int i,j;
for(i=0;i<n+1;i++)
{
for(j=0;j<n+1;j++)
{
x0=0.0;
x1 = orts(n-j-i);
x2 = cs(n,j);
x3 = cs(n-j,n-j-i);
if((i+j)<=n) x0=x1*x2*x3;
*(x+i*(n+1)+j)=x0;
}
}
}
Вычисление матрицы осуществляется с учетом размерности в конкретном параметрическом направлении.
Например, для кубической кривой Безье данная матрица примет вид
Матрица
является вершинами задающей полигональной
сетки.
Вычисление текущего состояния
осуществляется с выбранным шагом в
параметрических направлениях u
и w
Формирование векторов
на каждой итерации осуществляется при
помощи функции
void coeff(float x0,float *x, int n)
{
*(x+n-1)=1.0;
for (int i=n-2;i>=0;i--)
*(x+i)=*(x+i+1)*x0;
}
где
float x0 – текущее состояние параметрического направления u или w ;
float *x – вектор характеризует сплайн Безье степени n-1;
int n – задает размерность вектора в параметрическом направлении.
Формула (9) вычисляется справа налево. Для этого используем функции приведенные ниже.
Умножение матрицы на вектор
void matrvekt(float *mtr,float *vkt0, float *vkt1,int n,int m)
{
int i,j;
for(i=0;i<n;i++)
{
*(vkt1+i)=0.0;
for(j=0;j<m;j++)
*(vkt1+i)+=*(mtr+i*m+j)**(vkt0+j);
}
}
где
*mtr – указатель на матрицу
*vkt0 – указатель на входной вектор размерностью m
*vkt1 – указатель на входной вектор размерностью n
n и m – размерность матрицы
Функция ymnvekt(…)
вычисляет скалярное произведение
двух векторов, результатом которого
будет текущее состояние
float ymnvekt(float *vkt0, float *vkt1,int n)
{
float dd=0.0;
for(int i=0;i<n;i++)
dd+=*(vkt0+i)**(vkt1+i);
return dd;
}
Результатом выше изложенного и с
использованием функций приведенных
раннее, мною была запрограммирована
функция GLuint
drt(…),
по которой осуществляется построение
поверхности Безье, базирующаяся на
произвольной полигональной сетке
размерностью
:
GLuint drt(float *bx,float *by,float *bz, GLuint textur, int n,int m,int k)
{
GLuint drawlist = glGenLists(1);
float *b2x = new float[(k+1)*(k+1)],
*b2y = new float[(k+1)*(k+1)],
*b2z = new float[(k+1)*(k+1)];
float *x1 = new float[m*m],
*x2 = new float[n*n],
*b1x = new float[n*m],
*b1y = new float[n*m],
*b1z = new float[n*m],
*m1 = new float[m],
*vkt0 = new float[n],
*m2 = new float[n],
*prom0 = new float[m],
*prom1 = new float[m],
*u = new float[k+1],
*v = new float[(k+1)*(k+1)];
float a,b,c;
int i,j;
matr(x2,(n-1));
matr(x1,(m-1));
for(i=0;i<m;i++)
for(j=0;j<n;j++)
{
*(b1x+i*n+j) = *(bx+i*n+j);
*(b1y+i*n+j) = *(by+i*n+j);
*(b1z+i*n+j) = *(bz+i*n+j);
}
for(i=0;i<=k;i++)
{
*(u+i)=i/(1.0*k);
coeff(*(u+i),m1,m);
for (j=0;j<=k;j++)
{
*(v+i*(k+1)+j)=j/(1.0*k);
coeff(*(v+i*(k+1)+j),m2,n);
matrvekt(x2, m2, vkt0, n, n);
matrvekt(b1x, vkt0, prom0, m, n);
matrvekt(x1, prom0, prom1, m, m);
*(b2x+i*(k+1)+j)=ymnvekt(m1, prom1, m);
matrvekt(b1y, vkt0, prom0, m, n);
matrvekt(x1, prom0, prom1, m, m);
*(b2y+i*(k+1)+j)=ymnvekt(m1, prom1, m);
matrvekt(b1z, vkt0, prom0, m, n);
matrvekt(x1, prom0, prom1, m, m);
*(b2z+i*(k+1)+j)=ymnvekt(m1, prom1, m);
}
}
glNewList(drawlist,GL_COMPILE);
glBindTexture(GL_TEXTURE_2D,textur);
glEnable(GL_TEXTURE_2D);
for(i=0;i<k;i++)
for(j=0;j<k;j++)
{
glBegin(GL_QUADS);
glTexCoord2f(*(v+i*(k+1)+j), *(u+i));
glVertex3f(*(b2x+i*(k+1)+j), *(b2y+i*(k+1)+j), *(b2z+i*(k+1)+j));
glTexCoord2f(*(v+i*(k+1)+j+1),*(u+i));
glVertex3f(*(b2x+i*(k+1)+j+1), *(b2y+i*(k+1)+j+1), *(b2z+i*(k+1)+j+1));
glTexCoord2f(*(v+i*(k+1)+j+1), *(u+i+1));
glVertex3f(*(b2x+(i+1)*(k+1)+j+1), *(b2y+(i+1)*(k+1)+j+1), *(b2z+(i+1)*(k+1)+j+1));
glTexCoord2f(*(v+i*(k+1)+j), *(u+i+1));
glVertex3f( *(b2x+(i+1)*(k+1)+j), *(b2y+(i+1)*(k+1)+j), *(b2z+(i+1)*(k+1)+j));
glEnd();
}
glDisable(GL_TEXTURE_2D);
glEndList();
delete[] b2x;
delete[] b2y;
delete[] b2z;
delete[] u;
delete[] v;
delete[] x1;
delete[] x2;
delete[] m1;
delete[] m2;
delete[] vkt0;
delete[] prom0;
delete[] prom1;
delete[] b1x;
delete[] b1y;
delete[] b1z;
return drawlist;
}
где
*bx, *by, *bz – матрицы по координатам x, y, z
GLuint textur – текстура накладываемая на поверхность Безье
n и m – размерность матриц по координатам x, y, z
k
– количество шагов разбиения на участке
0..1
в параметрических направлениях
.
Данные процедуры проверены мною на полигональной сетке размерностью n=13 и m=11.
При построении поверхности Безье
рекомендуется ограничивать количество
шагов разбиения в параметрических
направлениях
4000.
Литература:
http://www.opengl.org.ru/lesson/nehe28.htm - Народный учебник по OpenGL Урок 28. Фрагменты поверхностей Безье
Д. Роджерс, Дж. Адамс Математические основы компьютерной графики – М. Мир 2001
Роджерс Д. - Алгоритмические основы машинной графики (1989)
Роджерс Д. Математические основы машинной графики.1980