Comments 38
Вот эта картинка — это одна заглавная картинка минус другая.Это немного нечестно. Из-за неравномерности sRGB значения в темном диапазоне слабо отличаются на глаз. Если инвертировать картинку (и протереть монитор от пыли), различия станут намного заметнее.
![](https://habrastorage.org/files/bd9/5a3/f6f/bd95a3f6f0874587ac5c6058e50a24dc.png)
Немного про оптимизации и 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
Мне интересен более простой способ вычисления барицентрических координат в клип-пространстве из координат 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);
Матрица А — это матрица, строками в которой являются координаты вершин треугольника, с заменой третьей координаты на 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;
}
Отправить деление 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 оптимизации, а не математические. Поэтому это скорее к вашему соавтору по части оптимизаций.
На картинке отмечен лишний пиксель при растеризации верхнего треугольника. Его нужно отдельно учитывать для целочисленного алгоритма. Если цикл на float — этой проблемы не будет. Кроме того цикл на float позволит заменить постоянное деление в for по X на умножение. Ну а на современных процессорах операции с float уже такие же быстрые как и целочисленные.
Ни в плане изящества кода, ни качества получаемого результата (он дает артефакты тонирования).
А по теме вот занимательное чтение, которое пробежало в конференции.
У меня два куба, текстуры интерполируются нормально, но посмотрите на вертикальную линию пересечения кубов — в случае, когда я использую bc_clip, она кривая (кубы вверху), а когда использую bc_screen — она прямая (кубы внизу).
![](https://habrastorage.org/files/0b2/ba0/222/0b2ba0222e4c4921bdf702c39b982b88.png)
Опять нужен сеанс дебага по фотографии.
![Без учёта перспективы Без учёта перспективы](https://habrastorage.org/getpro/habr/upload_files/57a/1e8/416/57a1e8416525b27ca298baf0af549d0c.png)
![На основе доперспективных барицентрических кооринат На основе доперспективных барицентрических кооринат](https://habrastorage.org/getpro/habr/upload_files/b22/8f2/49a/b228f249a8e39a66d8fa45deb680f309.png)
Явно стало лучше, но как-то не до конца. Как же так?
Опять уже сдался было, сел ужинать и, внезапно, понял, где, вероятно, корень проблемы.
Преобрзование растяжения у меня выполнялось по прежнему через x * minMid + xMid, а не матрицами, причём координату z я после перспективного искажения не преобразовывал вовсе. Стоило домножить её на 2 (r у меня был выставлен в -0.5) и хоть я не уверен пока (спешу опубликовать находку), почему именно 2 (эксперименты показывают, что это явно не -1/r), но суть ясна - надо сначала реализовать честную камеру на матрицах.
P. S. Кстати, подозрительно справа проступает синяя полоска, наверное из-за ещё каких-то ошибок округления или граничных случаев.
![](https://habrastorage.org/getpro/habr/upload_files/530/4c4/7c3/5304c47c372c7ad63d74dac3c324998c.png)
Да и верхняя горизонтальная граница стала зубчатая.
В упомянутых в статье исходниках вижу какое-то противоречие приводимым математическим выводам.
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)? Они как-то хитро сократились из-за последующего деления на сумму тех же компонент?
Наверное дело в
rA.z+1, rB.z+1, rC.z+1 нам известны, это координаты треугольника, переданные в растеризатор
Я то подумал, что "это" имеется в виду A.z, B.z, C.z, а имеются в виду целиком выражения.
Только теперь непонятно, каким образом это могут быть координаты перданного треугольника, если они выглядели вот так
![](https://habrastorage.org/getpro/habr/upload_files/98f/44e/09d/98f44e09d1ca137f3f8f69cffd8941e0.png)
Кажется у меня в комментарии выше всё отрисовалось правильно из-за удачной комбинации двух (или даже трёх) ошибок. Осталось понять, каких именно.
Но начинает казаться, что слишком стар я стал для всех этих "легко видеть, что...", и дальше кривого текстурирования мне не осилить. Где-то упоминалось, что весь курс часов на двадцать, я, кажется, на половину курса уже столько потратил.
Убил ещё 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 как будто бы помогало (когда я проверил на модели головы, проявились другие артефакты).
Но ведь pts[0][3] - это не A.z, это A.w. Который, в свою очередь, равен (после умножения на матрицу проекции) r*A.z+1? Передаётся-то четыре координаты на каждую точку, т.к. см. однородные координаты...
Краткий курс компьютерной графики: пишем упрощённый OpenGL своими руками, статья 4в из 6