Содержание
1. Введение 1.1. Предположения и соглашения 2. Основы 3D графики 2.1. Задание объектов и сцен 2.2. Проецирование 2.3. Матричные преобразования 2.4. Рисование одноцветного треугольника 2.5. Работа с произвольной камерой 3. Удаление невидимых частей 3.1. Отсечение нелицевых граней 3.2. Алгоритм художника 3.3. Z-буфер 3.4. Порталы 4. Текстурирование 4.1. Точное 4.2. Аффинное 4.3. Перспективно-корректное 4.4. Параболическое 5. Освещение 5.1. Модель освещения 5.2. Расчет нормали к объекту 5.3. Освещение по Ламберту 5.4. Освещение по Гуро 5.5. Освещение по Фонгу 6. Разное 6.1. Как совместить текстуру и освещение 6.1.1. 256-цветные режимы 6.1.2. 24/32-битные режимы 6.1.3. 15/16-битные режимы 6.2. Субпиксельная точность 6.3. Субтексельная точность 6.4. Поворот 3D вектора за шесть умножений 6.5. Билинейная фильтрация текстур 6.6. Алгоритм "бегущих кубиков" для полигонизации изоповерхностей
1. Введение
1.1. Предположения и соглашения
Предполагается, что читатель имеет какие-то основные сведения о компьютерной графике. Предполагается, что объект, состоящий из треугольных граней с наложенными на них текстурами при использовании освещения есть достаточно реалистичная и относительно простая модель. Также предполагается, что этот самый читатель способен воспринимать тот менторский тон, в котором, видимо, будет выдержана большая часть этого FAQ.
Рекомендуется наличие у читателя неких базовых познаний в стереометрии и линейной алгебре, а также программировании. Например, довольно сильно рекомендуется умение нарисовать одноцветный 2D треугольник.
Соглашения - следующие. Система 3D координат используется такая:
y
| z
| /
| /
-------O-------x
/ |
* |
/ |
Здесь буквами x, y, z обозначены положительные направления осей Ox, Oy и Oz соответственно. Также предполагается, что камера неподвижна и находится в точке (*) с координатами (0,0,-dist), ось зрения камеры направлена по оси Oz, а именно в точку (0,0,0) (т.е. camera target = (0,0,0)), ось Ox с точки зрения камеры направлена слева направо, ось Oy - снизу вверх, ось Oz - вглубь экрана. Размер экрана - xSize на ySize пикселов.
Проецирование на плоскость экрана в этом случае будет осуществляться по формулам
sx = xSize/2+x*dist/(z+dist),
sy = ySize/2-y*dist/(z+dist).
Здесь и далее используются обозначения
sx, sy - координаты проекции точки на экране,
x, y, z - 3D координаты точки,
dist - расстояние от камеры (она находится в точке (0,0,-dist)) до начала
координат,
u, v - координаты в текстуре (u - по горизонтали, v - по вертикали).
2. Основы 3D графики
2.1. Задание объектов и сцен
Покажем здесь достаточно распространенную схему задания 3D объектов и сцен.
Подобная схема, кстати, используется, в 3D Studio.
Каждая сцена представляет из себя следующее:
- набор объектов
- набор источников света
- набор текстур
- набор камер (хотя обычно используется одна)
Каждый объект задается следующим:
- набор вершин * вершина определяется своими 3D координатами и соответствующими ей координатами в текстуре
- набор граней * грань определяется тремя вершинами и текстурой (вообще говоря, не текстурой, а материалом: кроме текстуры могут быть заданы, например, коэффициенты рассеивания и отражения света)
- поведение объекта * то есть, расположение (то есть смещение, ось поворота, угол поворота, коэффициент масштабирования, и т.д.) в зависимости от номера кадра; обычно задается в нескольких ключевых точках и интерполируется между ними с помощью сплайнов
- положение
- ориентация (точка, в которую направлен этот источник, target)
- тип (фоновый/направленный/ненаправленный)
- цвет (обычно RGB)
- положение (location)
- направление (точнее, точкой, в которую направлена эта камера; target)
- угол зрения (FOV)
- угол поворота относительно своей оси (roll)
Каждый источник света задается следующим:
Каждая текстура представляет из себя прямоугольную 2D картинку, часто бывает фиксированных размеров (например, 64x64, 128x128, 256x256).
Каждая камера задается следующим:
2.2. Проецирование
Мы будем использовать только обычное перспективное проецирование на плоскость зрения "стандартной" камеры, определенной в пункте 1.1 (там же написаны и формулы проецирования). Случай произвольной камеры будет приводиться к случаю стандартной камеры.
2.3. Матричные преобразования
Вообще говоря, лучше всего немного почитать любую книжку по линейной алгебре. Здесь будет только краткий рассказ о 3D преобразованиях, о том, как их делать с помощью матриц, и о том, что же такое матрицы и как с ними работать.
Введем несколько терминов. n-мерный вектор, он же вектор размерности n, он же вектор размера n: упорядоченный набор n действительных чисел. Вообще говоря, практически то же самое, что и обычный 1D-массив. Матрица размера m на n (будет обозначаться как m*n, mxn): таблица размера m на n, в каждой клетке которой - действительное число. Это уже 2D-массив. Всего лишь. Вот пример матрицы 3x3:
[ 15 y*z 0.6 ]
[ 7 -3 91 ]
[ sin(x) 0.123 exp(t) ]
Займемся определением операций над векторами и матрицами. Вектор будем записывать в столбик и рассматривать его как матрицу размера n*1.
Скалярное произведение: определена для двух векторов одинаковых размеров. Результат - это число, он равен сумме произведений соответствующих элементов векторов. Пример:
[ 1 ] [ 4 ]
[ 2 ] * [ 5 ] = 1*4 + 2*5 + 3*6 = 32
[ 3 ] [ 6 ]
Векторное произведение: определена для (n-1) вектора одинакового размера n. Результат - вектор, причем, что интересно, перпендикулярный всем множителям. Результат меняется от перестановки мест множителей!!! Формально определяется как определитель матрицы, первая строка которой есть все базисные вектора, а все последующие - соответсвующий координаты всех множителей. Поскольку нам оно будет требоваться только для 3D пространства, мы определим векторное произведение двух 3D векторов явно:
[ Ax ] [ Bx ] | i j k | [ Ay*Bz-Az*Bx ]
AxB = [ Ay ] x [ By ] = | Ax Ay Az | = [ Az*Bx-Ax*Bz ]
[ Az ] [ Bz ] | Bx By Bz | [ Ax*By-Ay*Bx ]
Операция сложения двух матриц: определена для матриц одинаковых размеров. Каждый элемент суммы (то есть, каждое число в таблице) равняется сумме соответствующих элементов слагаемых-матриц. Пример:
[ 1 x 500 ] [ 8 a 3 ] [ 9 a+x 503 ]
[ 2 y 600 ] + [ 9 b 2 ] = [ 11 b+y 602 ]
[ 3 z 700 ] [ 10 c 1 ] [ 13 c+z 701 ]
Операция умножения матрицы на число: определена для любой матрицы и любого числа; каждый элемент результата равняется произведению соответствующего элемента матрицы-множителя и числа-множителя.
Операция умножения двух матриц: определена для двух матриц таких размеров a*b и c*d, что b = c.
Например, если b = c, но a != d, то при перестановке множителей операция вообще не определена. Так вот, результатом матрицы A размером a*b и матрицы B размером b*d будет матрица C размером a*d, в которой элемент, стоящий в строке i и столбце j равен произведению строки i матрицы A на столбец j матрицы B. Произведение строки на столбец определяется как сумма произведений соответствующих элементов строки и столбца. Чтобы было хоть чуть-чуть понятно, пример умножения строки на столбец (они должны быть равной длины, кстати; поэтому и такие ограничения на размеры матриц):
[ 4 ]
[ 1 2 3 ] * [ 5 ] = 1*4 + 2*5 + 3*6 = 32
[ 6 ]
А чтобы перемножить две матрицы, надо эту операцию проделать для каждого элемента. Вот пример:
[ 1 2 3 ] [ 0 3 ] [ 1*0+2*1+3*2 1*3+2*4+3*5 ]
[ 4 5 6 ] * [ 1 4 ] = [ 4*0+5*1+6*2 4*3+5*4+6*5 ] = ...
[ 7 8 9 ] [ 2 5 ] [ 7*0+8*1+9*2 7*3+8*4+9*5 ]
Умножение и сложение матриц обладают почти тем же набором свойств, что и обычные числа, хотя некоторые привычные свойства не выполняются (например, A*B != B*A); нам на самом деле понадобится знать, что произведение вида A*B*C*D*... не зависит от того, как расставить скобки. Или, если угодно, что
A*(B*C) = (A*B)*C.
Теперь забудем об этом на некоторое время и перейдем к преобразованиям. Любое движение (то есть преобразование пространства, сохраняющее расстояние между точками) в трехмерном пространстве, согласно теореме Шаля, может быть представлено в виде суперпозиции поворота и параллельного переноса, то есть последовательного выполнения поворота и параллельного переноса. Именно поэтому основная часть информация о поведении объекта - это его смещение, ось поворота и угол поворота. И именно поэтому нам достаточно знать, как сделать два преобразования - перенос и поворот.
Перенос точки (кстати, точки будут также рассматриваться как вектора с началом в начале координат и концом в собственно точке) с координатами (x,y,z) на вектор (dx,dy,dz) делается простым сложением всех координат. То есть результат - это (x+dx,y+dy,z+dz). Как бы сложили вектор-точку и вектор-перенос.
Поворот - занятие уже более интересное. Но тоже простое. Рассмотрим для примера поворот точки (x,y,z) относительно оси z. В этом случае z не меняется вообще, а (x,y) меняются также, как и при 2D повороте относительно начала координат.
Посмотрим, какие будут координаты у точки A' - результата поворота A(x,y) на угол alpha.
y * A' - результат
| - поворота
| -
| -
| - * A - исходная
| - --- точка
| - ---
| - ---
| ---
----O-----------------------x
|
|
Пусть r = sqrt(x*x+y*y). Пусть угол AOx равен phi, тогда из рисунка видно,
что cos(phi) = x/r, sin(phi) = y/r. Угол A'OA равен по условию alpha.
Отсюда
x' = r*cos(alpha+phi) = r*(cos(alpha)*cos(phi)-sin(alpha)*sin(phi)) =
= (r*cos(alpha))*cos(phi)-(r*sin(alpha))*sin(phi) =
= x*cos(phi)-y*sin(phi)
y' = r*sin(alpha+phi) = r*(cos(alpha)*sin(phi)+sin(alpha)*cos(phi)) =
= (r*cos(alpha))*sin(phi)+(r*sin(alpha))*cos(phi) =
= x*sin(phi)+x*cos(phi)
Для трехмерного случая, таким образом
x' = x*cos(phi)-y*sin(phi)
y' = x*sin(phi)+x*cos(phi)
z' = z
Аналогичные формулы получатся и для других осей поворота (то есть Ox, Oy). Поворот относительно произвольной оси, проходящей через начало координат, можно сделать с помощью этих поворотов - сделать поворот относительно Ox так, чтобы ось поворота стала перпендикулярна Oy, затем поворот относительно Oy так, чтобы ось поворота совпала с Oz, сделать собственно поворот, а затем обратные повороты относительно Oy и Ox. Можно даже вывести формулы для такого поворота и убедиться в том, что они очень громоздкие.
Вспомним про матрицы и вектора и посмотрим внимательно на выведенные формулы для поворота. Можно заметить, что
[ x' ] = [ cos(phi) -sin(phi) 0 ] [ x ]
[ y' ] = [ sin(phi) cos(phi) 0 ] [ y ]
[ z' ] = [ 0 0 1 ] [ z ]
То есть поворот на угол phi задается одной и той же матрицей, и с помощью этой матрицы (умножая ее на вектор-точку) можно получить координаты повернутой точки. Пока никакого выигрыша не видно - здесь умножение матрицы на вектор требует больше операций, чем расчет x' и y' по формулам.
Удобство матриц для нас заключается как раз в свойстве A*(B*C) = (A*B)*C. Пусть мы делаем несколько поворотов подряд, например 5 (как раз столько, сколько надо для поворота относительно произвольной оси), и пусть они задаюся матрицами A, B, C, D, E (A - матрица самого первого поворота, E - последнего). Тогда для вектора p мы получаем
p' = E*(D*(C*(B*(A*p)))) = E*D*C*B*A*p = (E*D*C*B*A)*p =
= (E*(D*(C*(B*A))))*p = T*p,
где
T = (E*(D*(C*(B*A))))
матрица преобразования, являющегося комбинацией пяти поворотов. Посчитав один раз эту матрицу, можно в дальнейшем без проблем применить довольно сложное преобразование из пяти поворотов к любому вектору с помощью всего одного умножения матрицы на вектор.
Таким образом, можно задать любой поворот матрицей, и любая комбинация поворотов также будет задаваться матрицей, которую можно довольно легко посчитать. Но есть еще параллельный перенос, есть еще масштабирование. Что делать с ними?
На самом деле, эти преобразования тоже легко записываются в виде матриц. Только вместо матриц 3x3 и 3-мерный векторов используются так называемые однородные 4-мерные координаты и матрицы 4x4. При этом вместо векторов вида
[ x ]
[ y ]
[ z ]
используются вектора вида
[ x ]
[ y ]
[ z ]
[ 1 ]
а вместо произвольных матриц 3x3 используются матрицы 4x4 такого вида:
[ a b c d ]
[ e f g h ]
[ i j k l ]
[ 0 0 0 1 ]
В случае, если d = h = l = 0, получается то же самое, что и для матрицы 3x3. Матрица параллельного переноса теперь определяется как
[ 1 0 0 dx ]
[ 0 1 0 dy ]
[ 0 0 1 dz ]
[ 0 0 0 1 ]
Матрицу масштабирования можно определить и для матриц 3x3, и для матриц 4x4:
[ kx 0 0 ] [ kx 0 0 0 ]
[ 0 ky 0 ] или [ 0 ky 0 0 ]
[ 0 0 kz ] [ 0 0 kz 0 ]
[ 0 0 0 1 ]
где kx, ky, kz - коэффициенты масштабирования по соответствующим осям.
Таким образом, получаем следующее. Любое нужное нам преобразование пространства можно задать матрицей 4x4 определенной структуры, разной для разных преобразований. Результат последовательного выполнений нескольких преобразований совпадает с результатом одного преобразования T, которое также задается матрицей 4x4, вычисляемой как произведение матриц всех этих преобразований. Важен порядок умножения, так как A*B != B*A. Результат применения преобразования T к вектору [ x y z ] считается как результат умножения матрицы T на вектор [ x y z 1 ]. Вот и все.
Осталось только на примере показать, почему A*B != B*A. Пусть A - матрица переноса, B - поворота. Если мы сначала перенесем объект, а потом повернем (это будет B*A), получим далеко не то, что будет, если сначала объект повернуть, а потом перенести (это уже A*B).
2.4. Рисование одноцветного треугольника
Без понимания того, как рисовать залитый одним цветом треугольник, дальше лезть в 3D графику явно не стоит. Поэтому вот объяснение.
Возьмем любой треугольник. Его изображение на экране - набор горизонтальных отрезков, причем из-за того, что треугольник - фигура выпуклая, каждой строке экрана соответствует не более одного отрезка. Поэтому достаточно пройтись по всем строкам экрана, в которых есть треугольник (то есть, от минимального до максимального значения y для вершин треугоольника) и нарисовать соответствующие горизонтальные отрезки. См. рис.
----------------------------
| A |
| **** |
| ******* |
| ********** |
| B************ |
| *********** |
| ********* |
|-----------@@@@@@@--------|
| ***** |
| *** |
| C |
----------------------------
Отсортируем вершины так, чтобы вершина A была верхней, C - нижней, тогда у нас min_y = A.y, max_y = C.y, и нам надо пройтись по всем линиям от min_y до max_y. Рассмотрим какую-то линию sy, A.y <= sy <= C.y. Если sy < B.y, то она пересекает стороны AB и AC; если sy >= B.y - то стороны BC и AC. Мы знаем координаты всех вершин, поэтому мы можем написать уравнения сторон и найти пересечение нужной стороны с прямой y = sy. Получим два конца отрезка. Так как мы не знаем, какой из них левый, а какой правый, сравним их координаты по x и обменяем значения, если надо. Рисуем этот отрезок, повторяем процедуру для каждой строки - и вуаля, трегуольник нарисован.
Остановимся более подробно на нахождении пересечения прямой y = sy (текущей строки) и стороны треугольника, например AB. Напишем уравнение прямой AB в форме x = k*y+b:
x = A.x+(y-A.y)*(B.x-A.x)/(B.y-A.y)
Подставляем сюда известное значение y = sy для текущей прямой:
x = A.x+(sy-A.y)*(B.x-A.x)/(B.y-A.y)
Вот, в общем-то, и все. Для других сторон пересечение ищется совершенно точно так же. А вот и пример кода.
// ...
// здесь сортируем вершины (A,B,C)
// ...
for (sy = A.y; sy <= C.y; sy++) {
x1 = A.x + (sy - A.y) * (C.x - A.x) / (C.y - A.y);
if (sy < B.y)
x2 = A.x + (sy - A.y) * (B.x - A.x) / (B.y - A.y);
else
x2 = B.x + (sy - B.y) * (C.x - B.x) / (C.y - B.y);
if (x1 > x2) { tmp = x1; x1 = x2; x2 = tmp; }
drawHorizontalLine(sy, x1, x2);
}
// ...
Надо, правда, защититься от случая, когда B.y = C.y - в этом (и только этом, потому как если C.y = A.y, то треугольник пустой и рисовать его не стоит, или можно рисовать горизонтальную линию; а если B.y = A.y, то sy >= A.y и до деления на B.y - A.y не дойдет) случае произойдет попытка деления на ноль. Код изменится совсем чуть-чуть:
// ...
// здесь сортируем вершины (A,B,C)
// ...
for (sy = A.y; sy <= C.y; sy++) {
x1 = A.x + (sy - A.y) * (C.x - A.x) / (C.y - A.y);
if (sy < B.y)
x2 = A.x + (sy - A.y) * (B.x - A.x) / (B.y - A.y);
else {
if (C.y == B.y)
x2 = B.x;
else
x2 = B.x + (sy - B.y) * (C.x - B.x) / (C.y - B.y);
}
if (x1 > x2) { tmp = x1; x1 = x2; x2 = tmp; }
drawHorizontalLine(sy, x1, x2);
}
// ...
Вот и все. Ну, горизонтальную линию, надеюсь, нарисовать сумеют все желающие.
2.5. Работа с произвольной камерой
Рассмотрим любую камеру как точку - центр проецирования и экран - плоский прямоугольник в 3D пространстве, на плоскость которого идет проецирование. Наша стандартная камера, например, задается точкой (0,0,-dist) и экраном с вершинами (-xSize/2,ySize/2), ..., (xSize/2,-ySize/2). Можно задать эту систему тремя векторами, задающими с точки зрения камеры направления вперед, вправо и вверх; вектор "вперед" соединияет центр проецирования и центр экрана, вектор "вправо" соединияет центр экрана и правую его границу, вектор "вверх", соответственно, центр экрана и верхнюю его границу. Обозначим эти вектора как p, q и r соответственно, а центр проецированя за s. Вот пример для стандартной камеры.
y
| z
| /
=====^=======+
@ / + <--------- экран
@ <-----+----------- r, "вверх"
@ / +
------------------O@@@@@@@>-------x
@ | ^-+----------- q, "вправо"
@ <---------+----------- p, "вперед"
@ | +
@ ===========+
* |
/ |
|
|
Здесь (для стандартной камеры; обозначим ее вектора как Sp, Sq, Sr, Ss)
Sp = p = (0,0,dist)
Sq = q = (xSize/2,0,0)
Sr = r = (ySize/2,0,0)
Ss = s = (0,0,-dist)
Любые три взаимно перпендикулярных вектора и точка - центр координат задают в 3D пространстве систему координат. Так что объект мы можем рассматривать в системе обычных координат (x,y,z), в системе координат стандартной камеры (Sp,Sq,Sr) или в системе (p,q,r), соответствующей какой-то произвольной камере. В любом случае, если (a,b,c) - координаты точки в системе координат камеры (точнее, в системе координат с центром в точке s и базисом (p,q,r)), то координаты проекции точки на экране равны
screenX = xSize/2 + xSize/2 * a/c
screenY = ySize/2 - ySize/2 * b/c
В случае стандартной камеры переход от обычной системы координат к системе координат камеры очевиден:
a = x / (xSize/2)
b = y / (ySize/2)
c = (z + dist) / dist
Подставив это в формулы для screenX, screenY, получим как раз те самые формулы для проекции на стандартную камеру.
Поскольку со стандартной камерой нам достаточно удобно и понятно работать, для произвольной камеры мы должны сделаеть такое преобразование пространства, что как бы совместит произвольную камеру и стандартную камеру. То есть, такое преобразование, что вектора p, q, r перейдут в Sp, Sq, Sr, а точка s в точку Ss.
Посчитаем матрицу для *обратного* преобразования; оно должно переводить Sp, Sq, Sr, Ss в p, q, r, s. Преобразование, переводящее Ss в s (и наоборот) - это обычный паралелльный перенос; остается написать преобразование перевода Sp, Sq, Sr в p, q, r. Пусть у нас есть координаты p, q, r в системе координат (x,y,z):
p = (px,py,pz)
q = (qx,qy,qz)
r = (rx,ry,rz)
Для Sp, Sq, Sr координаты (в этой же системе) известны и равны следующему:
Sp = (0,0,dist)
Sq = (xSize/2,0,0)
Sr = (ySize/2,0,0)
Пусть T - искомая матрица перевода,
[ a b c ]
T = [ d e f ], a..i - какие-то неизвестные.
[ g h i ]
Поскольку T переводит Sp, Sq, Sr в p, q, r; то есть
p = T*Sp
q = T*Sq
r = T*Sr
то, подставляя, например, p и Sp, получаем:
[ px ] [ a b c ] [ 0 ] [ c*dist ]
[ py ] = [ d e f ] [ 0 ] = [ f*dist ], откуда
[ pz ] [ g h i ] [ dist ] [ i*dist ]
c = px/dist
f = py/dist
i = pz/dist.
Аналогично находим все остальные элементы матрицы T:
[ qx*2/xSize rx*2/ySize px/dist ]
T = [ qy*2/xSize ry*2/ySize py/dist ]
[ qz*2/xSize rz*2/ySize pz/dist ]
Но нас интересует обратное к этому преобразование. Оно задается обратной матрицей к T, то есть такой матрицей T1, что
[ 1 0 0 ]
T * T1 = T1 * T = [ 0 1 0 ]
[ 0 0 1 ]
Обратная матрица, вообще говоря, существует далеко не всегда, да и вычисление ее в общем случае - достаточно неприятная задача. Однако в данном случае из-за специального вида матрицы T (конкретнее, из-за того, что T - ортогональная матрица) она не только всегда существует, но и считается очень просто:
[ qx*2/xSize rx*2/ySize px/dist ] [ qx1 rx1 px1 ]
T = [ qy*2/xSize ry*2/ySize py/dist ] = [ qy1 ry1 py1 ]
[ qz*2/xSize rz*2/ySize pz/dist ] [ qz1 rz1 pz1 ]
[ qx1/lq qy1/lq qz1/lq ]
T1 = [ rx1/lr ry1/lr rz1/lr ]
[ px1/lp py1/lp pz1/lp ]
где
lp = px1*px1 + py1*py1 + pz1*pz1
lq = qx1*qx1 + qy1*qy1 + qz1*qz1
lr = rx1*rx1 + ry1*ry1 + rz1*rz1
Сделав сначала параллельный перенос, совмещающий s и Ss, а потом полученное преобразование, как раз и получим преобразование, переводящее произвольную камеру в стандартную.
Теперь надо выяснить, как, собственно посчитать координаты p, q, r через имеющиеся у нас характеристики: положение, направление, угол зрения и угол поворота. 3D Studio (и мы вслед за ней) рассчитывает эти вектора по такому алгоритму:
1. Считаем p = target - location
2. Если p.x == 0 и p.z == 0, то q = (0, 0, 1); иначе q = (p.z, 0, -p.x)
3. Считаем r = crossProduct(p, q) - векторное произведение p на q
4. Считаем lp = length(p) - длина p
5. Приводим r и q к длине 2*lp*tan(FOV/2)
Здесь мы не учитываем поворот камеры вокруг своей оси, его удобнее сделать после перехода к стандартной камере - в этом случае получаем обычный поворот относительно оси z на угол roll.
Таким образом, окончательная матрица перевода должна представлять собой произведение матрицы параллельного переноса, матрицы T и матрицы поворота вокруг оси z на угол roll:
FinalCameraMatrix = RollMatrix * T * MoveMatrix
Расчет матриц RollMatrix и MoveMatrix очевиден.
3. Удаление невидимых частей
3.1. Отсечение нелицевых граней
Пусть у нас есть объект, внутри которого камера заведомо не окажется. Обычно такие объекты составляют большую часть или всю сцену. Тогда для каждой грани мы можем увидеть только одну ее сторону - лицевую. Грань - плоскость, она делит все 3D пространство на два полупространства. Так вот, лицевую сторону видно только из одного полупространства, из того, в которое "смотрит" нормаль к этой грани, направленная *из* объекта. Проверив, в какое полупространство попадает камера, можно сразу определить, является ли грань лицевой (то есть, может ли камера увидеть лицевую сторону этой грани) и надо ли ее рисовать.
Пусть грань имеет вершины v1, v2, v3 и нормаль (nx,ny,nz). Тогда уравнение плоскости, в которой она лежит, будет иметь вид
nx*x+ny*y+nz*z+d = 0.
d находим из того факта, что v1 в плоскости лежит:
nx*v1.x+ny*v1.y+nz*v1.z+d = 0,
d = -(nx*v1.x+ny*v1.y+nz*v1.z).
Функция nx*x+ny*y+nz*z+d принимает положительные значение по одну сторону от плоскости, отрицательные по другую и равна нулю на самой плоскости. Точка (v1.x+nx,v1.y+ny,v1.z+nz) лежит, очевидно, в том полупространстве, откуда грань видно. Значение функции (назовем ее функцией видимости) в этой точке равно
nx*(v1.x+nx)+ny*(v1.y+ny)+nz*(v1.z+nz)+d =
nx*(v1.x+nx)+ny*(v1.y+ny)+nz*(v1.z+nz)-(nx*v1.x+ny*v1.y+nz*v1.z) =
nx*nx+ny*ny+nz*nz > 0.
Таким образом, если значение функции в какой-то точке больше нуля, то грань из этой точки потенциально видна. Нас интересует видимость грани из нашей камеры, а камера у нас зафиксирована в точке (0,0,-dist). Таким образом, получаем, что грани, для которых
-nz*dist-(nx*v1.x+ny*v1.y+nz*v1.z) < 0,
или, что равносильно,
nz*dist+nx*v1.x+ny*v1.y+nz*v1.z > 0,
заведомо не видны и время на их обработку и отрисовку тратить не стоит.
Отдельный вопрос - как считать нормали к граням. Точнее, как выбрать одну из двух нормалей, которая будет смотреть из объекта. Обычно эта проблема решается еще на этапе построения 3D моделей - например, пакет 3D Studio записывает вершины граней в порядке A-B-C так, чтобы векторное произведение BA*CA и было нормалью. Еще один способ - выбрать внутренную точку для объекта (либо вручную, либо взять его центр тяжести, либо еще как-нибудь - методов может быть придумано сколько угодно) и использовать ее: если для этой точки функция видимости положительна, то есть грань якобы видна, то надо поменять знак nx, ny и nz.
Осталось отметить, что для выпуклых объектов этот метод полностью решает задачу об удалении невидимых частей. Для невыпуклых же он позволяет быстро и просто сократить число граней, подлежащих дальнейшей проверке на видимость и собственно отрисовке.
3.2. Алгоритм художника
Пусть имеется некий набор граней (т.е. сцена), который требуется нарисовать. Отсортируем грани по удаленности от камеры и отрисуем все грани, начиная с самых удаленных. Довольно распространенная характеристика удаленности для грани ABC - это среднее значение z, mid_z = (A.z+B.z+C.z)/3. Вот и весь алгоритм. Просто, и обычно достаточно быстро.
Существует, правда, несколько проблем.
Во-первых, при некотором расположении граней этот алгоритм вообще не может дать правильного результата - в каком порядке грани не рисуй, получится неправильно. Вот стандартный пример.
+----+ +----+
| | | |
+---| |--------------
| | | |
+---| |--------------
| | | |
| | | |
+-------------| |---+
| | | |
+-------------| |---+
| | | |
+----+ +----+
Во-вторых, при некотором расположении граней и использовании среднего значения z как характеристики удаленности алгоритм тоже дает неправильный результат. Пример чуть ниже. В этом случае горизонтальную грань надо отрисовывать второй, но по среднему значению z она лежит дальше и таким образом получается, что ее надо отрисовывать первой. Возможные пути решения этой проблемы - какие-то изменения характеристики удаленности, либо моделирование, не вызывающее таких ситуаций.
y
| |
| | <- грань #1 (параллельна Oxy)
| |
| --------------------------------------
| ^^^ грань #2 (параллельна Oxz)
|
|
|
|
--*------O-------------------------------------------------z
^-- камера
И наконец, при использовании этого алгоритма отрисовываются вообще все грани сцены, и при большом количестве загораживающих друг друга граней мы будем тратить большую часть времени на отрисовку невидимых в конечном итоге частей. То есть совершенно впустую. В таком случае целесообразно использовать какие-то другие методы (например BSP и PVS, порталы, и т.д).
3.3. Z-буфер
Заведем буфер (собственно z-буфер) размером с экран, и забьем его каким-то большим числом, настолько большим, что координаты z для точек сцены заведомо меньше. Например, если z - fixedpoint 16:16, то можно использовать максимально возможное значение, то есть 0x7FFFFFFF. Для каждой рисуемой точки считаем значение z; если оно больше, чем значение в z-буфере (точка закрыта какой-то другой точкой), или меньше, чем -dist (точка находится за камерой), то переходим к следующей точке. Если меньше, то рисуем точку на экране (или в видеобуфере), а в z-буфер записываем текущее значение z. Вот кусок кода для лучшего понимания:
// ...
if (z > -dist && z < zbuffer[screenX][screenY]) {
screen[screenX][screenY] = color;
zbuffer[screenX][screenY] = z;
}
// ...
Имеет смысл считать значения не z, а z1 = 1/(z+dist), так как эта величина изменяется по грани линейно, и линейная интерполяция дает точные результаты (более подробно это описано в 4.3). Тогда условия чуть изменяются - точка загорожена другой, если значение z1 *меньше* значения в z-буфере; и точка находится за камерой, если z1 < 0. Код тоже чуть изменится:
// ...
if (z1 > 0 && z1 > zbuffer[screenX][screenY]) {
screen[screenX][screenY] = color;
zbuffer[screenX][screenY] = z1;
}
// ...
Вот и все.
Это самый простой метод удаления невидимых частей, причем всегда дающий полностью правильные результаты. Единственная проблема - скорость, z-буфер, как самый простой, является и самым медленным методом.
3.4. Порталы
Этот метод предназначен для т.н. indoor environments - в буквальном переводе "комнатная среда"; это как бы внутренность какого-то здания. В качестве примера могут служить Doom, Quake.
В этом случае сцена составляется из комнат. Комнаты связаны друг с другом какими-то проходами - дверями, окнами, чем угодно. Hа каждый такой проход вешается обычный полигон, называемый порталом. Для отрисовки сцены достаточно нарисовать текущую комнату вместе с тем, что из нее видно. Но то, что из нее видно, может быть видно *только* через порталы.
Таким образом, для отрисовки всей сцены мы просто берем и рисуем текущую комнату, используя любые методы удаления невидимых частей (например, алгоритм художника и отсечение нелицевых граней). А если в процессе отрисовки нам встречается полигон, помеченный как портал, то мы рисуем комнату, которую видно через этот портал, причем отсекаем ее уже по проекции этого портала.
Как обычно, кусок кода для уяснения:
// ...
void drawRoom(int roomID, polygon *clipArea) {
int i;
for (i = 0; i < rooms[roomID].numFaces; i++)
drawFace(&rooms[roomID].face[i], clipArea);
}
void drawFace(polygon *face, polygon *clipArea) {
if (face->isPortal) {
if (isVisible(face, clipArea)
drawRoom(face->destinationRoom, clipArea);
} else drawClippedPolygon(face, clipArea);
}
// ...
drawRoom(currentRoom, fullScreen);
// ...
Здесь функция isVisible(face, clipArea) возвращает false, если грань face и область отсечения clipArea не пересекаются; fullScreen - полигон размером с экран (или любую нужную область отрисовки).
В природе есть Alpha 2 - готовый engine, использующий порталы и написанный на TMT Pascal, причем, что приятно, доступный в исходниках. Скачать его можно с http://www.geocities.com/CapeCanaveral/5402.
Осталось упомянуть, что несмотря на то, что в изложенной форме метод и требует возможности отрисовки произвольной грани, отсеченной в произвольную грань, можно обойтись и отсечением по bounding box пpоекции поpтала.
4. Текстурирование
4.1. Точное
Задача текстурирования формулируется таким образом: есть грань - согласно предположениям, треугольная - с наложенной на нее текстурой. То есть каждая точка грани окрашена цветом соответствующей ей точки в текстуре. Текстура накладывается линейным образом. Есть точка экрана с координатами на экране (sx,sy), принадлежащая проекции грани. Требуется найти ее цвет, то есть цвет соответствующей этой точке экрана точки текстуры. А для этого надо найти координаты текстуры для этой точки - точнее, для той точки, проекцией которой на экран является наша (sx,sy).
Пусть вершины грани есть точки A(Ax,Ay,Az), B(Bx,By,Bz) и C(Cx,Cy,Cz), а соответствующие им точки текстуры - (Au,Av), (Bu,Bv) и (Cu,Cv). Найдем координаты текстуры для точки, проекцией которой является (sx,sy).
Для точек (x,y,z), проекцией которых является (sx,sy) имеем:
sx = xSize/2+x*dist/(z+dist),
sy = ySize/2-y*dist/(z+dist).
Для упрощения формул будем использовать обозначения
i = sx-XSize/2,
j = YSize/2-sy,
Z = z+dist.
Тогда эти формулы примут видM
i = x*dist/Z,
j = y*dist/Z,
или, что равносильно:
i*Z = x*dist,
j*Z = y*dist.
Рассмотрим точку D, принадлежащую грани. Для нее D = A+a*(B-A)+b*(C-A), так как она лежит в грани. D однозначно задается парой (a,b). Для нее координаты текстуры (из того, что текстура накладывается линейно) будут такие:
Du = Au+a*(Bu-Au)+b*(Cu-Au),
Dv = Av+a*(Bv-Av)+b*(Cv-Av).
Пусть проекция D на экран как раз и имеет координаты (sx,sy), тогда для нее выполнены написанные выше соотношения:
i*DZ = Dx*dist,
j*DZ = Dy*dist.
Подставим сюда координаты D, выраженные через координаты A, B и неизвестные коэффициенты a, b:
i*(Az+a*(Bz-Az)+b*(Cz-Az)+dist) = dist*(Ax+a*(Bx-Ax)+b*(Cx-Ax)),
j*(Az+a*(Bz-Az)+b*(Cz-Az)+dist) = dist*(Ay+a*(By-Ay)+b*(Cy-Ay)),
т.к. мы договорились, что AZ = Az+dist, и, кстати, поэтому BZ-AZ = Bz-Az, это можно чуть укоротить:
i*(AZ+a*(BZ-AZ)+b*(CZ-AZ)) = dist*(Ax+a*(Bx-Ax)+b*(Cx-Ax)),
j*(AZ+a*(BZ-AZ)+b*(CZ-AZ)) = dist*(Ay+a*(By-Ay)+b*(Cy-Ay)).
Выделим a и b:
a*(i*(BZ-AZ)-dist*(Bx-Ax))+b*(i*(CZ-AZ)-dist*(Cx-Ax)) = dist*Ax-i*AZ,
a*(j*(BZ-AZ)-dist*(By-Ay))+b*(j*(CZ-AZ)-dist*(Cy-Ay)) = dist*Ay-j*AZ.
Получили систему двух линейных уравнений, из которой можно найти a и b, а по ним немедленно вычисляются u и v. Введем вектор M = BA (Mx = Bx-Ax и т.д.) и вектор N = CA, обозначим d = dist (все это для краткости записи). Тогда эти два уравнения перепишутся в виде:
a*(i*Mz-d*Mx)+b*(i*Nz-d*Nx) = d*Ax-i*AZ,
a*(j*Mz-d*My)+b*(j*Nz-d*Ny) = d*Ay-j*AZ,
и решая систему (например, по формулам Крамера) получаем:
(dAx-iAZ)(jNz-dNy)-(dAy-jAZ)(iNz-dNx)
a = ------------------------------------- =
(iMz-dMx)(jNz-dNy)-(iNz-dNx)(jMz-dMy)
djAxNz-ddAxNy-ijAZNz+idAZNy-diAyNz+ddAyNx+ijNzAZ-djAZNx
= ------------------------------------------------------- =
ijMzNz-idMzNy-djMxNz+ddMxNy-ijMzNz+idNzMy+djNxMz-ddNxMy
dj(AxNz-AZNx)+dd(AyNx-AxNy)+id(AZNy-AyNz)
= ----------------------------------------- =
id(MyNz-MzNy)+dj(NxMz-MxNz)+dd(MxNy-NxMy)
i(AZNy-AyNz)+j(AxNz-AZNx)+d(AyNx-AxNy)
= --------------------------------------
i(MyNz-MzNy)+j(NxMz-MxNz)+d(MxNy-NxMy)
аналогично получаем
i(AZMy-AyMz)+j(AxMz-AZMx)+d(AyMx-AxMy)
b = --------------------------------------
i(MyNz-MzNy)+j(NxMz-MxNz)+d(MxNy-NxMy)
Знаки умножения здесь везде опущены, иначе читать это все совсем невозможно.
Осталось немного, подставить Mx = Bx-Ax и так далее, а потом подставить эти a и b в формулы для Du и Dv. В процессе подстановок что-то сократится, что-то упростится, но формулы для Du и Dv все равно получатся длиной в несколько строк.
Но из этих формул видно кое-что интересное. А именно, то, что
u = (C1*sx+C2*sy+C3) / (C4*sx+C5*sy+C6),
v = (C7*sx+C8*sy+C9) / (C4*sx+C5*sy+C6),
причем C1, ..., C9 - просто какие-то коэффициенты, зависящие от грани. То есть можно посчитать эти коэффициенты один раз на грань, а потом считать u, v по написанным пару строк выше формулам. Но все равно получается как минимум одно деление на точку. Это слишком медленно. Поэтому все методы текстурирования на самом деле используют приближенные вычисления, хотя и именно по этим формулам.
Стоит отметить тот факт, что на самом деле здесь не грань не обязательно должна быть треугольной - можно взять любые три вершины многоугольной грани. То же самое справедливо и для остальных описанных методов текстурирования.
4.2. Аффинное
Этот метод текстурирования основан на приближении u, v линейным функциями. Итак, пусть u - линейная функция, u = k1*sx+k2*sy+k3. Можно посчитать k1, k2, k3 исходя из того, что хотя бы в вершинах грани u должно совпадать с точным значением - это даст нам три уравнения, из которых быстро и просто находятся эти коэффициенты, и потом считать u по этой формуле. Но это все равно медленно - два умножения на пиксел.
Будем рисовать грань по строкам - это общепринято, довольно просто, и не доводит до умопомешательства кэш-память процессора. Вершины грани заранее отсортируем по sy (например, A.sy <= B.sy <= C.sy). Для каждой строки можно посчитать начальное значение x, u (также, как и для x, ведь u по любой прямой меняется тоже линейно), а также длину этой строки.
A
/ \
/ \
*--------*
/ B
/ /
/ /
/ /
/ /
/ /
//
C
В нарисованном случае, например,
x_start = A.sx+(current_sy-A.sy)*(C.sx-A.sx)/(C.sy-A.sy),
u_start = A.u+(current_sy-A.sy)*(C.u-A.u)/(C.sy-A.sy),
x_end = A.sx+(current_sy-B.sy)*(B.sx-A.sx)/(B.sy-A.sy).
length = x_end - x_start.
Какие вершины использовать в этих формулах - это уже проблемы рисования треугольника, а не текстурирования. Лично я просто храню x_start, x_end, u_start, на каждом переходе вниз на строчку прибавляю к
delta_x_start = (C.sx-A.sx)/(C.sy-A.sy),
и аналогично высчитываемые приращения для x_end, u_start. Вот только надо аккуратно следить за тем, какая сторона правая, какая - левая, и на каком мы сейчас промежутке находимся - то ли AB, то ли BC, и соответственно изменять приращения. Впрочем, все это - уже обыкновенное рисование треугольника. В примерах просто сделано решение "в лоб" - проверяем, какой участок - AB или BC - пересекает текущая строка, считаем x/u/v на обоих концах, считаем длину строки и берем соответствующие левому концу (то есть меньшему x) значения u и v.
Так вот. Посчитали начало строки, длину строки, u в начале строки. Осталось заметить, что раз уж u = k1*sx+k2*sy+k3, то при переходе к следующему пикселу строки у нас u изменяется на k1 (так же известный как du/dsx). Это число и надо как-то посчитать. Например, так:
A
/ \
/ \ x_start = A.sx+(B.sy-A.sy)*(C.sx-A.sx)/(C.sy-A.sy),
/ \ x_end = B.sx,
*-----------B u_start = A.u+(B.sy-A.sy)*(C.u-A.u)/(C.sy-A.sy),
/ / u_end = B.u,
/ / du_dsx = (u_start-u_end)/(x_start-x_end).
/ /
/ /
/ /
//
C
du/dsx - просто число, оно не меняется на всем треугольнике, поэтому просто считаем его там, где удобно, и берем посчитанное значение.
v (из тех же соображений) считается абсолютно точно так же, надо только во всех приведенных формулах u заменить на v, и все.
Теперь осталось только взять и нарисовать эту строку:
// ...
u = u_start;
v = v_start;
for (current_sx = x_start; current_sx <= x_end; current_sx++) {
putpixel(current_sx, current_sy, texture[(int)v][(int)u]);
u += du_dsx;
v += dv_dsx;
}
// ...
Пройдясь по всем строкам грани - т.е. пробежавшись current_sy по значениям от A.sy до C.sy (вершины отсортированы!), получим текстурированную грань.
Voila!
4.3. Перспективно-корректное
Этот метод основан на приближении u, v кусочно-линейным функциями. Кратко говоря, при рисовании каждая строка разбивается на куски (обычно несколько кусков длиной 8/16/32 пикселов и один оставшийся произвольной длины), в начале и конце каждого куска считаются точные значения u, v, а на куске они интерполируется линейно.
Точные значения u и v, в принципе, можно считать по формулам из 4.1, но обычно используют более простой путь. Он основан на том факте, что значения 1/Z, u/Z и v/Z зависят от sx, sy ЛИНЕЙНО. Доказательство этого факта пока опущено. Таким образом, достаточно для каждой вершины посчитать 1/Z, u/Z, v/Z и линейно их интерполировать - точно так же, как интерполируются u и v в 4.2. Причем, так как эти значения зависят от sx, sy строго линейно, то интерполяция дает не сильно приближенные результаты, а абсолютно точные!
Сами же точные значения u, v считаются как
u = (u/Z) / (1/Z),
v = (v/Z) / (1/Z).
Дальше все становится совсем просто. При рисовании треугольника, на ребрах интерполируем не u и v, как в 4.2, а 1/Z, u/Z, v/Z. Кроме того, заранее считаем d(u/Z)/dsx, d(v/Z)/dsx, d(1/Z)/dsx (то есть, изменений этих самых u/Z, v/Z, 1/Z соотвествующих шагу по dsx на 1) так, как считали du/dsx - это будет нужно для быстрого вычисления точных значений u, v. Каждую линию рисуем кусками по 8/16/32 пикселов (на самом деле, кусками любой длины; просто если длина - степень двойки, то при вычислении du/dx и dv/dx для текущего куска можно деление на длину куска заменить сдвигом вправо) и, если надо, рисуем оставшийся хвостик. Для расчета точных значения u, v в конце каждого куска пользуемся посчитанными (ага!) значения d(u/Z)/dsx, d(v/Z)/dsx, d(1/Z)/dsx; раз значения u/Z, v/Z, 1/Z в начале куска известны, меняются они линейно и длина куска известна (либо 16 пикселов, либо длина остатка), то в конце куска они считаются все это до боли просто:
uZ_b = uZ_a + length * duZ_dsx; // расчет u/Z, v/Z, 1/Z в конце куска
vZ_b = vZ_a + length * dvZ_dsx;
Z1_b = Z1_a + length * dZ1_dsx;
Все вместе выглядеть это будет примерно так:
// ...
current_sx = x_start;
length = x_end - x_start;
// расчет u/Z, v/Z, 1/Z, u, v в начале самого первого куска
uZ_a = uZ_start;
vZ_a = vZ_start;
Z1_a = Z1_start; // это 1/Z
u_a = uZ_a / Z1_a;
v_a = vZ_a / Z1_a;
// рисуем куски по 16 пикселов
while (length >= 16) {
uZ_b = uZ_a + 16 * duZ_dsx; // расчет u/Z, v/Z, 1/Z, u, v в конце куска
vZ_b = vZ_a + 16 * dvZ_dsx;
Z1_b = Z1_a + 16 * dZ1_dsx;
u_b = uZ_b / Z1_b;
v_b = vZ_b / Z1_b;
u = u_a; // начинаем текстурирование с начала куска
v = v_a;
du = (u_b - u_a) / 16; // можно сделать >> 4, используя fixedpoint
dv = (v_b - v_a) / 16; // можно сделать >> 4, используя fixedpoint
// рисуем 16 пикселов старым добрым "аффинным" методом
len = 16;
while (len--) {
putpixel(current_sx, current_sy, texture[(int)v][(int)u]);
u += du;
v += dv;
current_sx++;
}
length -= 16;
// конец куска становится началом следующего куска
uZ_a = uZ_b;
vZ_a = vZ_b;
Z1_a = Z1_b;
u_a = u_b;
v_a = v_b;
}
// дорисовываем "хвост" линии, если он непуст
if (length != 0) {
uZ_b = uZ_a + length * duZ_dsx;
vZ_b = vZ_a + length * dvZ_dsx;
Z1_b = Z1_a + length * dZ1_dsx;
u_b = uZ_b / Z1_b;
v_b = vZ_b / Z1_b;
u = u_a; // начинаем текстурирование с начала куска
v = v_a;
du = (u_b - u_a) / length;
dv = (v_b - v_a) / length;
// рисуем остаток пикселов старым добрым "аффинным" методом
while (length--) {
putpixel(current_sx, current_sy, texture[v][u]);
u += du;
v += dv;
current_sx++;
}
}
// ...
Как и в 4.2, пройдемся подобным куском кода по всем строкам грани, не забыв вместо "// ..." вставить интерполяцию всяких там [u/v/1]Z_start, содранную с интерполяции u_start.. и - о чудо, текстурированная с учетом перспективы грань!
Осталось сказать еще пару слов о кое-какой оптимизации.
Во-первых, два деления при расчете u и v в цикле прорисовки можно (и нужно) заменить на одно - посчитать tmp = 1/Z, дальше u = uZ * tmp, v = vZ * tmp.
Во-вторых, немного переставив местами блоки расчета очередной пары точных значений u и v и прорисовки очередного куска линии, можно добиться того, что это самое одно деление, нужное для расчета u и v для *следующего* куска будет находиться сразу перед прорисовкой *текущего* куска. А в этом случае деление может исполняться в сопроцессоре одновременно с отрисовкой куска линии в процессоре. То есть единственная медленная операция будет считаться на полную халяву! Получим перспективно-корректное текстурирование, которое (теоретически) будет работать ненамного медленнее аффинного.
В-третьих, деление на length при дорисовке хвостика длиной от 1 до 15 пикселов можно заменить на умножение на 1/length, заранее посчитав табличку для значений 1/length.
И наконец, мелкие треугольники можно текстурировать аффинным методом, а большие - методом с коррекцией. Размер треугольника можно определять хотя бы по длине самой длинной горизонтальной линии:
x_start = A.sx+(B.sy-A.sy)*(C.sx-A.sx)/(C.sy-A.sy),
x_end = B.sx,
longest_length = x_end - x_start,
4.4. Параболическое
Этот метод основан на приближении u, v квадратичными функциями - параболами, то есть. Для каждой строки строится приближающие u, v квадратичные функции, дальше с их помощью они интерполируются по строке. Для этого нам понадобятся точные значения u, v в трех точках - начале, середине и конце строки. Их считаем точно так же, как в 4.3.
Итак, пусть у нас есть точные значения u в начале, середине и конце строки, то есть на расстоянии 0, length/2 и length пикселов от начала этой строки, обозначим их как ua, ub, и uc соответственно. Мы пытаемся приблизить u квадратичной функцией, то есть полагаем, что
u = A*x*x + B*x + C,
где x - расстояние от текущей точки до начала строки. Тогда, подставив в формулу ua, ub, uc и соответствующие им x, получаем:
ua = C
ub = A*(length*length)/4 + B*length/2 + C
uc = A*(length*length) + B*length + C
Т.о. C = ua, для A и B имеем систему уравнений:
A*(length*length)/4 + B*length/2 = ub - ua
A*(length*length) + B*length = uc - ua
Умножим первое уравнение на четыре, вычтем из него второе:
4*A*(length*length)/4 + 4*B*length/2 -
A*(length*length - B*length = 4*(ub - ua) - (uc - ua),
B*length = 4*(ub - ua) - (uc - ua),
B = (4*(ub - ua) - (uc - ua)) / length.
Умножим первое уравнение на два, вычтем его из второго:
A*(length*length) + B*length -
2*A*(length*length)/4 - 2*B*length/2 = (uc - ua) - 2*(ub - ua),
A*(length*length)/2 = (uc - ua) - 2*(ub - ua),
A = (2*(uc - ua) - 4*(ub - ua)) / (length*length).
Получили формулы для A, B, C. Найдем теперь du и ddu. Для текущей точки x имеем:
du(x) = u(x+1) - u(x),
du = (A*(x+1)*(x+1)+B*(x+1)+C) - (A*x*x+B*x+C) =
= A*(2*x+1) + B
ddu(x) = du(x+1) - du(x),
ddu = (A*(2*(x+1)+1)+B) - (A*(2*x+1)+B) = 2*A
Т.о., начальные значения u, du, ddu будут равны
u = C
du = A + B
ddu = 2*A
А по известным начальным значениям u, du, ddu последовательно вычисляем значения u в любой точке:
u(0), du(0), ddu - известны
u(1) = u(0) + du(0), du(1) = du(0) + ddu
u(2) = u(1) + du(1), du(2) = du(1) + ddu
u(3) = u(2) + du(2), du(3) = du(2) + ddu
...
Для v все делается полностью аналогично.
Таким образом, рисование строки будет выглядеть примерно так:
// ...
// считаем u, v для начала, середины и конца строки
ua = uz_start / z1_start;
va = vz_start / z1_start;
ub = (uz_start + uz_end) / (z1_start + z1_end);
vb = (vz_start + vz_end) / (z1_start + z1_end);
uc = uz_end / z1_end;
vc = vz_end / z1_end;
// считаем начальное du и ddu
pa = 2 * ((uc - ua) - 2 * (ub - ua)) / (length * length);
pb = (4 * (ub - ua) - (uc - ua)) / length;
pc = ua;
u = pc;
du = pa + pb;
ddu = 2 * pa;
// считаем начальное dv и ddv
pa = 2 * ((vc - va) - 2 * (vb - va)) / (length * length);
pb = (4 * (vb - va) - (vc - va)) / length;
pc = v_a;
v = pc;
dv = pa + pb;
ddv = 2 * pa;
// рисуем кусок
while (length--) {
*dest++ = texture[v][u];
u += du;
v += dv;
du += ddu;
dv += ddv;
}
// ...
По сравнению с перспективно-корректным текстурированием имеем более медленный внутренний цикл, но меньшее для длинных строк количество делений. Расчет ua, va и иже с ними можно сделать с помощью трех делений, деления на length и на (length*length) можно заменить на умножение на 1/length и 1/(length*length), беря эти значения из заранее посчитанной таблички. Т.о., на строку приходится три деления независимо от ее длины. Качество более-менее приемлемое; а для коротких строк можно использовать обычную линейную интерполяцию, точно так же, как и в случае с перспективно-корректным текстурированием. Получаем вполне конкурентоспособный метод текстурирования.
5. Освещение
5.1. Модель освещения
Освещенность произвольно взятой точки P, появившуюся из-за источника света, излучающего во все стороны (omnilight) в общем случае будем вычислять по уравнению Фонга:
ambient = Ka,
diffuse = Kd * cos(N, L),
specular = Ks * pow(cos(R, V), Ns),
intensity = ambient + amp * (diffuse + specular).
N
|
L | R V
\ | / /
\ | / /
\ | / /
\ | / /
\|//
-------------P--------------
Здесь использованы следующие обозначения:
Ka - коэффициент фоновой интенсивности (характеристика окружающей среды)
Kd - коэффициент рассеяния (характеристика поверхности)
Ks - коэффициент отражения (характеристика поверхности)
Ns - коэффициент вида отражения (характеристика поверхности)
amp - "мощность" источника света
P - рассматриваемая точка
N - нормаль к поверхности изображаемого объекта в точке P
L - вектор, проведенный из источника света в точку P (луч света)
V - вектор, проведенный из точки P в "точку зрения" камеры
R - отраженный луч света (отражение L относительно N)
ambient - "фоновая" освещенность
diffuse - "рассеянная" освещенность
specular - "отраженная" освещенность
intensity - освещенность (суммарная)
cos(A,B) - косинус угла между векторами A и B
Это уравнение даже поддается относительно простому объяснению - освещенность как бы складывается из фонового уровня освещенности, рассеянного (по всем направлениям - а значит, и по направлению глаза) в этой точке света от источника и отраженного (тоже в глаз) света от источника.
Как обычно, расчеты по этой формуле дадут довольно реалистичный результат, но считать все это для каждой точки грани слишком медленно.
Да, косинус угла между векторами считается как
cos(A,B) = A*B/(|A|*|B|),
где A*B - скалярное произведение векторов, |A| и |B| - их длины. Поэтому имеет смысл все векторы перед использование привести к длине 1, тогда косинус угла между векторами будет равен их скалярному произведению. Небольшая оптимизация.
5.2. Расчет нормали к объекту
Во всех формулах для освещенности у нас так или иначе будет фигурировать вектор N - нормаль к объекту в точке P. Сразу возникает вопрос, а как же этот вектор считать.
Обычно придерживаются такой логики. Модель у нас состоит из плоских граней, но эта сетка плоских граней приближает какой-то искривленный объект. Нормаль к этому искривленному объекту меняется в каждой точке, а для плоских граней она постоянна для всех точек грани, и резко меняется при переходе на другую грань. Поэтому нормаль к объекту обычно приближают следующим образом: считают нормали в вершинах, а нормаль в какой-то точке грани линейно интерполируют между вершинами; то есть линейно интерполируют по грани все три координаты нормали.
Нормаль в вершине рассчитываются как сумма приведенных к длине 1 нормалей ко всем граням, к которым принадлежит эта вершина. То есть. Сначала считаем нормали ко всем граням и приводим их к длине 1. Далее, для каждой вершины надо перебрать все грани, и если очередная грань содержит эту вершину, то к нормали в вершине прибавляется нормаль к этой грани. Первоначально все нормали к вершинам полагаются равными нулю. Для больших моделей этот процесс может быть довольно долгим, но достаточно провести его заранее один раз и сохранить все посчитанные нормали к вершинам.
Для вящей понятности приведу кусок псевдокода:
// ...
for (i = 0; i < numberOfVertics; i++) {
vertexNormal[i].x = 0;
vertexNormal[i].y = 0;
vertexNormal[i].z = 0;
}
for (i = 0; i < numberOfVertics; i++) {
for (j = 0; j < numberOfFaces; j++) {
if (face[j].vertex0 == i ||
face[j].vertex1 == i ||
face[j].vertex2 == i)
{
vertexNormal[i].x += faceNormal[j].x;
vertexNormal[i].y += faceNormal[j].y;
vertexNormal[i].z += faceNormal[j].z;
}
}
}
// ...
5.3. Освещение по Ламберту
Это, видимо, самое сильное упрощение формулы, которое можно придумать. Делаются такие предположения:
- V не сильно зависит от P, т.о. V принимается постоянным для всей грани
- L не сильно зависит от P, т.о. L принимается постоянным для всей грани
- Ks = 0 (то есть грань не отражает свет, а только рассеивает)
- нормаль к объекту N равна нормали к грани n в любой точке грани P
В этом случае формула принимает вид
intensity = ambient + amp * cos(n, L),
где n - нормаль к грани.
Причем в этой формуле, по предположениям, нет величин, зависящих от P. Так что освещенность равна константе для всей грани, то есть все точки грани освещены одинаково. Тогда достаточно посчитать ее лишь один раз на грань. Осталось упомянуть, что освещение по Ламберту обычно принято называть flat shading.
5.4. Освещение по Гуро
Здесь делаются такие предположения:
- Ks = 0 (то есть грань не отражает свет, а только рассеивает)
- освещенность меняется по грани линейно
Ну, а раз освещенность меняется линейно - так же, как и координаты текстуры u, v - можно использовать точно такую же процедуру, как для старого доброго текстурирования, только вместо двух координат текстуры u, v нам надо будет интерполировать одну величину - освещенность. В вершинах она считается по той же самой формуле (с учетом Ks = 0, конечно),
-
intensity = ambient + amp * cos(N, L).
Здесь за N в вершине берется как раз нормаль к объекту в вершине (vertex normal), посчитанная так, как описано в 5.2.
То есть, для написания процедуры освещения по Гуро достаточно написать три строчки расчета освещенности (обозначать ее будем c) в вершинах грани по приведенной здесь формуле и заменить в своей процедуре текстурирования u на c, выкинув при этом v (оно нам не нужно, нас интересует только c). Чтобы добавить освещение по Гуро в текстурирование, надо, соответсвенно, добавить полностью скопированный с расчета u расчет c. Что делать с полученным по цвету найденного через (u, v) в текстуре пиксела и освещенности с, описано в соответсвующем пункте.
Заметим, что освещение по Гуро - это не линейная интерполяция освещенности по изображению грани на экране, как это обычно принято и понимать и делать. Освещенность меняется линейно по трехмерной грани, а вовсе не по ее проекции. То есть, вообще говоря, здесь мы тоже должны учитывать перспективу. Правда, линии одинаковой освещенности все равно будут прямыми, но выглядеть картинка, освещенная по Гуро с учетом перспективы, будет все же покорректнее.
5.5. Освещение по Фонгу
Здесь принято делать, как минимум, такие предположения:
Ks = 0 (то есть грань не отражает свет, а только рассеивает)
Да-да, здесь нет никакой ошибки. Практически все обычно используемые (в demo по меньшей мере) методы т.н. "освещения по Фонгу" НЕ учитывают отраженной компоненты освещенности.
Здесь будет рассказано о самом, наверное, популярном методе освещения по Фонгу, который сводит освещение к текстурированию по определенной текстуре.
Этот метод базируется на таких добавочных предположении:
- L - константа (как бы точечный источник, удаленный бесконечно далеко)
- длина единичной нормали к объекту при интерполяции между вершинами грани НЕ меняется
То есть, утверждается, что если в вершинах нормаль к объекту имеет длину 1, то при интерполяции этой нормали между вершинами по какой-то грани мы будем получать в каждой точке нормаль той же самой длины 1. На самом деле это вовсе даже не так, но для несильно разнящихся углов наклона нормалей приблизительно верно.
Так как Ks = 0, а длина N по предположению равна 1 на всей грани, имеем:
intensity = ambient + amp * (N * L).
Рассмотрим упрощенный случай, когда вектор L = (0,0,1). Общий случай можно без особых вычислительных затрат привести к этому упрощенному, как - будет рассказано чуть позже. Так вот, в этом случае
intensity = ambient + amp * (N.x * L.x + N.y * L.y + N.z * L.z) =
= ambient + amp * (N.x * 0 + N.y * 0 + N.z * 1) =
= ambient + amp * N.z =
= ambient + amp * sqrt(1 - (N.x * N.x + N.y * N.y)).
То есть интенсивность выражается через N.x, N.y, а эти величины меняются линейно. N.x и N.y у нас - числа с плавающей запятой от -1 до 1 (т.к. длина вектора равна 1), интерполировать их - занятие медленно, да корень считать раз в пиксел тоже не хочется. Поэтому вместо интерполяции N.x и N.y обычно интерполируют, например, 128*(N.x+1) и 128*(N.y+1), причем уже в целых числах. Тогда все возможные значения таким образом отмасштабировнных N.x, N.y - это 0, 1, ..., 255. Поэтому можно заранее посчитать табличку значений intensity для каждой пары отмасштабировнных N.x, N.y.
То есть, мы линейно интерполируем 128*(N.x+1) и 128*(N.y+1) (эти значения меняются тоже линейно, раз N.x, N.y меняются линейно) и по ним по таблице определяем интенсивность. Это и есть текстурирование, только в качестве текстуры используется таблица освещенности размером 256x256 (или любым другим), а в качестве координат текстуры u, v для каждой вершины берутся отмасшатбированные координаты нормали в этой вершине.
Таблица, согласно всего вышеупомянутого, считается так:
// ...
for (i = 0; i < 256; i++) {
for (j = 0; j < 256; j++) {
r =
pow((i - 128) / 256.0, 2) + // это N.x*N.x
pow((j - 128) / 256.0, 2); // это N.y*N.y
if (r > 1) r = 1; // длина N меньше 1, поэтому r > 1 быть не может
phongTable[i][j] = amp * sqrt(1 - r);
}
}
// ...
Правда, обычно используют другую - нелинейную - таблицу, видимо, хоть для какой-то компенсации всяких ошибок линеаризации... Результаты выглядят действительно получше. Считается она так:
// ...
for (i = 0; i < 256; i++)
for (j = 0; j < 256; j++)
phongTable[i][j] =
amp * pow(sin(i * PI / 256) * sin(j * PI / 256), 4);
// ...
Для полного комплекта осталось только привести кусочек кода по вычислению координат в этой таблице:
// ...
len = N.x * N.x + N.y * N.y + N.z * N.z;
N.x /= len; // на случай, если длина N не равна 1
N.y /= len;
N.z /= len;
u = (1 + N.x) * 128; // собственно расчет координат
v = (1 + N.y) * 128;
// ...
Теперь вернемся к вопросу о том, как привести случай с произвольным вектором освещения к только что рассмотренному, где L = (0,0,1). Здесь все вроде бы просто. Просто применим к нормалям в вершинах любой поворот, совмещающий наш произвольный вектор света с вектором (0,0,1). Скалярное произведение при этом не изменяется, поэтому так делать можно. Ну, а после этого поворота уже имеем только что расписанный упрощенный случай.
Этот поворот нормалей в каждой вершине не требует практически никаких затрат по следующей причине. Поворот сам по себе, конечно, достаточно медленная процедура и процессорное время отъедает. Но при движении и вращении камеры и объекта мы все равно должны будем соответственно поворачивать нормали. Так вот, эти два поворота можно совместить в один. Если использовать матрицы, все это делается совсем просто - достаточно перемножить (в нужном порядке!) матрицу собственного поворота объекта, матрицу перехода от произвольной камеры к нашей "стандартной" камере и матрицу перехода от произвольного вектора света к "стандартному" вектору света (0,0,1). Т.е. добавится расчет этой матрицы перехода и одно матричное умножение на объект, а это уже мелочь.
6. Разное
6.1. Как совместить текстуру и освещение
Прежде всего, немного теории. Нам понадобится знать то, что конечный цвет пиксела с составляющими (r,g,b) и освещенного цветом (light_r,light_g,light_b) считается как
result_r = r * light_r / max(light_r);
result_g = g * light_g / max(light_g);
result_b = b * light_b / max(light_b);
Здесь max(light_r) - это максимально возможное значение для light_r. Не максимальное по всем тем градациям освещенности, что мы используем, а вообще максимальное для всех возможных градаций. Например, если у нас значения для light_r могут гулять от 0 до 255, а текущий источник света имеет цвет (30, 40,50), и соответсвенно градации освещенности будут равны (k*30,k*40,k*50), где 0 <= k <= 1, то max(light_r) = 255, а не 30. Немного путаное объяснение, но вообще это известная и понятная почти всем вещь.
Метод 1: заранее посчитать таблицу, переводящую пару (цвет, освещенность) в цвет. Можно не менять палитру, а искать для каждой пары (цвет, освещенность) наилучшим образом приближающий ее цвет; а можно (как описано в demo.design FAQ) взять палитру и сгенерировать из нее true color картинку 256x256, в которой пиксел (x, y) нарисован цветом (уже true color!), соответствующим цвету x и освещенности y, а потом перевести ее чем-нибудь типа Image Alchemy в 256 цветов. После чего использовать палитру получившейся картинки, а в качестве таблицы (colorTable) будет сама получившаяся картинка:
outputColor = colorTable[intensity][color].
Метод 2: если нас устроит использовать немного цветов и градаций освещения, то тогда в палитру можно впихнуть все возможные градации всех используемых цветов. Тогда определение нужного индекса в палитре - это одно умножение, а лучше - сдвиг, и одно сложение. Пример: пусть у нас есть 8 цветов и 32 градации освещенности. Палитру заполняем так: 32 градации первого цвета, второго, ..., восьмого. Тогда (для этого примера)
outputColor = (color << 5) + intensity.
Здесь все делается теми же самыми таблицами. Только таблица переводит не цвет в цвет, а компоненту цвета в компоненту цвета. То есть, создаем таблицы redShades[numShades], greenTable[numShades], blueTable[numShades], а потом для каждой компоненты каждого пиксела и нужной градации освещенности по этой таблице определяем выходное значение компоненты:
r = redShades[intensity],
g = greenShades[intensity],
b = blueShades[intensity].
Каждая компонента в этих режимах - это отдельный байт, поэтому никаких проблем не возникает.
Метод 1: тупой, но действенный. Использовать большую таблицу и занести в нее все возможные комбинации цвета и градации освещения. Таблица получится совсем не маленькая, размером 65536*32 = 2 мегабайта. Я написал здесь 32, потому как в этих режимах на компоненту отводится по 5 бит (за исключением 6-битной зеленой компоненты в 16-битном режим), и делать больше градаций освещенности, чем 32, бессмысленно.
Метод 2: делать все так же, как в 24/32-битных режимах. Проблемы возникнут из-за того, что придется с муками выдирать нужные несколько бит компоненты из пиксела. Таблицы для компонент лучше заранее сделать со всеми нужными сдвигами, т.е. значения элементов таблиц должны быть такого вида:
000bbbbb - синий, 8 бит
00000gggggg00000 - зеленый, 16 бит
rrrrr000 - красный, 8 бит
Тогда конечный цвет считается примерно так:
outputColor =
(redTable[(color >> 10) & 0x2F] << 8) +
greenTable[(color >> 5) & 0x1F] +
blueTable[color & 0x1F].
На ассемблере это делается, видимо, побыстрее - и покрасивее. Примерно так:
; ...
mov bx,color
shr bx,10
and bx,02Fh
mov ah,redTable[bx]
mov bx,color
and bx,01Fh
mov al,blueTable[bx]
mov bx,color
shr bx,5 ; можно заменить на
and bx,01Fh ; shr bx,4
shl bx,1 ; and bx,02Eh
or ax,greenTable[bx]
mov outputColor,ax
; ...
Метод 3: рисовать все в 24/32-бита, освещение соответсвенно с текстурой совмещать по пункту 6.1.2, а потом непосредственно при выводе на экран делать преобразование из 24/32-бит в 15/16. Или использовать PTC и предоставить делать нужное преобразование именно ему. PTC - это такая графическаю система для C++, взять ее можно на http://www.gaffer.org/ptc.
6.2. Субпиксельная точность
Субпиксельная точность означает следующее: только те пикселы, центры которых лежат внутри данного многоугольника, должны быть нарисованы. На самом деле, в этом определении можно использовать любую (зафиксированную) точку внутри пиксела, единственным последствием этого будет сдвиг всего экрана меньше, чем на 0.5 пиксела.
Необходимость в субпиксельной точности появляется только из-за того, что дисплеи дискретны, состоят из пикселов. Чем меньше становятся пикселы (то есть чем выше разрешение), тем меньше важна эта точность. Однако разницу можно почувствовать даже на разрешениях порядка 1280x1024, а тем более при обычных для 3D engine 320x200 или 640x480.
Реализовать субпиксельную точность на редкость просто. Представьте себе, что мы рисуем грань. Мы начинаем рисовать ее с какого-то (нецелого) start_y, так как при преобразованиях (например, проецировании) у нас в общем случае будут получаться вовсе не целые числа. Обычно их просто берут и округляют. А для достижения субпиксельной точности достаточно как бы "сдвинуть" грань вверх (или вниз, разницы особо большой нет) до первого целого значения start_y, и соответственно изменить координаты начала и конца строки развертки sx_start, sx_end. Обычно sx_start и sx_end просто линейно интерполируются с помощью градиентов dx_start, dx_end. Так вот в этом случае достаточно просто изменить sx_start, sx_end, а интерполируемые значения будут изменены автоматически:
// ...
sx_start += dx_start * (ceil(start_y) - start_y);
sx_end += dx_end * (ceil(start_y) - start_y);
// ...
То есть мы сдвигаем начальные координаты x для обоих концов соответственно тому, насколько начальное y отстоит от конца своей строки развертки. Ту же самую операцию надо сделать и со всем остальными переменными, которые мы будем интепролировать по ребрам (например u, v, интенсивность для Гуро):
// ...
u_start += du_start * (ceil(start_y) - start_y);
u_end += du_end * (ceil(start_y) - start_y);
// ...
Ну и, разумеется, рисовать начинать надо со строки ceil(start_y).
Вот и все. На скорость работы это не влият вообще, а грани сразу перестают мелко и противно дрожать при перемещении, особенно сильно это заметно при маленьких поворотах (в доли градуса).
6.3. Субтексельная точность
Субтексельная точность сильно смахивает на субпиксельную. Необходимость в ней появляется из-за следующего факта: рисовать строку мы начинаем с какого-то нецелого sx, но с целого пиксела. Из-за этого (sx может "гулять" в пределах одного пиксела почти на единицу!) возникает дрожание текстур примерно на тот самый пиксел. Устраняется это точно таким же сдвигом, как и в субпиксельной точности, то есть перед отрисовкой строки корректируются начальные значения величин, интерполируемых по строке (u, v, освещенность...):
Это уже немного замедляет работу (на величину порядка 5-10%), но улучшение
качества картинки того стоит. Кроме того, субтексельную точность вполне можно
совместить с 2D-отсечением и получить отсечение, которое не съедает никаких
ресурсов и не замедляет работу. Вообще.
6.4. Поворот 3D вектора за шесть умножений
Обычно поворот 3D вектора делают умножением матрицы поворота на этот вектор.
Эта операция требует 9 умножений и 6 сложений. Но с использованием небольшого
precalculation (предварительного расчета) ее можно несколько ускорить.
Пусть нам надо умножить какую-то строку матрицы (a,b,c) на вектор (x,y,z).
Результат должен быть равен
r = a*x+b*y+c*z.
То есть как раз по 3 умножения и 2 сложения на одну строку. Но с другой
стороны,
r = a*x+b*y+c*z = (a*x+b*y+a*b+x*y)+c*z-a*b-x*y =
= (a+y)*(b+x)+c*z-a*b-x*y.
Проще эта формула явно не выглядит. Но дело в том, что x*y - это постоянная
величина, так как x, y - это координаты вершины неповернутого объекта, а они
обычно не меняются. А a*b достаточно посчитать при расчете матрицы поворота,
и это тоже постоянная величина для каждой матрицы. Т.о.
r = (a+y)*(b+x)+c*z-c1-c2.
В результате имеем 2 умножения и 4 сложения на одну строку, то есть те самые
6 умножений и 12 сложений на вектор. Выиграли 3 умножения ценой 6 сложений.
6.5. Билинейная фильтрация текстур
Да-да, это именно тот метод, с помощью которого смазывают текстуры всякие
ускорители типа 3Dfx. Итак, пусть у нас есть какая-то текстура. Текстура у
нас - это 2D картинка, а 2D картинка в свою очередь - набор замеров цвета
через какие-то промежутки. В реальной же жизни цвет не меняется скачком
через каждый, например, миллиметр, а является какой-то непрерывной функцией
от положения, причем меняется довольно плавно. При обычном текстурировании
мы получаем координаты в текстуре, округляем их до ближайшего целого числа
и выбираем нужный цвет из текстуры. То есть мы как бы берем значение цвета
в самой близкой к рисуемой точке сетки замеров цвета, поэтому у нас цвет
резко меняется, оставаясь непрерывным между узлами сетки, поэтому возникает
эффект больших квадратов.
При билинейной фильтрации цвет всего-навсего линейно интерполируется между
узлами сетки замеров. То есть. Пусть в текущей точке у нас получились
координаты текстуры u, v - какие-то нецелые, вообще говоря, числа. Тогда
по целым частям u, v определяется, между какими узлами сетки (если угодно,
между какими пикселами текстуры) находится наша точка, а по дробным - как
близко она находится к каждому из узлов. Вот картинка.
1------------2
| |
a * b
| |
| |
| |
3------------4
Здесь 1, 2, 3, 4 - "окружающие" точку пикселы текстуры (они же узлы сетки
замера цвета). Пусть iu, iv - целые части координат текстуры точки u, v; fu,
fv - дробные части. Тогда 1, 2, 3, 4 имеют координаты в текстуре (iu,iv),
(iu+1,iv), (iu,iv+1), (iu+1,iv+1). Проинтерполируем какую-то компоненту
цвета (R, G или B) по прямым 1-3 и 2-4:
a.c = 1.c + (3.c - 1.c) * fv;
b.c = 2.c + (4.c - 2.c) * fv;
то есть
a.c = c[iu][iv] + (c[iu][iv+1] - c[iu][iv]) * fv;
b.c = c[iu+1][iv] + (c[iu+1][iv+1] - c[iu+1][iv]) * fv;
Теперь проинтерполируем цвет по прямой ab в нашей точке:
c = a.c + b.c * fu;
Проинтерполировав по этой формуле каждую компоненту цвета, получим наконец
готовый результат - цвет точки, но уже с учетом билинейной фильтрации.
Здесь у нас получилось по три умножения на компоненту. То есть в сумме девять
умножений на пиксел. Можно, конечно, честно считать по этим формулам, делая
девять умножений для каждого пиксела. Но можно заменить все умножения на
выборки по таблице. u, v обычно - это fixedpoint; fu, fv - тоже (кстати, в
случае с fixedpoint целые и дробные части вычисляются ровно одним and'ом).
Пусть мы используем 24-битный цвет и 16:16 fixedpoint; тогда одна компонента
цвета занимает 8 бит, а дробную часть можно одним сдвигом перевести в 24:8
fixedpoint. Получаем 256 возможных значений для любой компоненты цвета и 256
возможных значений для дробной части, то есть - табличка 256x256. Если цвет
15/16-битный, или используется более грубое (скажем, до пяти бит) округление
дробной части, то табличка становится еще меньше. Памяти, конечно, не жалко,
но кэш-память пока не резиновая, так что чем меньше lookup-таблица, тем оно
лучше для скорости. Вот и все.
Осталось только упомянуть, что лучше занести в табличку не байты, а слова,
для данного примера это будет 8:8 fixedpoint, и складывать все результаты
тоже как слова, а потом сдвигом переводить обратно в целые числа. Иначе
(особенно в случае 15/16-битных режимов) будет заметен небольшой шум на
текстуре, появляющийся из-за ошибок округления.
6.6. Алгоритм "бегущих кубиков" для полигонизации изоповерхностей
Сооветствующий английский термин - marching cubes algorithm; приведен здесь,
так сказать, just to avoid any confusion. Алгоритм предназначен для быстрого
построения полигональной модели изоповерхности трехмерного скалярного поля,
заданного значениями на равномерной сетке. Выражаясь менее академичным и,
надеюсь, более понятным языком, этот алгоритм нужен для быстрого построения
набора треугольных граней, достаточно хорошо приближающего изоповерхность,
то есть такую поверхность, на которой определенная функция равна какой-то
константе - изоуровню. Например, сфера радиуса 1 с центром в нуле - это
изоповерхность функции f=1/(x*x+y*y+z*z) для изоуровня 1. Т.н. "3D капли",
столь любимые ныне, тоже являются изоповерхностью, правда, немного более
сложной функции (да, функция задана в трехмерном пространстве):
n
f = sum c[i]/((x[i]-x)*(x[i]-x)+(y[i]-y)*(y[i]-y)+(z[i]-z)*(z[i]-z)),
i=1
где
n - количество источников поля (капелек);
c[i] - интенсивность каждого источника (радиус капельки);
x[i], y[i], z[i] - координаты каждого источника (центра капельки).
Работает алгоритм следующим образом.
Рассматриваем какой-то параллелипипед, внутри которого заведомо находится
изоповерхность (или та ее часть, которую мы хотим нарисовать). Разбиваем
его сеткой, ну, как бы разрезаем его на несколько меньших параллелипипедов,
и считаем значения функции (поля) в узлах сетки, то есть вершинах этих самых
маленьких параллелипипедов. Для большей ясности можно заменить длинное слово
"параллелипипед" на слово "куб" и представить себе кубик Рубика. Теперь
проходим по всем кубикам (дальше я буду везде использовать слово "кубик",
надоело "параллелипипед" писать). Смотрим на значения функции в вершинах
кубика. Если все они больше (или меньше) изоуровня - значит, кубик находится
целиком над (или под) изоповерхностью, внутри этого кубика поверхности нет и
мы его просто отбрасываем. Если же часть больше, а часть меньше, то некоторые
ребра кубика пересекаются с изоповерхностью. Линейной интерполяцией приближаем
эти точки пересечения и в зависимости от того, какие вершины находятся над
изоповерхностью, а какие под, генерируем несколько треугольных граней. Все
вершины этих граней - это как раз точки пересечения поверхности с ребрами.
Как генерировать грани и пересечения каких ребер с поверхностью считать,
смотрим по таблице, это зависит лишь от того, какие вершины находятся над
поверхностью, а какие - под. Вершин восемь, состояний два - над и под. Это
дает нам 256 возможных расположений, так что таблица не такая уж и большая.
Индекс в таблице тоже генерируется совсем просто: если вершина находится
над поверхностью, устанавливаем соответствующий этой вершине бит индекса,
иначе - сбрасываем.
Вот простенький пример.
4--------3
/ | / |
1--------2 |
| | * |
| 8----|---7
| / | *
5---*----6
Пусть точка 6 находится под поверхностью, остальные - над. Тогда поверхность
проходит через точки (*), приближаем их линейной интерполяцией и проводим
через них треугольную грань. Какие ребра пересекаются с поверхностью, какие
точки пересечния и как соединять - узнаем из таблиц по индексу; в данном
случае индекс равен 11011111b=0DFh (установлены все биты, кроме 6го).
Алгоритм по шагам:
* посчитать значения функции в узлах сетки (вершинах кубиков)
* для каждого кубика:
- посчитать индекс этого кубика в таблице. Каждый бит соответствует
какой-то из вершин, если в ней значение функции больше изоуровня,
то надо установить этот бит, иначе - сбросить
- взять из таблицы по индексу кубика слово, каждый бит в котором
определяет, пересекает ли соответствующее ему (биту) ребро кубика
изоповерхность, или нет. Если да, то посчитать (простая линейная
интерполяция между соответствующими ребру вершинами) координаты
точки пересечения
- взять из таблицы по индексу кубика число генерируемых граней и
собственно тройки номеров ребер, (уже посчитанные) пересечения
которых с изоповерхностью и будут вершинами нужных граней
Неоднократно упоминавшиеся магические таблицы приведены в примерчике, там
же приведен кусок кода, который довольно легко переделть под свой engine.
Пример взят с http://www.mhri.edu.au/~pdb/modelling/polygonise, сам по себе
он компилироваться НЕ будет, но все там написано правильно и никаких ошибок
в таблицах нет (так что ищите ошибку в своем коде).