由法线分布得到高度场(表面)算法实现

算法来自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} yf(x+1,y)yf(x,y)+xf(x,y+1)xf(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++;
                }
            }//省略另一段}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值