Краткий курс компьютерной графики: пишем упрощённый OpenGL своими руками, статья 4в из 6

[homm]

Вот эта картинка — это одна заглавная картинка минус другая.

Это немного нечестно. Из-за неравномерности sRGB значения в темном диапазоне слабо отличаются на глаз. Если инвертировать картинку (и протереть монитор от пыли), различия станут намного заметнее.

[haqreu]

Да, голова не исчезла полностью, оставила немного перхоти. Но это не повод копья ломать :)

[splav_asv]

Тут скорее дело в том, что у головы близкие цвета в соседних точках текстуры.
В итоге после изменения контрастности получилось это:

[haqreu]

Невооружённым глазом это не видно. В отличие от проблем на полу.

[haqreu]

 

[MrShoor]

Хех, спасибо за статью, и за упоминание меня.
Немного про оптимизации и triangle. bc_screen не нужен. Он у вас сейчас используется только чтобы определить принадлежность точки треугольнику. Если избавится от bc_screen — то и вычисления bc_clip будут оптимальнее. Я когда прошлый раз смотрел ваш код, не обратил внимание на то, как считается barycentric.
Вообще красивее это сделать через сканлайны, тогда вообще не нужно будет проверять лежит точка в треугольнике или нет, и мы гарантированно будем рисовать только то, что лежит в треугольнике. Сканлайны не так страшны в реализации, как кажется на первый взгляд.
Любой треугольник можно разбить максимум на 2 горизонтальной прямой, проходящей через одну вершину:

Это дает нам возможность создать 2 цикла по Y, от верхней вершины до разделительной линии. И от разделительной линии до нижней вершины. Поэтому первый этап — сортировка 3-х вершины по Y, чтобы определиться с границами первого цикла.
Второй этап — сортировка угловых коэффициентов линий треугольника, чтобы идти от левого ребра к правому.
Буквально несколько месяцев назад набросал коллеге код такого целочисленного сканлайна в качестве примера. Сейчас нашел этот код в логах скайпа:

Собственно код (паскаль)

procedure Exchange(var p0, p1: TPoint); inline;
var tmp: TPoint;
begin
  tmp := p1;
  p1 := p0;
  p0 := tmp;
end;

procedure DrawTrinangle(canvas: TCanvas; p0,p1,p2: TPoint);
var
  j, i: Integer;
  dirLeft, dirRight: TPoint;
begin
  if p0.Y > p1.Y then Exchange(p0, p1); //сортируем вершины по Y
  if p1.Y > p2.Y then Exchange(p1, p2);
  if p0.Y > p1.Y then Exchange(p0, p1);

  dirLeft.X := p2.X - p0.X;
  dirLeft.Y := p2.Y - p0.Y;

  dirRight.X := p1.X - p0.X;
  dirRight.Y := p1.Y - p0.Y;

  if dirRight.X*dirLeft.Y < dirLeft.X*dirRight.Y then Exchange(dirRight, dirLeft); //сортируем угловые коэффициенты верхнего треугольника

  for j := p0.Y to p1.Y - 1 do
    for i := (j-p0.Y)*dirLeft.X div dirLeft.Y + p0.X to (j-p0.Y)*dirRight.X div dirRight.Y + p0.X do
      canvas.Pixels[i, j] := clBlue;

  dirLeft.X := p0.X - p2.X;
  dirLeft.Y := p0.Y - p2.Y;

  dirRight.X := p1.X - p2.X;
  dirRight.Y := p1.Y - p2.Y;

  if dirRight.X*dirLeft.Y > dirLeft.X*dirRight.Y then Exchange(dirRight, dirLeft); //сортируем угловые коэффициенты нижнего треугольника

  for j := p1.Y to p2.Y - 1 do
    for i := (j-p2.Y)*dirLeft.X div dirLeft.Y + p2.X to (j-p2.Y)*dirRight.X div dirRight.Y + p2.X do
      canvas.Pixels[i, j] := clBlue;
end;

Там есть мелкие «недочеты» (типа одного лишнего пикселя, из-за целочисленности). Но в целом как видно совсем небольшой и не сложный код, дающий ощутимый прирост к производительности, особенно на больших тонких и длинных треугольниках.

upd. Но это все оптимизации, я думаю с этими вещами лучше к gbg

[haqreu]

Ну собственно сканлайны у меня и были закодированы во второй статье. Только после этого всё равно придётся делать вычисления нормальных бар. координат, если мы хотим нормальной проекции.

Мне интересен более простой способ вычисления барицентрических координат в клип-пространстве из координат 2д пикселя. Вы упоминаете, что bc_clip можно проще вычислять. Я не вижу куда ещё упрощать, вы что можете предложить?

            Vec3f bc_screen  = barycentric(pts2[0], pts2[1], pts2[2], P);
            Vec3f bc_clip    = Vec3f(bc_screen.x/pts[0][3], bc_screen.y/pts[1][3], bc_screen.z/pts[2][3]);
            bc_clip = bc_clip/(bc_clip.x+bc_clip.y+bc_clip.z);

[gbg]

Барицентрические координаты можно получать, умножая матрицу (A^-1)^T на вектор (x,y,1)
Матрица А — это матрица, строками в которой являются координаты вершин треугольника, с заменой третьей координаты на 1.
Матрицу нужно обратить один раз.

Исходник примера

#include <iostream>
 
#include <iomanip>
#include <vector>
 
//#include <vec-simple.h>
 
#ifndef VECSIMPLE_H
#define VECSIMPLE_H
#ifndef VEC_H
#define VEC_H
#include <ostream>
 
#define ASSERT(n) (n);
using namespace std;
 
template <size_t Dim,typename Number> class vec
{
 
public:
typedef Number NumberT;
static const size_t DimN=Dim;
Number items[Dim];
 
vec()
{
 
}
 
vec(const vec<Dim,Number >& src)
{
for(size_t i=Dim;i--;)
{
items[i]=src[i];
}
}
 
vec(const Number* src)
{
for(size_t i=Dim;i--;)
{
items[i]=src[i];
}
}
 
size_t maxPos() const
{
size_t ret=0;
for(size_t i=Dim;--i;)
{
if(items[i]>items[ret])
{
ret=i;
}
}
return(ret);
}
 
ostream& print(ostream& out) const
{
out<<"{ ";
for(size_t i=0;i<Dim;i++)
{
out<<setw(6)<<items[i]<<" ";
}
out<<"} ";
return(out);
}
static vec<Dim,Number > fill(const Number& val=0)
{
vec<Dim, Number> ret;
for(size_t i=Dim;i--;)
{
ret[i]=val;
}
return(ret);
}
 
bool operator!=(const vec<Dim,Number>&v)
{
for(size_t i=Dim;i--;)
{
if(v[i]!=items[i])
{
return(true);
}
}
return(false);
}
 
Number& operator [](size_t index)
{
return(items[index]);
}
const Number& operator [](size_t index) const
{
return(items[index]);
}
};
 
template<size_t Dim,typename Number>vec<Dim,Number > operator+(vec<Dim,Number > lhs, const vec<Dim,Number >& rhs)
{
for(size_t i=Dim;i--;)
{
lhs[i]+=rhs[i];
}
return(lhs);
 
}
 
template<size_t Dim,typename Number>vec<Dim,Number > operator-(vec<Dim,Number > lhs, const vec<Dim,Number >& rhs)
{
for(size_t i=Dim;i--;)
{
lhs[i]-=rhs[i];
}
return(lhs);
}
 
template<size_t Dim,typename Number>vec<Dim,Number > operator*(vec<Dim,Number > lhs, const Number& rhs)
{
for(size_t i=Dim;i--;)
{
lhs[i]*=rhs;
}
return(lhs);
}
 
 
 
template<size_t Dim,typename Number>vec<Dim,Number > operator/(vec<Dim,Number > lhs, const Number& rhs)
{
for(size_t i=Dim;i--;)
{
lhs[i]/=rhs;
}
return(lhs);
}
 
 
template<size_t Dim,typename Number> Number operator*(const vec<Dim,Number >&lhs, const vec<Dim,Number >& rhs)
{
Number ret=0;
for(size_t i=Dim;i--;)
{
ret+=lhs[i]*rhs[i];
}
return(ret);
}
 
 
template<size_t len,size_t Dim, typename Number> vec<len,Number > proj(const vec<Dim,Number> &v,size_t start=0)
{
return(vec<len,Number >(&v.items[start]));
}
 
template<size_t len,size_t Dim, typename Number> vec<len,Number > dive(const vec<Dim,Number> &v,size_t start=0)
{
vec<len,Number> ret(&v.items[start]);
for(size_t i=Dim+start;i<len;i++)
{
ret[i]=1;
}
return(ret);
}
 
template<size_t Dim,typename Number> ostream& operator<<(ostream& os,const vec<Dim,Number >& v)
{
return(v.print(os));
}
 
template<size_t DimRows,size_t DimCols,typename Number> class mat;
 
template<size_t DimCols,size_t DimRows,typename Number> struct dt
{
static Number det(const mat<DimRows,DimCols,Number>& src)
{
Number ret=0;
for(size_t i=DimCols;i--;)
{
ret+=src[0][i]*src.algAdd(0,i);
}
return(ret);
}
};
 
template<typename Number> struct dt<1,1,Number>
{
static Number det(const mat<1,1,Number>& src)
{
return(src[0][0]);
}
};
 
template<size_t DimRows,size_t DimCols,typename Number> class mat
{
vec<DimCols,Number> rows[DimRows];
 
 
public:
typedef Number NumberT;
static size_t shift(size_t in,const size_t& val)
{
return(in<val ? in : ++in);
}
 
mat()
{
}
mat(const mat<DimRows,DimCols,Number >& src)
{
for(size_t i=DimCols;i--;)
{
for(size_t j=DimRows;j--;)
{
const Number t=src[i][j];
rows[i][j]=t;
}
}
}
 
ostream& print(ostream& out) const
{
for(size_t i=0;i<DimRows;i++)
{
out<<rows[i]<<"\n";
}
return(out);
}
vec<DimCols,Number > minimums()
{
vec<DimCols,Number > ret=rows[0];
for(size_t i=DimRows;--i;)
{
for(size_t j=DimCols;j--;)
{
ret[j]=min(ret[j],rows[i][j]);
}
}
return(ret);
}
 
vec<DimCols,Number > maximums()
{
vec<DimCols,Number > ret=rows[0];
for(size_t i=DimRows;--i;)
{
for(size_t j=DimCols;j--;)
{
ret[j]=max(ret[j],rows[i][j]);
}
}
return(ret);
}
vec<DimCols,Number>& operator[] (size_t index)
{
return(rows[index]);
}
 
const vec<DimCols,Number>& operator[] (size_t index) const
{
return(rows[index]);
}
static mat<DimCols,DimRows,Number> ones()
{
mat<DimCols,DimRows,Number> ret;
for(size_t i=DimRows;i--;)
{
for(size_t j=DimCols;j--;)
{
ret[i][j]=(i==j);
}
 
}
return(ret);
}
 
Number det() const
{
return(dt<DimCols,DimRows,Number>::det(*this));
}
 
mat<DimRows-1,DimCols-1,Number> minor(size_t row,size_t col) const
{
mat<DimRows-1,DimCols-1,Number> ret;
for(size_t i=DimRows-1;i--;)
{
for(size_t j=DimCols-1;j--;)
{
ret[i][j]=rows[ret.shift(i,row)][ret.shift(j,col)];
}
}
return(ret);
}
 
 
Number algAdd(size_t row,size_t col) const
{
return(minor(row,col).det()*( (row+col)%2 ? -1 : 1));
}
 
mat<DimRows,DimCols,Number> Adjacent()const
{
mat<DimRows,DimCols,Number> ret;
for(size_t i=DimRows;i--;)
{
for(size_t j=DimCols;j--;)
{
ret[i][j]=algAdd(i,j);
}
}
return(ret);
}
 
mat<DimRows,DimCols,Number> invertT()const
{
mat<DimRows,DimCols,Number> ret=Adjacent();
return(ret/(ret[0]*rows[0]));
}
 
void setCol(const Number& val,size_t col)
{
for(size_t i=DimRows;i--;)
{
rows[i][col]=val;
}
}
 
};
 
template<size_t Dim,typename Number>vec<Dim,Number > operator*(const mat<Dim,Dim,Number >& lhs, const vec<Dim,Number>& rhs)
{
vec<Dim,Number> ret;
for(size_t i=Dim;i--;)
{
ret[i]=lhs[i]*rhs;
}
return(ret);
}
 
template<size_t DimCols,size_t DimRows,typename Number>mat<DimCols,DimRows,Number > operator/(mat<DimCols,DimRows,Number> lhs, const Number& rhs)
{
for(size_t i=DimRows;i--;)
{
lhs[i]=lhs[i]/rhs;
}
return(lhs);
}
 
template<size_t DimRows,size_t DimCols,typename Number> ostream& operator<<(ostream& os,const mat<DimRows,DimCols,Number>& v)
{
return(v.print(os));
}
 
typedef mat<2,2,float > tempMat;
 
#endif // VEC_H
 
 
#endif // VECSIMPLE_H
 
 
 
using namespace std;
 
typedef vec<3,float> pixel;
 
class screen
{
vector<pixel > buffer;
size_t width;
size_t height;
public:
screen(size_t width,size_t height):width(width),height(height)
{
buffer.assign(width*height,pixel());
}
template<typename t> void putpixel(const t& location,const pixel& px)
{
 
buffer[location[0]+location[1]*width]=px;
}
 
ostream& print(ostream& out) const
{
for(size_t i=0;i<height;i++)
{
out<<"|";
for(size_t j=0;j<width;j++)
{
const size_t mp=buffer[i+j*width].maxPos();
static const char syms[]={' ','.','\'','-','+','*'};
static const char syms2[]={'R','G','B'};
if(buffer[i+j*width][mp] > static_cast<pixel::NumberT>(sizeof(syms)-1)/static_cast<pixel::NumberT>(sizeof(syms)))
{
cout<<syms2[mp];
}
else
{
const size_t idx=buffer[i+j*width][mp]*float(sizeof(syms)-1);
out<<syms[idx];
}
}
out<<"|\n";
}
return(out);
}
};
 
typedef vec<3,float > v3i;
typedef vec<2,float > v2i;
typedef mat<3,3,float> m3i;
 
 
template<typename V>class allPointsOfSquare
{
vec<V::DimN,typename V::NumberT> topLeft;
vec<V::DimN,typename V::NumberT> bottomRight;
vec<V::DimN,typename V::NumberT> pos;
public:
allPointsOfSquare(const V& topLeft,const V& bottomRight):topLeft(topLeft),bottomRight(bottomRight),pos(topLeft)
{
}
const v2i& operator *() const
{
return(pos);
}
bool next()
{
const bool ret=(pos!=bottomRight);
for(size_t i=V::DimN;i--;)
{
pos[i]++;
if(pos[i]>bottomRight[i])
{
pos[i]=0;
}
else
{
break;
}
}
return(ret);
}
};
 
class IShader
{
public:
virtual ~IShader()
{
 
}
 
virtual pixel shade(const v3i& a)=0;
};
 
class baryShader:public IShader
{
public:
virtual ~baryShader()
{
 
}
 
virtual pixel shade(const v3i&a)
{
return(a);
}
};
 
 
void fillTria(mat<3,3,float > coord,IShader* shader,screen& scr)
{
//находим углы прямоугольника, в котором лежит треугольник.
//это соответственно минимумы и максимумы сторон, его содержащих
 
v3i topLeft=coord.minimums();
v3i bottomRight=coord.maximums();
 
coord.setCol(1,2);
 
mat<3,3,float > bcm=coord.invertT();
 
mat<3,2,float > directions;
 
for(size_t i=3;i--;)
{
directions[i]=proj<2>(coord[ (i+1) % 3 ])-proj<2>(coord[i]);
}
 
for(allPointsOfSquare<v2i> sweep=allPointsOfSquare<v2i>(
proj<2>(topLeft),
proj<2>(bottomRight))
;
sweep.next()
;
)
{
size_t i=0;
for(;i<3;i++)
{
v2i curDirection=(*sweep)-proj<2>(coord[i]);
 
tempMat crss;
crss[0]=proj<2>(curDirection);
crss[1]=proj<2>(directions[i]);
 
if(crss.det()>0)
{
break;
}
}
if(i==3)
{
v3i a=bcm*dive<3>(*sweep);
 
scr.putpixel(*sweep,shader->shade(a));
 
}
}
}
 
 
 
 
ostream& operator<<(ostream& os,const screen& v)
{
return(v.print(os));
}
 
 
 
 
int main()
{
screen scr(80,80);
// cout<<scr;
mat<3,3,float > t;
t[0][0]=5;
t[0][1]=7;
t[0][2]=3;
 
 
t[1][0]=70;
t[1][1]=13;
t[1][2]=13;
 
 
 
t[2][0]=20;
t[2][1]=50;
t[2][2]=21;
cerr<<"started";
 
 
baryShader bs;
fillTria(t,&bs,scr);
cerr<<"filled";
 
cout<<scr;
 
return 0;
}
 

он же на ideone

[MrShoor]

Внутри barycentric выкинуть проверку: if (std::abs(u[2])>1e-2)
Отправить деление bc_screen.x/pts[0][3], bc_screen.y/pts[1][3], bc_screen.z/pts[2][3] внутрь barycentric, что позволит заменить деление на умножение. Так же отправить туда bc_clip = bc_clip/(bc_clip.x+bc_clip.y+bc_clip.z);
Ну то есть это уже чисто CPU оптимизации, а не математические. Поэтому это скорее к вашему соавтору по части оптимизаций.

[haqreu]

А, ну чуть укоротить код — это немного не то, что мне надо. При таких раскладах я оставлю bc_screen как есть, т.к. у него есть чёткая семантика, иначе становится совсем непонятно, что и когда мы считаем.

[gbg]

Этот алгоритм сразу вызывает вопрос — не наделает ли он дырок при построении смежных треугольников.

[MrShoor]

Не наделает:

На картинке отмечен лишний пиксель при растеризации верхнего треугольника. Его нужно отдельно учитывать для целочисленного алгоритма. Если цикл на float — этой проблемы не будет. Кроме того цикл на float позволит заменить постоянное деление в for по X на умножение. Ну а на современных процессорах операции с float уже такие же быстрые как и целочисленные.

[gbg]

Спасибо, вот и еще один кандидат на бенчмарк нарисовался.

[MrShoor]

Только обязательно переделать на float. А то лишние пиксели — это грубо говоря ошибка, да и деление на каждую строку сканлайна из-за целочисленности избыточно.

[dom1n1k]

Скайнлайны — это традиционный и самый известный метод, но он нифига не подходит под определение красивого.
Ни в плане изящества кода, ни качества получаемого результата (он дает артефакты тонирования).

[MrShoor]

Про артефакты тонирования — полный бред. Ознакомьтесь с методом, прежде чем так говорить.

[haqreu]

Ох. Давайте воздержимся от выражений типа «полный бред», по крайней мере, в публичной дискуссии, хорошо?

А по теме вот занимательное чтение, которое пробежало в конференции.

[valeriyk]

Нужно ли использовать коррекцию перспективы для Z-буфера?
У меня два куба, текстуры интерполируются нормально, но посмотрите на вертикальную линию пересечения кубов — в случае, когда я использую bc_clip, она кривая (кубы вверху), а когда использую bc_screen — она прямая (кубы внизу).

[MrShoor]

Да. Нужно. Вот тут объяснение на пальцах.

[haqreu]

Вы имеете в виду со сканирующей линией? Точно так же. Просто вместо билинейной интерполяции посчитайте барицентрические координаты текущего пикселя и сделайте поправку, описанную в этой статье.

[haqreu]

Вполне можно, нужно просто скорректировать координаты концов сканирующей линии.

[valeriyk]

Кстати, было бы неплохо получить адд-он (здесь или на гитхабе) про то, как негру добавить глаза!

[haqreu]

Эээ. Модель глаз лежит рядом с негром, достаточно её отрендерить вместе с головой, нет?

[valeriyk]

ааа, так это модель двух глаз сразу?

[haqreu]

Да.

[ferocactus]

Опять нужен сеанс дебага по фотографии.

Без учёта перспективы

На основе доперспективных барицентрических кооринат

Явно стало лучше, но как-то не до конца. Как же так?

[ferocactus]

Опять уже сдался было, сел ужинать и, внезапно, понял, где, вероятно, корень проблемы.

Преобрзование растяжения у меня выполнялось по прежнему через x * minMid + xMid, а не матрицами, причём координату z я после перспективного искажения не преобразовывал вовсе. Стоило домножить её на 2 (r у меня был выставлен в -0.5) и хоть я не уверен пока (спешу опубликовать находку), почему именно 2 (эксперименты показывают, что это явно не -1/r), но суть ясна - надо сначала реализовать честную камеру на матрицах.

P. S. Кстати, подозрительно справа проступает синяя полоска, наверное из-за ещё каких-то ошибок округления или граничных случаев.

Да и верхняя горизонтальная граница стала зубчатая.

[ferocactus]

В упомянутых в статье исходниках вижу какое-то противоречие приводимым математическим выводам.

vec3 bc_screen = barycentric(pts2, vec2(x, y));
vec3 bc_clip   = vec3(bc_screen.x/pts[0][3], bc_screen.y/pts[1][3], bc_screen.z/pts[2][3];););
bc_clip = bc_clip/(bc_clip.x+bc_clip.y+bc_clip.z);

Куда же пропали "r *" и "+1" из знаменателя? Разве там не должно быть:
bc_screen.x/(r * (pts[0][3] + 1)? Они как-то хитро сократились из-за последующего деления на сумму тех же компонент?

[ferocactus]

Наверное дело в

 rA.z+1, rB.z+1, rC.z+1 нам известны, это координаты треугольника, переданные в растеризатор

Я то подумал, что "это" имеется в виду A.z, B.z, C.z, а имеются в виду целиком выражения.
Только теперь непонятно, каким образом это могут быть координаты перданного треугольника, если они выглядели вот так

Кажется у меня в комментарии выше всё отрисовалось правильно из-за удачной комбинации двух (или даже трёх) ошибок. Осталось понять, каких именно.

Но начинает казаться, что слишком стар я стал для всех этих "легко видеть, что...", и дальше кривого текстурирования мне не осилить. Где-то упоминалось, что весь курс часов на двадцать, я, кажется, на половину курса уже столько потратил.

[ferocactus]

Убил ещё 2 часа, но это полный тупик.

Я заметил, что под координатам переданными в растеризатор имеются в виду исходные кординаты clip пространства (после цепочки преобразований, включая перспективное), но до проецирования из 4d в 3d. Так что это не то, что на картинке в предыдущем комментарии, но и никак не rA.z+1, rB.z+1, rC.z+1. Видимо, в коде автора, либо ошибка, которая чудом не влияет на результат, либо им таки были проделаны некие преобразования, показывающие эквивалентность применения rA.z+1, rB.z+1, rC.z+1 и A.z, B.z, C.z.

В то же время, не ясно, почему растяжение и смещение по x, y не с помощью матриц, выполненное после проецирования, ломает алгоритм, и почему домножение на 2 как будто бы помогало (когда я проверил на модели головы, проявились другие артефакты).

[ferocactus]

Потратив ещё часов 5, реализовал честную цепочку преобразований на матрицах, и всё заработало.

Это лишь подчёркивает вопрос - куда же у автора пропали "r *" и "+1" из знаменателей при вычислении барицентрических координат?

[haqreu]

Но ведь pts[0][3] - это не A.z, это A.w. Который, в свою очередь, равен (после умножения на матрицу проекции) r*A.z+1? Передаётся-то четыре координаты на каждую точку, т.к. см. однородные координаты...

[ferocactus]

Спасибо. Хитро, как я и ожидал. Значит, и мне надо будет сэкономить на его вычислении.