算法来自Microfacet BSDFs Generated from NDFs and Explicit Microgeometry
具体算法
自己实现代码时候的思路与论文代码的对比:
- 初始化梯度场
在这里,最开始的时候直接是将法线向量m的第一个分量当作x轴方向梯度,第二个分量作为y轴方向梯度。(正确情况(画个图):gx= − m . x m . z -\frac{m.x}{m.z} −m.zm.x,gy= − m . y m . z -\frac{m.y}{m.z} −m.zm.y)
directions[i*3+0] = static_cast<float>( 0.0f );
directions[i*3+1] = static_cast<float>( -s.x/s.z );
directions[i*3+2] = static_cast<float>( -s.y/s.z );
-
梯度限制
我没对梯度设限,论文代码是设置绝对值不超过30 -
Schwartz error
根据论文以及其他资料,我认为的施瓦茨误差的目的是确保曲面的连续性于是我所设计的施瓦茨误差是
abs( ∂ f ( x + 1 , y ) ∂ y − ∂ f ( x , y ) ∂ y + ∂ f ( x , y + 1 ) ∂ x − ∂ f ( x , y ) ∂ x \frac{\partial f(x+1,y)}{\partial y}-\frac{\partial f(x,y)}{\partial y}+\frac{\partial f(x,y+1)}{\partial x}-\frac{\partial f(x,y)}{\partial x} ∂y∂f(x+1,y)−∂y∂f(x,y)+∂x∂f(x,y+1)−∂x∂f(x,y))(基于二元函数连续的充分条件,求导顺序可交换)
论文代码中的施瓦茨误差是为了保证变化后的法线分布(gxc,gyc)和原图(gx,gy)基本一致,为
abs(gx-gxc)+abs(gy-gyc) -
对FC算法生成的高度场(曲面)的处理
自己的代码,有考虑过高度过高带来的影响,最后将高度控制在(0,2),这导致了梯度,法线被压缩到一个较小的范围,效果很不理想。
MinMaxAbValue(realimg, maxv, minv);
cv::Mat maximg = cv::Mat::zeros(realimg.size(), CV_32FC1);
maximg.setTo(maxv);
realimg += maximg;
realimg /= maxv;
return realimg;
论文代码
double vmean=0.0;
for(size_t p=0; p < size_t(surface.GetSize()); p++)
vmean+=surface(p);
vmean/=double(size_t(surface.GetSize())) ;
surface-=vmean;
- 在迭代中进行高度场限制
我没控制生成的高度场
1.论文对高度场有加过滤器控制高度
if(it < (n_of_iterations -2))
{
CFFT2D< double, double > coreFFT;
CFFT2D< CComplex< double >, double > coreiFFT;
CMatrix< CComplex< double > > SURFACE;
SURFACE = coreFFT.FFT(surface);
SURFACE = coreFFT.Shift(SURFACE);
for(size_t r=0; r < SURFACE.GetSizeRow(); r++)
for(size_t c=0; c < SURFACE.GetSizeCol(); c++)
{
SURFACE(r,c)*=filterhf(r,c);
}
SURFACE = coreFFT.iShift(SURFACE);
SURFACE = coreiFFT.iFFT(SURFACE);
for(size_t p=0; p < surface.GetSize(); p++)
surface(p)=SURFACE(p);
}
2.论文在迭代次数超过50次时,从开始到离结束还有50次的这个阶段都进行高度场最高高度的控制。
if(height_map_limit != MAX_REAL_NUMBER && it < (n_of_iterations-50))
{
//Supprime les points
CImageb imlim,imlimm;
imlim.Zero(surface.GetSize());
for(size_t i=0; i < surface.GetSize(); i++)
{
if(surface(i) > height_map_limit)
{
surface(i)=height_map_limit;
imlim(i)=true;
}
}
imlimm=MMorpho::Erosion(imlim,elem);
imlim-=imlimm;
CMatrixd kernel(3,3);
kernel.SetOne();
kernel/=double(kernel.GetFullSize());
CMatrixd surfacefiltered=surface.Convolve(kernel,MATRIX_CONV_SAME);
for(size_t i=0; i < surface.GetSize(); i++)
if(imlimm(i))
surface(i)=surfacefiltered(i);
}
- 对梯度进行划分
自己实现代码的时候没有想过要对梯度进行划分(在一块图中每个格子是一个像素,每个像素都有gxgy两个梯度(这些像素是否包含原图中的位置坐标?这些像素是标号的,所以包含位置信息),将这些像素进行以梯度为归类标准的重新划分,再得到一张格子比较大大图,但是这里面每个格子中放了数量不等的像素(这些像素还记录了他们所在大格子的位置信息))
float gradient_max = ceil(max_abs_gradient);
float gradient_limit = BBDef::Min(gradient_max,30.f);
size_t gradient_grid_size = 300;
double grid_step = (gradient_limit*2.0) / double(gradient_grid_size);
CMatrixd gx(nrow,ncol), gy(nrow,ncol);
for(size_t p=0; p < size_t(nrow*ncol); p++)
{
float vgx=Real::clamp(directions[p*3+1],-gradient_limit,gradient_limit);
float vgy=Real::clamp(directions[p*3+2],-gradient_limit,gradient_limit);
gx(p)= double(vgx);
gy(p)= double(vgy);//p可以得到原图位置信息,gx(p),gy(p)梯度信息可以得到大格子的位置信息
}
//得到大格子位置信息
double vgx=gx(pp);
double vgy=gy(pp);
size_t nrow = BBDef::Clamp(size_t((gradlim-vgy)/stepgrid),size_t(0),gridsize-1);
size_t ncol = BBDef::Clamp(size_t((gradlim-vgx)/stepgrid),size_t(0),gridsize-1);
- 寻找最佳匹配
误差小于threshold时候可以直接把原来梯度继承过来,当误差大于所设的threshold的时候,也就是经过一变换后梯度变化较大,需要去寻找更加匹配的像素。
自己的代码一是在全局搜索,且我是拿变换后的梯度去匹配原图梯度。
论文代码是在一个大格子所有误差过大的像素里面寻找最佳匹配,且拿原图梯度匹配变换梯度,这两种匹配方法肯定有总体误差的区别,我们的目的是让变换后梯度和原图尽量一致,那么由此看来就应该拿原图梯度匹配变换梯度。
理由如下:
把这个问题换成A选B的问题,A相当于原图梯度场,B相当于变换后的梯度场。A是主导(目的是和原图最匹配),所以是A选B不是B选A。 - Symmetric padding
是用来解决论文中计算梯度遗留的一个问题,图片边缘的梯度是没有的。假如不解决这个问题高度场就会随着迭代尺寸变小。具体操作看得我脑袋大(吼吼吼我发现了这是干啥了matlab里面也有这个函数函数用法)
论文代码
//!Fonction pour realiser un padding sur la matrice.
/*!
\param padsize : La taille a ajouter dans toutes les directions.
\param padtype : Le type de padding (0: zero ou valeur 1: circular, 2: replicate, 3: symmetric).
\param paddir : La direction du padding (0: avant, 1: apres, 2: les deux).
\param padValue : La valeur pour le mode valeur (padtype=0).
*/
void PadArray(const CSizeND<N>& padsize, size_t padtype=ARRAY_PADDING_VALUE, size_t paddir=2, const T& padValue=T(/*0*/))
{
BB_ASSERT(padtype<4 && paddir<3);//Test des parametres
//Copie de l'array actuel(Copy of the current array)
CArrayND<T,N> aold=(*this);
//Change la taille(Change the size)
CSizeND<N> newsize(aold.GetSize());
CSizeND<N> oldsize(aold.GetSize());
if(paddir<2)
newsize+=padsize;
else
newsize+=padsize*size_t(2);
this->Allocation(newsize);
//Gestion des pointeurs(Pointers management)
CBlock1D<size_t,N> posp;
for(size_t n=0; n<N; n++)
{
posp[n]=1;
for(size_t nn=n+1; nn<N; nn++)
posp[n]*=newsize[nn];
posp[n]*=padsize[n];
}
CBlock1D<size_t,N> newsizec;
newsizec[0]=newsize[N-1];
CBlock1D<size_t,N> oldsizec;
oldsizec[0]=oldsize[N-1];
for(size_t n=1; n<N; n++)
{
newsizec[n]=newsizec[n-1]*newsize[N-1-n];
oldsizec[n]=oldsizec[n-1]*oldsize[N-1-n];
}
//Recopie suivant toutes les directions
T* pdataNew=this->m_Value;
T* pdataOld=aold.m_Value;
internal_padding_copy(pdataOld,oldsize,pdataNew,newsize,posp,newsizec,oldsizec,padsize,0,paddir,padtype,padValue);
};
};
//!Fonction pour la recopie des anciennes valeurs pour le padding.
/*!
\param pdataold Le pointeur sur les donnees a copier.
\param sizeold La taille des donnees a copier.
\param pdatanew Le pointeur sur les nouvelles donnees.
\param sizepad La taille du padding pour chaque dimension (le nb de zero a mettre dans chaque dimension).
\param dim La dimension actuelle.
\param paddir : La direction du padding (0: avant, 1: apres, 2: les deux).
*/
void internal_padding_copy(T* &pdataold, CSizeND<N>& sizeold,
T* &pdatanew, CSizeND<N>& sizenew, CBlock1D<size_t,N>& sizepadc,
CBlock1D<size_t,N>& sizenewc, CBlock1D<size_t,N>& sizeoldc,
const CSizeND<N>& sizepad,
size_t dim, size_t paddir, size_t padtype, const T& padValue)
{//省略一段
else if(padtype==ARRAY_PADDING_SYMMETRIC)//Recopie par symetrie
{
for(size_t i=0; i<sizepadc[dim]; i++)
{
//Cherche la distance par rapport a la 1er valeur
//suivant toutes les dimensions -dim.
int dist[N];
for(size_t n=0; n<N; n++)
dist[n]=0;
dist[0] = int(i%sizenewc[0]-sizepad[N-1]);
if(dist[0]<0)
dist[0]=-dist[0]-1;
dist[0]=(dist[0]/int(sizeold[N-1]))%2==0? dist[0]%int(sizeold[N-1]) : int(sizeold[N-1])-1-dist[0]%int(sizeold[N-1]);
int j=int(i);
for(int n=1; n<int(N)-int(dim); n++)
{
j/=int(sizenew[N-n]);
dist[n]=j%int(sizenew[N-n-1])-int(sizepad[N-n-1]);
if(dist[n]<0)
dist[n]=-dist[n]-1;
dist[n]=(dist[n]/int(sizeold[N-n-1]))%2==0 ? dist[n]%int(sizeold[N-n-1]) : int(sizeold[N-n-1])-1-dist[n]%int(sizeold[N-n-1]);
dist[n]*=int(sizeoldc[n-1]);
}
//Calcul la position dans la matrice old
int pos=0;
for(size_t n=0; n<N; n++)
pos+=dist[n];
(*pdatanew)=(*(pdataold+pos));
pdatanew++;
}
}//省略另一段}