Математика и 3D Графика

Введение 3D

Начнем с координатной системы . В декартовой системе координат мы используем для точки 2 координаты : (x, y) . В 3-мерной мы добавляем 3-ю координату z . Мы будем использовать следующие направления осей :

ось X - направлена вправо .
ось Y - направлена вверх .
ось Z - положительная полуось направлена вглубь экрана .

Вектор

Что такое вектор ? Рассмотрим 2D-вектор P (4,5). Это отрезок , соединяющий начало координат с данной точкой , имеющий направление и длину . Эту длину условно обозначим как | P |. Длина вектора вычисляется по формуле :
| P | = sqrt( x² + y² )>

Рассмотрим 3D вектор: P(4, -5, 9). Его длина :

| P | = sqrt( x² + y² + z² )>

Матрицы

Традиционное начало : матрица - это 2-мерный массив , или таблица чисел . В основном мы будем рассматривать матрицы размером 4 * 4. Почему ? Потому что пространство трехмерно , а нам еще потребуется дополнительная строка и столбец . В 2D нам нужна матрица 3x3 . Вот пример :

4x4 идентичная матрица
1	0	0	0
0	1	0	0
0	0	1	0
0	0	0	1

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


10	-7	22	45
sin(a)	cos(a)	34	32
-35	28	17	6
45	-99	32	16

Операции над векторами и матрицами .

Начнем со сложения 2-х векторов - результат есть вектор :
( x₁ , y₁ , z₁ ) + ( x₂ , y₂ , z₂ ) = ( x₁ + x₂ , y₁ + y₂ , z₁ + z₂ )>
Умножение вектора на скаляр :
k · ( x, y, z ) = ( kx, ky, kz )>
Умножение 2-х векторов - по научному dot product, есть не вектор , а число : ( x₁ , y₁ , z₁ ) · ( x₂ , y₂ , z₂ ) = x₁x₂ + y₁y₂ + z₁z₂ >Для нахождения косинуса угла между 2-мя векторами нужно этот самый dot product разделить на произведение модулей этих векторов :
cos (V ^ W) =
V · W

| V | | W |
>Это может быть использовано при нахождении угла между плоскостью и источником света .

Далее , еще одно произведение 2-х векторов - т.н. cross product.

( x₁ , y₁ , z₁ ) X ( x₂ , y₂ , z₂ ) = ( y₁z₂ - z₁y₂ , z₁x₂ - x₁z₂ , x₁y₂ - y₁x₂ ) >
Что может быть использовано при вычислении нормали к плоскости .

Перейдем к сложению 2-х матриц . Правило простое : Пусть у нас i - номер строки , j - номер столбца . Нужно в первой матрице взять элемент с индексом (i, j) и сложить его с соответствующим элементом (i, j) второй матрицы , и все дела . Теперь о произведении 2-х матриц : Пусть у нас есть 2 матрицы - M и N . Так вот : M x N * НЕ *> есть то же самое , что и N x M.

Вот формула для умножения матриц 4x4
If A=(a_ij)_4x4 and B=(b_ij)_4x4, then
A x B=

₄ S> a_1jb_j1^j=1	₄ S> a_1jb_j2^j=1	₄ S> a_1jb_j3^j=1	₄ S> a_1jb_j4^j=1
₄ S> a_2jb_j1^j=1	₄ S> a_2jb_j2^j=1	₄ S> a_2jb_j3^j=1	₄ S> a_2jb_j4^j=1
₄ S> a_3jb_j1^j=1	₄ S> a_3jb_j2^j=1	₄ S> a_3jb_j3^j=1	₄ S> a_3jb_j4^j=1
₄ S> a_4jb_j1^j=1	₄ S> a_4jb_j2^j=1	₄ S> a_4jb_j3^j=1	₄ S> a_4jb_j4^j=1

if AxB=(c_ik)_4x4 , то можно написать :
c_ik = S>_{4, j=1}a_ijb_jk>Как говорится , сложно описать музыку словами : т.е. мы берем , допустим , первую строку 1-й матрицы и умножаем почленно на первый столбец 2-й матрицы , сумма чего будет элемент с индексом (1,1) в результирующей матрице .
Попробуем теперь умножить 3D-вектор на матрицу 4*4 - и получим 3D-вектор : if B=(b_ij)_4x4: ( a₁, a₂, a₃ ) x B = (S>a_ib_i1 + b_4,1, S>a_ib_i2 + b_4,2, S>a_ib_i3 + b_4,3 )>

Нужно последовательно перемножить 3 координаты вектора на каждый из столбиков матрицы с добавлением элемента , стоящего в 4-й строке .

Transformations

Мы уже встречались с формулами типа :
t_{( tx, ty )}: ( x, y ) ==> ( x + tx, y + ty )>
Это есть т.н. операция "translation" , или линейное изменение координат в 2D . Далее , масштабирование :
s_{( k )}: ( x, y ) ==> ( kx, ky )>
Вращение :
r( q> ): ( x, y ) ==> ( x cos(q>) - y sin(q>), x sin(q>) + y cos(q>) )>
Это все для 2D, а для 3D нужно добавить координату z .При масштабировании умножаем z на k.

Пусть мы имеем вектор ( x, y, z ) и несколько матриц trasformation, по одной для каждого типа . Мы будем умножать матрицу на вектор , и получим результирующий вектор . Вот эти 3D transformation-матрицы :

Матрица для простого 3D-translation (tx, ty, tz)
1	0	0	0
0	1	0	0
0	0	1	0
tx	ty	tz	1

Матрица для 3D-масштабирования (sx, sy, sz)
sz	0	0	0
0	sy	0	0
0	0	sx	0
0	0	0	1

Матрица для 3D-вращения вокруг оси x q>
0	0	0
cos(q>)	sin(q>)	0
-sin(q>)	cos(q>)	0
0	0	1

Матрица для 3D-вращения вокруг оси y q>
cos(q>)	0	-sin(q>)	0
0	1	0	0
sin(q>)	0	cos(q>)	0
0	0	0	1

Матрица для 3D-вращения вокруг оси z q>
cos(q>)	sin(q>)	0	0
-sin(q>)	cos(q>)	0	0
0	0	1	0
0	0	0	1

Плоскости & и всякие к ним нормали ...

Плоскость - это ориентированная поверхность , уравнение которой можно записать так :
Ax + By + Cz + D = 0>
где A, B, C мы обзовем нормалями к плоскости, и D - расстояние от плоскости до начала координат. Нормаль - это вектор , перпендикулярный к плоскости. Вычислим эту нормаль с помощью "cross products" двух ребер на плоскости . Для идентификации 2-х ребер нам понадобятся 3 точки . Если P₁ есть первая точка, P₂ вторая, и P₃ третья, то :
Edge₁ = P₁ - P₂>
и
Edge₂ = P₃ - P₂>
и вычисляем нормаль : Normal = Edge₁ X Edge₂>
Запишем уравнение плоскости в виде : D = - (Ax + By + Cz)>
или
D = - (A·P₁.x + B·P₁.y + C·P₁.z)>
или :
D = - Normal · P₁>
Для вычисления A, B, C : A = y₁ ( z₂ - z₃ ) + y₂ ( z₃ - z₁ ) + y₃ ( z₁ - z₂ ) B = z₁ ( x₂ - x₃ ) + z₂ ( x₃ - x₁ ) + z₃ ( x₁ - x₂ ) C= x₁ ( y₂ - y₃ ) + x₂ ( y₃ - y₁ ) + x₃ ( y₁ - y₂ ) D = - x₁ ( y₂z₃ - y₃z₂ ) - x₂ ( y₃z₁ - y₁z₃ ) - x₃ ( y₁z₂ - y₂z₁ )>

3D Transformations

Координаты

Приступим к коду . Начнем с основных структур . Какие типы координат нам нужны ?

Экранные координаты: 2D координаты дисплея
Локальные координаты: 3D-координаты , привязанные к самому обьекту

Мировые координаты: 3D-координаты всего мирового пространства
Aligned Coordinates: Координата наблюдателя или камеры

Вот базовые структуры :

typedef struct
{
   short x, y;
}_2D;

typedef struct
{
   float x, y, z;
}_3D;

Теперь определим структуру для вершины или "vertex" . Вершина - это общая точка для 2-х прилегающих ребер (edge) , она может быть описана как система из 3-х 3D-векторов :

typedef struct
{
   _3D Local;
   _3D World;
   _3D Aligned;
}Vertex_t;

Матрицы размера 4x4 мы будем хранить в формате floating-point-чисел.

float matrix[4][4];

Несколько функций для работы с матрицами - копирование :

void MAT_Copy(float source[4][4], float dest[4][4])
{
    int i,j;
    for(i=0; i<4; i++)
      for(j=0; j<4; j++)
         dest[i][j]=source[i][j];
}

Умножение 2-х матриц :

void MAT_Mult(float mat1[4][4], float mat2[4][4], float dest[4][4])
{
   int i,j;
   for(i=0; i<4; i++)
      for(j=0; j<4; j++)
         dest[i][j]=mat1[i][0]*mat2[0][j]+
                    mat1[i][1]*mat2[1][j]+
                    mat1[i][2]*mat2[2][j]+
                    mat1[i][3]*mat2[3][j];
}

Умножение вектора на матрицу :

void VEC_MultMatrix(_3D *Source,float mat[4][4],_3D *Dest)
{
    Dest->x=Source->x*mat[0][0]+
            Source->y*mat[1][0]+
            Source->z*mat[2][0]+
                      mat[3][0];
    Dest->y=Source->x*mat[0][1]+
            Source->y*mat[1][1]+
            Source->z*mat[2][1]+
                      mat[3][1];
    Dest->z=Source->x*mat[0][2]+
            Source->y*mat[1][2]+
            Source->z*mat[2][2]+
                      mat[3][2];
}

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

float SinTable[256], CosTable[256];

Далее можно использовать вот такие макросы :

#define SIN(x) SinTable[ABS((int)x&255)]
#define COS(x) CosTable[ABS((int)x&255)]

Инициализация таблиц :

void M3D_Init()
{
   int d;
   for(d=0; d<256; d++)
   {
       SinTable[d]=sin(d*PI/128.0);
       CosTable[d]=cos(d*PI/128.0);
   }
}

Создание transformation-матриц

Эти матрицы могут быть записаны так :

float mat1[4][4], mat2[4][4];

void MAT_Identity(float mat[4][4])
{
    mat[0][0]=1; mat[0][1]=0; mat[0][2]=0; mat[0][3]=0;
    mat[1][0]=0; mat[1][1]=1; mat[1][2]=0; mat[1][3]=0;
    mat[2][0]=0; mat[2][1]=0; mat[2][2]=1; mat[2][3]=0;
    mat[3][0]=0; mat[3][1]=0; mat[3][2]=0; mat[3][3]=1;
}

void TR_Translate(float matrix[4][4],float tx,float ty,float tz)
{
   float tmat[4][4];
   tmat[0][0]=1;  tmat[0][1]=0;  tmat[0][2]=0;  tmat[0][3]=0;
   tmat[1][0]=0;  tmat[1][1]=1;  tmat[1][2]=0;  tmat[1][3]=0;
   tmat[2][0]=0;  tmat[2][1]=0;  tmat[2][2]=1;  tmat[2][3]=0;
   tmat[3][0]=tx; tmat[3][1]=ty; tmat[3][2]=tz; tmat[3][3]=1;
   MAT_Mult(matrix,tmat,mat1);
   MAT_Copy(mat1,matrix);
}

void TR_Scale(float matrix[4][4],float sx,float sy, float sz)
{
   float smat[4][4];
   smat[0][0]=sx; smat[0][1]=0;  smat[0][2]=0;  smat[0][3]=0;
   smat[1][0]=0;  smat[1][1]=sy; smat[1][2]=0;  smat[1][3]=0;
   smat[2][0]=0;  smat[2][1]=0;  smat[2][2]=sz; smat[2][3]=0;
   smat[3][0]=0;  smat[3][1]=0;  smat[3][2]=0;  smat[3][3]=1;
   MAT_Mult(matrix,smat,mat1);
   MAT_Copy(mat1,matrix);
}

void TR_Rotate(float matrix[4][4],int ax,int ay,int az)
{
   float xmat[4][4], ymat[4][4], zmat[4][4];
   xmat[0][0]=1;        xmat[0][1]=0;        xmat[0][2]=0;        xmat[0][3]=0;
   xmat[1][0]=0;        xmat[1][1]=COS(ax);  xmat[1][2]=SIN(ax);  xmat[1][3]=0;
   xmat[2][0]=0;        xmat[2][1]=-SIN(ax); xmat[2][2]=COS(ax);  xmat[2][3]=0;
   xmat[3][0]=0;        xmat[3][1]=0;        xmat[3][2]=0;        xmat[3][3]=1;

   ymat[0][0]=COS(ay);  ymat[0][1]=0;        ymat[0][2]=-SIN(ay); ymat[0][3]=0;
   ymat[1][0]=0;        ymat[1][1]=1;        ymat[1][2]=0;        ymat[1][3]=0;
   ymat[2][0]=SIN(ay);  ymat[2][1]=0;        ymat[2][2]=COS(ay);  ymat[2][3]=0;
   ymat[3][0]=0;        ymat[3][1]=0;        ymat[3][2]=0;        ymat[3][3]=1;

   zmat[0][0]=COS(az);  zmat[0][1]=SIN(az);  zmat[0][2]=0;        zmat[0][3]=0;
   zmat[1][0]=-SIN(az); zmat[1][1]=COS(az);  zmat[1][2]=0;        zmat[1][3]=0;
   zmat[2][0]=0;        zmat[2][1]=0;        zmat[2][2]=1;        zmat[2][3]=0;
   zmat[3][0]=0;        zmat[3][1]=0;        zmat[3][2]=0;        zmat[3][3]=1;

   MAT_Mult(matrix,ymat,mat1);
   MAT_Mult(mat1,xmat,mat2);
   MAT_Mult(mat2,zmat,matrix);
}

Как создать перспективу

Интересный вопрос , как для программистов , так и для художников . Вот формула для преобразования 3D в 2D : P( _f ):(x, y, z)==>( fx / z + XOrigin, fy / z + YOrigin )>
Здесь f "фокусное расстояние" , т.е. расстояние от камеры до экрана - обычно в диапазоне от 80 до 200 . XOrigin и YOrigin - координаты центра дисплея . Код :

#define FOCAL_DISTANCE 200
void Project(vertex_t * Vertex)
{
   if(!Vertex->Aligned.z)
        Vertex->Aligned.z=1;
   Vertex->Screen.x = FOCAL_DISTANCE * Vertex->Aligned.x / Vertex->Aligned.z + XOrigin;
   Vertex->Screen.y = FOCAL_DISTANCE * Vertex->Aligned.y / Vertex->Aligned.z + YOrigin;
}

Transforming-обьекты

Итак , перед тем , как выполнять 3D-преобразования , нужно :

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

Полигоны

Базовая структура

Нам нужна 2D и 3D структуры :

typedef struct
{
   _2D Points[20];
    int PointsCount;
   int Texture;
}Polygon2D_t;

typedef struct
{
    int Count;
   int * Vertex;
   int Texture;

   Vertex_t P,M,N;
}Polygon_t;

Если представить , например , куб , мы увидим , что три грани имеют одну общую вершину . Зачем одну вершину хранить в 3 разных местах памяти ? Для этого мы создадим структуру "ОБЬЕКТ" и будем хранить там вершины с помощью индексов :

typedef struct
{
   int VertexCount;
   int PolygonCount;
   Vertex_t * Vertex;
   Polygon_t * Polygon;
   _3D Scaling;
   _3D Position;
   _3D Angle;
   int NeedUpdate;
}Object_t;

Треугольник

Прорисовка треугольника есть частный и простейший случай прорисовки полигона . Потому начнем с полигона , имеющего 3 вершины :

void POLY_Draw(Polygon2D_t *Polygon)
{
    _2D P1,P2,P3;
   int i;

   P1 = Polygon->Points[0];
   for(i=1; i < Polygon->PointsCount-1; i++)
   {
       P2=Polygon->Points[i];
       P3=Polygon->Points[i+1];
       POLY_Triangle(P1,P2,P3,Polygon->Texture);
   }
}

Сканирование треугольника

Для этого нам надо знать начальную и конечную координату х для каждой горизонтальной линии . Начнем со следующих макросов :

#define MIN(a,b) ((a<b)?(a):(b))
#define MAX(a,b) ((a>b)?(a):(b))
#define MaxPoint(a,b) ((a.y > b.y) ? a : b)
#define MinPoint(a,b) ((b.y > a.y) ? a : b)

#define MaxPoint3(a,b,c) MaxPoint(MaxPoint(a,b),MaxPoint(b,c))
#define MidPoint3(a,b,c) MaxPoint(MinPoint(a,b),MinPoint(a,c))
#define MinPoint3(a,b,c) MinPoint(MinPoint(a,b),MinPoint(b,c))

Теперь сама функция :

void POLY_Triangle(_2D   p1,_2D p2,_2D p3,char c)
{
   _2D p1d,p2d,p3d;
   int xd1,yd1,xd2,yd2,i;
   int Lx,Rx;

Произведем сортировку 3-х вершин :

   p1d = MinPoint3(p1,p2,p3);
   p2d = MidPoint3(p2,p3,p1);
   p3d = MaxPoint3(p3,p1,p2);

Обратите внимание на порядок параметров в функции MaxPoint3 - он существенен . Далее проверка :

   if(p2.y < p1.y)
   {
      p1d=MinPoint3(p2,p1,p3);
      p2d=MidPoint3(p1,p3,p2);
   }

Дальше геометрия :

   xd1=p2d.x-p1d.x;
   yd1=p2d.y-p1d.y;
   xd2=p3d.x-p1d.x;
   yd2=p3d.y-p1d.y;

   if(yd1)
      for(i=p1d.y; i<=p2d.y; i++)
      {

Вычисляем координаты начала и конца сканирующей горизонтальной линии :

          Lx = p1d.x + ((i - p1d.y) * xd1) / yd1;
          Rx = p1d.x + ((i - p1d.y) * xd2) / yd2;

И прорисовываем ее :

          if(Lx!=Rx)
             VID_HLine(MIN(Lx,Rx),MAX(Lx,Rx),i,c);
      }

Так , верхнюю половину треугольника мы прорисовали , займемся нижней :

   xd1=p3d.x-p2d.x;
   yd1=p3d.y-p2d.y;

   if(yd1)
      for(i   = p2d.y;   i <= p3d.y;   i++)
      {
         Lx =   p1d.x   + ((i   - p1d.y)   * xd2) / yd2;
         Rx =   p2d.x   + ((i   - p2d.y)   * xd1) / yd1;
         if(Lx!=Rx)
            VID_HLine(MIN(Lx,Rx),MAX(Lx,Rx),i,c);
      }
}

Осталась малость - собственно прорисовка :

void VID_HLine(int x1, int x2, int y, char c)
{
   int x;
   for(x=x1; x<=x2; x++)
      putpixel(x, y, c);
}

Клиппинг методом Sutherland-Hodgman

Обзор

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

typedef struct
{
    int Count;
   _3D Vertex[20];
}CPolygon_t;

Функция клиппинга :

void M3D_Project(CPolygon_t *Polygon,Polygon2D_t *Clipped,int focaldistance)
{
   int v;
    for(v=0; v<Polygon->Count; v++)
   {
      if(!Polygon->Vertex[v].z)Polygon->Vertex[v].z++;
         Clipped->Points[v].x=Polygon->Vertex[v].x*focaldistance/
            Polygon->Vertex[v].z+Origin.x;
         Clipped->Points[v].y=Polygon->Vertex[v].y*focaldistance/
            Polygon->Vertex[v].z+Origin.y;
   }
   Clipped->PointsCount=Polygon->Count;
}

Z-Clipping

Определим макросы :

WORD ZMin=20;
#define INIT_ZDELTAS dold=V2.z-V1.z;  dnew=ZMin-V1.z;
#define INIT_ZCLIP INIT_ZDELTAS if(dold)m=dnew/dold;

Затем определим функцию , которая будет иметь 3 параметра : указатель на полигон и 2 координаты для ребра :

void CLIP_Front(CPolygon_t *Polygon,_3D V1,_3D V2)
{
   float dold,dnew, m=1;
   INIT_ZCLIP

Далее мы определяем на пересечение данного ребра с граничными координатами :

Если внутри :

   if ( (V1.z>=ZMin) && (V2.z>=ZMin) )
      Polygon->Vertex[Polygon->Count++]=V2;

Если выходит за пределы :

   if ( (V1.z>=ZMin) && (V2.z<ZMin) )
   {
      Polygon->Vertex[Polygon->Count   ].x=V1.x + (V2.x-V1.x)*m;
      Polygon->Vertex[Polygon->Count   ].y=V1.y + (V2.y-V1.y)*m;
      Polygon->Vertex[Polygon->Count++ ].z=ZMin;
   }

Если совпадает :

   if ( (V1.z<ZMin) && (V2.z>=ZMin) )
   {
      Polygon->Vertex[Polygon->Count   ].x=V1.x + (V2.x-V1.x)*m;
      Polygon->Vertex[Polygon->Count   ].y=V1.y + (V2.y-V1.y)*m;
      Polygon->Vertex[Polygon->Count++ ].z=ZMin;
      Polygon->Vertex[Polygon->Count++ ]=V2;
   }

Теперь нам нужно написать Z-Clipping-процедуру . Запишем макрос для вершины полигона :

#define Vert(x) Object->Vertex[Polygon->Vertex[x]]

Сама функция :

void CLIP_Z(Polygon_t *Polygon,Object_t *Object,CPolygon_t *ZClipped)
{
   int d,v;
   ZClipped->Count=0;
   for (v=0; v<Polygon->Count; v++)
   {
      d=v+1;
      if(d==Polygon->Count)d=0;
      CLIP_Front(ZClipped, Vert(v).Aligned,Vert(d).Aligned);
   }
}

Выберем разрешение 320*200. Определим макросы :

#define INIT_DELTAS dx=V2.x-V1.x;  dy=V2.y-V1.y;
#define INIT_CLIP INIT_DELTAS if(dx)m=dy/dx;

Выберем левый верхний и правый нижний углы экрана :

_2D TopLeft, DownRight;

Определим функции :

_2D P2D(short xshort y)
{
   _2D Temp;
   Temp.x=x;
   Temp.y=y;
   return Temp;
}
_3D P3D(float x,float y,float z)
{
   _3D Temp;
   Temp.x=x;
   Temp.y=y;
   Temp.z=z;
   return Temp;
}

Проинициализируем :

TopLeft=P2D(0, 0);
DownRight=P2D(319, 199);

И далее - 4 функции для edge clipping :

void CLIP_Left(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
   float   dx,dy, m=1;
   INIT_CLIP

   // ************OK************
   if ( (V1.x>=TopLeft.x) && (V2.x>=TopLeft.x) )
   Polygon->Points[Polygon->PointsCount++]=V2;
   // *********LEAVING**********
   if ( (V1.x>=TopLeft.x) && (V2.x<TopLeft.x) )
   {
      Polygon->Points[Polygon->PointsCount].x=TopLeft.x;
      Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(TopLeft.x-V1.x);
   }
   // ********ENTERING*********
   if ( (V1.x<TopLeft.x) && (V2.x>=TopLeft.x) )
   {
      Polygon->Points[Polygon->PointsCount].x=TopLeft.x;
      Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(TopLeft.x-V1.x);
      Polygon->Points[Polygon->PointsCount++]=V2;
   }
}
void CLIP_Right(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
   float dx,dy, m=1;
   INIT_CLIP
   // ************OK************
   if ( (V1.x<=DownRight.x) && (V2.x<=DownRight.x) )
      Polygon->Points[Polygon->PointsCount++]=V2;
   // *********LEAVING**********
   if ( (V1.x<=DownRight.x) && (V2.x>DownRight.x) )
   {
      Polygon->Points[Polygon->PointsCount].x=DownRight.x;
      Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(DownRight.x-V1.x);
   }
   // ********ENTERING*********
   if ( (V1.x>DownRight.x) && (V2.x<=DownRight.x) )
   {
      Polygon->Points[Polygon->PointsCount].x=DownRight.x;
      Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(DownRight.x-V1.x);
      Polygon->Points[Polygon->PointsCount++]=V2;
   }
}
/*
=================
CLIP_Top
=================
*/
void CLIP_Top(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
   float   dx,dy, m=1;
   INIT_CLIP
   // ************OK************
   if ( (V1.y>=TopLeft.y) && (V2.y>=TopLeft.y) )
      Polygon->Points[Polygon->PointsCount++]=V2;
   // *********LEAVING**********
   if ( (V1.y>=TopLeft.y) && (V2.y<TopLeft.y) )
   {
      if(dx)
         Polygon->Points[Polygon->PointsCount].x=V1.x+(TopLeft.y-V1.y)/m;
      else
         Polygon->Points[Polygon->PointsCount].x=V1.x;
      Polygon->Points[Polygon->PointsCount++].y=TopLeft.y;
   }
   // ********ENTERING*********
   if ( (V1.y<TopLeft.y) && (V2.y>=TopLeft.y) )
   {
      if(dx)
         Polygon->Points[Polygon->PointsCount].x=V1.x+(TopLeft.y-V1.y)/m;
      else
         Polygon->Points[Polygon->PointsCount].x=V1.x;
      Polygon->Points[Polygon->PointsCount++].y=TopLeft.y;
      Polygon->Points[Polygon->PointsCount++]=V2;
   }
}
void CLIP_Bottom(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
   float dx,dy, m=1;
   INIT_CLIP
   // ************OK************
   if ( (V1.y<=DownRight.y) && (V2.y<=DownRight.y) )
      Polygon->Points[Polygon->PointsCount++]=V2;
   // *********LEAVING**********
   if ( (V1.y<=DownRight.y) && (V2.y>DownRight.y) )
   {
      if(dx)
         Polygon->Points[Polygon->PointsCount].x=V1.x+(DownRight.y-V1.y)/m;
      else
         Polygon->Points[Polygon->PointsCount].x=V1.x;
      Polygon->Points[Polygon->PointsCount++].y=DownRight.y;
   }
   // ********ENTERING*********
   if ( (V1.y>DownRight.y) && (V2.y<=DownRight.y) )
   {
      if(dx)
         Polygon->Points[Polygon->PointsCount].x=V1.x+(DownRight.y-V1.y)/m;
      else
         Polygon->Points[Polygon->PointsCount].x=V1.x;
      Polygon->Points[Polygon->PointsCount++].y=DownRight.y;
      Polygon->Points[Polygon->PointsCount++]=V2;
   }
}

Введем глобальные переменные :

polygon2D_t TmpPoly;

 Функция клиппинга :

void CLIP_Polygon(Polygon2D_t *Polygon,Polygon2D_t *Clipped)
{
   int v,d;

   Clipped->PointsCount=0;
   TmpPoly.PointsCount=0;

   for (v=0; v<Polygon->PointsCount; v++)
   {
      d=v+1;
      if(d==Polygon->PointsCount)d=0;
      CLIP_Left(&TmpPoly, Polygon->Points[v],Polygon->Points[d]);
   }
   for (v=0; v<TmpPoly.PointsCount; v++)
   {
      d=v+1;
      if(d==TmpPoly.PointsCount)d=0;
      CLIP_Right(Clipped, TmpPoly.Points[v],TmpPoly.Points[d]);
   }
   TmpPoly.PointsCount=0;
   for (v=0; v<Clipped->PointsCount; v++)
   {
      d=v+1;
      if(d==Clipped->PointsCount)d=0;
      CLIP_Top(&TmpPoly, Clipped->Points[v],Clipped->Points[d]);
   }
   Clipped->PointsCount=0;
   for (v=0; v<TmpPoly.PointsCount; v++)
   {
      d=v+1;
      if(d==TmpPoly.PointsCount)d=0;
      CLIP_Bottom(Clipped, TmpPoly.Points[v],TmpPoly.Points[d]);
   }
}

Удаление невидимых поверхностей

Удаление невидимых поверхностей - сокращенно HSR - есть сердцевина любого движка .

Алгоритм художника

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

BSP trees

Быстр , но труден в реализации

Z-Buffering

При числе полигонов более 5000 становится быстрее , чем алгоритм художника . Легок в реализации Требует много памяти

Впредь мы и будем говорить о Z-Buffering. Определим глобальные переменные :

float A,B,C,D;
BOOL backface;

Начнем с перевода координат в формат float :

void ENG3D_SetPlane(Polygon_t *Polygon,Object_t *Object)
{
   float x1=Vert(0).Aligned.x;
   float x2=Vert(1).Aligned.x;
   float x3=Vert(2).Aligned.x;
   float y1=Vert(0).Aligned.y;
   float y2=Vert(1).Aligned.y;
   float y3=Vert(2).Aligned.y;
   float z1=Vert(0).Aligned.z;
   float z2=Vert(1).Aligned.z;
   float z3=Vert(2).Aligned.z;

Вычислим константы для идентификации плоскости :

   A=y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2);
   B=z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2);
   C=x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2);
   D=-x1*(y2*z3-y3*z2)-x2*(y3*z1-y1*z3)-x3*(y1*z2-y2*z1);

Проверим , в какую сторону плоскость обращена видимой стороной :

   backface=D<0;
}

Z-Buffering

Z-Buffering заключается в организации 2-мерного массива и хранения в нем координаты z для каждой точки , выводимой на экран . Необходимо в начале продекларировать глобальный массив , выделив для него память :

typedef long ZBUFTYPE;
ZBUFTYPE *zbuffer;
zbuffer=(ZBUFTYPE *)malloc(sizeof(ZBUFTYPE)*MEMORYSIZE);

Затем проинициализируем этот z-буфер минимально возможным числом :

   int c;
   for(c=0; c<MEMORYSIZE; c++)
      zbuffer[c]=-32767;

Как определить z координату ? Координата z нужна для вершин , а не для каждой точки полигона .
u = f · x / z>
и
v = f · y / z>
где u есть координата x на экране, минус XOrigin, и v есть координата y на экране, минус YOrigin. Уравнение плоскости :
Ax + By + Cz + D = 0>
Поскольку :
x = uz / f>
и
y = vz / f>
Подставив это в уравнение плоскости , получим :
A(uz / f) + B(vz / f) + Cz = -D >
отсюда определим z :
z(A(u / f) + B(v / f) + C) = -D>
Или : z = -D / (A(u / f) + B(v / f) + C)>
Переворачиваем дробь : 1 / z = -(A(u / f) + B(v / f) +C) / D>
1 / z = -(A / (fD))u - (B / (fD))v - C / D>
Или :
1 / z = -(A / (fD))u₁ - (B / (fD))v - C / D>
Инкремент для каждого пиксела : -(A / (fD))>
Код :

   #define FIXMUL (1<<20)

   int offset=y*MODEINFO.XResolution+x1;
   int i=x1-Origin.x, j=y-Origin.y;
   float z_,dz_;
   ZBUFTYPE z,dz;

   //Initialize 1/z value (z: 1/z)
   dz_=((A/(float)Focal_Distance)/-D);
   z_=((dz_*i)+( (B*j/(float)Focal_Distance) + C) /-D);
   dz=dz_*FIXMUL;
   z=z_*FIXMUL;

И для каждого пиксела мы делаем :

   if(z>ZBuffer[offset])
   {
      zbuffer[offset]=z;
      SCREENBUFFER[offset]=color;
   }
   z+=dz;

3D Texture Mapping

Первую вещь , которую нужно сделать - это завести текстурные координаты . Инициализация текстурного массива :

#define MAXTEXTURES 16
bitmap_t Textures[MAXTEXTURES];

Пусть образ будет храниться в .pcx-файле , и его размер 64*64 . Будем использовать текстурные координаты в структуре polygon_t :

vertex_t P,M,N;

Мы проинициализируем их после создания полигонов . P - это начало координат текстуры , M - горизонтальный конец текстуры , N - вертикальный конец текстуры .

void TEX_Setup(Polygon_t * Polygon, Object_t *Object)
{
   Polygon->P.Local=P3D(Vert(1).Local.x,
                        Vert(1).Local.y,
                        Vert(1).Local.z);
   Polygon->M.Local=P3D(Vert(0).Local.x,
                        Vert(0).Local.y,
                        Vert(0).Local.z);
   Polygon->N.Local=P3D(Vert(2).Local.x,
                        Vert(2).Local.y,
                        Vert(2).Local.z);
}

Создадим для них world transforming и align transforming - функции :

void TR_Object(Object_t *Object, float matrix[4][4])
{
   int v,p;
   for(v=0; v<Object->VertexCount; v++)
      VEC_MultMatrix(&Object->Vertex[v].Local,matrix,&Object->Vertex[v].World);
   for(p=0; p<Object->PolygonCount; p++)
   {
      VEC_MultMatrix(&Object->Polygon[p].P.Local,matrix,&Object->Polygon[p].P.World);
      VEC_MultMatrix(&Object->Polygon[p].M.Local,matrix,&Object->Polygon[p].M.World);
      VEC_MultMatrix(&Object->Polygon[p].N.Local,matrix,&Object->Polygon[p].N.World);
   }
}
void TR_AlignObject(Object_t *Object, float matrix[4][4])
{
   int v,p;
   for(v=0; v<Object->VertexCount; v++)
      VEC_MultMatrix(&Object->Vertex[v].World,matrix,&Object->Vertex[v].Aligned);
   for(p=0; p<Object->PolygonCount; p++)
   {
      VEC_MultMatrix(&Object->Polygon[p].P.World,matrix,&Object->Polygon[p].P.Aligned);
      VEC_MultMatrix(&Object->Polygon[p].M.World,matrix,&Object->Polygon[p].M.Aligned);
      VEC_MultMatrix(&Object->Polygon[p].N.World,matrix,&Object->Polygon[p].N.Aligned);
   }
}

Волшебные числа

Определим горизонтальную и вертикальную координату для каждого пиксела в нашем битмап-образе , который мы собираемся выводить на экран : Эти текстурные координаты назовем u и v. Уравнения :
u = a * TEXTURE_SIZE / c>
и v = b * TEXTURE_SIZE / c>
Эти магические координаты эквивалентны следующему :
a = O_x + V_x j + H_x i>
b = O_y + V_y j + H_y i>
c = O_z + V_z j + H_z i>
Эти O, H, V числа и есть магические числа . Они вычисляются так : O_x = N_xP_y - N_yP_x
H_x = N_yP_z - N_zP_y
V_x = N_zP_x - N_xP_z
O_y = M_xP_y - M_yP_x
H_y = M_yP_z - M_zP_y
V_y = M_zP_x - M_xP_z
O_z = M_yN_x - M_xN_y
H_z = M_zN_y - M_yN_z
V_z = M_xN_z - M_zN_x
>

Perspective Correct Texture Mapping

Вычисление чисел O,H, V требует небольшой модификации :

   //Used to fix errors that happen when the numbers get too big
   #define FIX_FACTOR 1/640

   //Initialize texture vectors
   P=Polygon->P.Aligned;
   M=VEC_Difference(Polygon->M.Aligned,Polygon->P.Aligned);
   N=VEC_Difference(Polygon->N.Aligned,Polygon->P.Aligned);

   P.x*=Focal_Distance;
   P.y*=Focal_Distance;

   M.x*=Focal_Distance;
   M.y*=Focal_Distance;

   N.x*=Focal_Distance;
   N.y*=Focal_Distance;

Функция VEC_Difference :

_3D VEC_Difference(_3D Vector1, _3D Vector2)
{
    return P3D(Vector1.x-Vector2.x,Vector1.y-Vector2.y,Vector1.z-Vector2.z);
}

Определимся :

_3D O, H, V;

ENG3D_SetPlane:

   H.x=(N.y*P.z-N.z*P.y)*FIX_FACTOR;
   V.x=(N.z*P.x-N.x*P.z)*FIX_FACTOR;
   O.x=(N.x*P.y-N.y*P.x)*FIX_FACTOR;

   H.z=(M.z*N.y-M.y*N.z)*FIX_FACTOR;
   V.z=(M.x*N.z-M.z*N.x)*FIX_FACTOR;
   O.z=(M.y*N.x-M.x*N.y)*FIX_FACTOR;

   H.y=(M.y*P.z-M.z*P.y)*FIX_FACTOR;
   V.y=(M.z*P.x-M.x*P.z)*FIX_FACTOR;
   O.y=(M.x*P.y-M.y*P.x)*FIX_FACTOR;

Подкорректируем функцию сканирования VID_HLine для TEX_HLine для texture mapping :

   a=-((long)O.x+((long)V.x*(long)j)+((long)H.x*(long)i))*64L;
   b= ((long)O.y+((long)V.y*(long)j)+((long)H.y*(long)i))*64L;
   c= ((long)O.z+((long)V.z*(long)j)+((long)H.z*(long)i));
   long Hx,Hy,Hz;
   int u,v;
   BYTE color=0;
   BYTE *mapping=Textures[texture].Picture;

Умножим H.x и H.y на 64 , чтобы не делать этого для каждого пиксела :

   Hx=H.x*-64;
   Hy=H.y*64;
   Hz=H.z;

И для каждого пиксела выполняем :

         if(c)
         {
            u=a/c;
            v=b/c;
            color=mapping[((v&63)<<6)+(u&63)];
            if(color)
            {
               zbuffer[offset]=z;
               SCREENBUFFER[offset]=LightTable[light][color];
            }
         }
         a+=Hx;
         b+=Hy;
         c+=Hz;

3D Shading

Вычисление нормалей

Несколько функций :

float VEC_DotProduct(_3D Vector1, _3D Vector2)
{
    return (Vector1.x*Vector2.x+Vector1.y*Vector2.y+Vector1.z*Vector2.z);
}
_3D VEC_CrossProduct(_3D Vector1, _3D Vector2)
{
   return P3D(Vector1.y*Vector2.z-Vector1.z*Vector2.y,
              Vector1.z*Vector2.x-Vector1.x*Vector2.z,
              Vector1.x*Vector2.y-Vector1.y*Vector2.x);
}
void VEC_Normalize(_3D * Vector)
{
   float m=sqrt(Vector->x*Vector->x+Vector->y*Vector->y+Vector->z*Vector->z);
   Vector->x/=m;
   Vector->y/=m;
   Vector->z/=m;
}

Для 3D shading, нужно добавить ::

   //Compute the normal of the plane
   PlaneNormal=VEC_CrossProduct(P3D(x2-x1,y2-y1,z2-z1),P3D(x3-x1,y3-y1,z3-z1));
   VEC_Normalize(&PlaneNormal);

Как известно , косинус угла между 2-мя векторами равен т.н. "dot product". Для нахождения освещенности плоскости , добавим "ambient light" кмировой координате источника света , помноженный на максимальную световую интенсивность :

Глобальные определения :

   WORD AmbientLight=20;
   #define MAXLIGHT 32
   static Vertex_t LightSource;
   WORD light;

//Подсчитаем световую интенсивность как "dot product" нормали и источника света
   light=MAXLIGHT*VEC_DotProduct(PlaneNormal,LightSource.World)+AmbientLight;
   if(light>MAXLIGHT)light=MAXLIGHT;
   if(light<1)light=1;

Инициализация светового источника :

   //Initialize Light Source
   LightSource.Local=P3D(0,0,0);
   MAT_Identity(matrix);
   TR_Translate(matrix, 10,10,10);
   TR_Rotate(matrix, 0,128-32,128);
   VEC_MultMatrix(&LightSource.Local,matrix,&LightSource.World);
   VEC_Normalize(&LightSource.World);

Содержание

Смотрите также: