这是真实感图像渲染系列的第五篇文章。提前声明,这篇文章讲述的算法之和PPM在思想上相似,细节上由于笔者也没有仔细研究论文,所以可能会有比较大的出入。
算法概述
在光线追踪算法 中,我们提到了光线追踪实际上无法成功捕获到一部分光路。随之而来的是透明物体下方出现错误的阴影。
渐进式光子映射算法(Progressive Photon Mapping),尝试分别从摄像机和光源两边一起投放射线。两边的射线在某漫反射面上相聚。RayTracing和PPM的关系挺像深度优先搜索和双向深度优先搜索的关系。
摄像机出发有若干视线,这些视线落在漫反射面上形成视点。由于视线本身存在归属的像素,所以每一个视点也存在归属像素。假设该视点位置得到了光强,那么该光强将累加至对应像素。
与此同时,光源处随机发出若干光子,光子经过折射,反射,以及随机漫反射后,也会在若干位置碰撞漫反射界面。我们希望光子会对碰撞位置周围半径R的所有视点有一定的光强贡献。称之为“送光”。
从摄像机出发视点的计算,和RayTracing视线的计算相似,唯一的区别是RayTracing在视线碰撞漫反射面后,通过采样光源得到此处的光强,而PPM在得到的所有视点后,使用一颗KD树进行储存便于后续的“送光”。
从光源处发的光子,其起始位置、方向皆随机选取。与视线相比,唯一的不同是光线并不会终止于一个漫反射面,在漫反射面上,光子仍然会随机选取方向反射,并且,在对于该漫反射面碰撞点周围R距离内的所有视点,都进行“送光”。
对于整体稳定的考虑
PPM算法中,光子是随机发射的,由于光线追踪本身速度并不快,所以有限时间内,实际可以发射光子数是有限的。另一方面,如果我们场景是一个“点光源”和一个“地平面”,很容易发现,点光源随机发射的可以达到地面可见区域的光子更少。我们现在需要考虑的,就是如何合理利用这些足够少的光子,使得渲染出来的图片尽量均匀。
下图是R比较小导致送光不均匀的例子,观察地面上“斑驳”的光影:
光子本身需要给半径R周围的所有视点送光,一方面,提高R有助于使得送光范围更大,少数离散的光子也能渲染出均匀的图片。另一方面,将送光光强定义为半径R从中心到两边线性递减的函数,可以消除边界处光强骤减的问题。
对于局部高亮的考虑
过大的R导致局部细节无法突出,下图是R较大的情形,可以发现橙色透明玻璃球的聚焦并不算成功。
因此最优的方法应该是将渲染分成若干不同的round,R随着渲染的进行而逐渐变小。这样,前期图像趋向于亮度均匀,后期,细节被强调出来。
对于整体光强的考虑
每个光子究竟需要携带多少光强这是一个非常难以衡量的问题。由之前的讨论,我们意识到了光子的发射需要区分round,不同round对应的R不同,换句话来说,前面的round和后面的round同等重要。但是,倘若前面round和后面round光子拥有相同的光强,这会下面的问题:
由于像素本身颜色是对应所有视点接收到的光强总和,在RayTracing算法中,我们限制光子的“折射”“反射”皆存在系数,最终的光强一定是小于等于(1,1,1)的。在PhotonMapping由于是随机化算法,我们是无法保证这一点的,沿用RayTracing的光强算法,我们得到的像素光强将是一个无限接近于(0,0,0)的值,这并不是我们想要的。同时,如果将所有像素光强乘以一个值,也会导致问题,比如一块正常的蓝色(0.1,0.1,1),如果不幸被乘以了系数10,最终得到的是(1,1,10) = (1,1,1),蓝色就变成白色了。这个问题可能导致下图所示的结果,不是我们想要的。
解决思路很巧妙,核心是调整每一轮的光强,我是用的方式是将每一轮光强设置为调和级数。即最终图片光强为. 实际上这是一个巧妙地想法,因为后期的round期望对整体亮度不产生改变,而对局部亮度改变,恰好调和级数之和不收敛,但是单纯取出奇数或偶数项的话,
他是收敛的假装他就突然收敛了。意味着倘若一个位置在之后的每一次都有很多光子碰撞,则该位置倾向于变量,而倘若之后光子碰撞频率较低的话,该位置亮度倾向于收敛。
实际上,究竟多少光强是一件很灵活的事,例如我发现高亮部分不明显,那么我可以定义每一轮初始光子光强随轮数增强,增强的系数等同于半径衰减系数等等。即而
.
KD树的建立
我工程中使用的KD树是基于维度轮换分割的,写起来非常方便。这里单独讲一讲KD树编写中的一些技巧:
- KD树建树期间,按照指定维度分割出较小的一半和较大的一半。使用std::sort函数排序固然可以,但是存在更加高效的算法:std::nth_elements(begin,begin+k,end)可以提取序列第k大,将小于第k大数的放在前面,其余放在后面。将k设为n/2,则KD树的分割部分直接通过调用nth_elements完成。
- KD树通过对子树维护点集包围盒,实现快速判断一个询问点是否在包围盒中。实际上,很多人即使得到了包围盒以及一个点的坐标,仍然无法得到一个高效的算法判断点距离包围盒的距离。这里介绍一个比较快捷的方式:
- 令包围盒用
定义,询问点
- 取
,
,
- 答案
满足:
- 令包围盒用
渲染效果图
源代码
progressive_photon_mapping.h
#ifndef PROGRESSIVE_PHOTON_MAPPING_H #define PROGRESSIVE_PHOTON_MAPPING_H #include <vector> #include "../core/light.h" #include "../core/camera.h" #include "../core/object.h" #include "../core/color.h" #include "scene.h" #include "json/json.h" #include "render.h" class ViewPoint { private: Vector C,N; Color color; double strength; int x,y; public: ViewPoint(){} ViewPoint(Vector pos,Vector N,Color color,double stgh,int x,int y); friend class ProgressivePhotonMapping; }; class ProgressivePhotonMapping : private Scene, public Render{ private: struct KDTreeNode { public: ViewPoint value; KDTreeNode* lch; KDTreeNode* rch; int split_dim; Vector bd_max,bd_min; }; template<int dim> class ViewPointComparer { public: bool operator()(const ViewPoint& p1,const ViewPoint& p2) { if (dim == 0) return p1.C.getX() < p2.C.getX(); if (dim == 1) return p1.C.getY() < p2.C.getY(); if (dim == 2) return p1.C.getZ() < p2.C.getZ(); } }; KDTreeNode* view_root; int max_depth; int start_rows; int start_cols; int bazier_quality; int total_round; int photon_num; double total_brightness; double round_decay; double initial_r; int max_jump; std::vector<ViewPoint> view_pts; Color* bg_pic; public: //INIT ProgressivePhotonMapping(); ~ProgressivePhotonMapping(); virtual void accept(const Json::Value& val); void run(); void RayTracing(const Vector& rayO, const Vector& rayD,Color rayC,int xx,int yy,int depth,double lambda,double dist); void photonTracing(const Vector& rayO, const Vector& rayD, const Color& rayC,int depth, double r, double lambda, double dist,int diff_count); void buildKDTree(KDTreeNode* &now,std::vector<ViewPoint> &lst,int l = -1,int r = -1,int dim = 0); void releaseKDTree(KDTreeNode* &node); void queryKDTree(KDTreeNode* node, vector<const ViewPoint*> &result, const Vector& pos, double r); }; #endif
progressive_photon_mapping.cpp
#include "progressive_photon_mapping.h" #include <string> #include "../core/object.h" #include <iostream> #include <ctime> #include <algorithm> using namespace std; ViewPoint::ViewPoint(Vector pos,Vector N,Color color,double stgh,int x,int y) : C(pos),N(N),color(color),strength(stgh),x(x),y(y) {} ProgressivePhotonMapping::ProgressivePhotonMapping() { view_root = NULL; } ProgressivePhotonMapping::~ProgressivePhotonMapping() { } void ProgressivePhotonMapping::accept(const Json::Value& val) { Render::accept(val); Scene::accept(val,rx,ry); max_depth = val["max_depth"].asInt(); start_rows = val["start_rows"].asInt(); start_cols = val["start_cols"].asInt(); bazier_quality = val["bazier_quality"].asInt(); total_round = val["total_round"].asInt(); photon_num = val["photon_num"].asInt(); total_brightness = val["total_brightness"].asDouble(); round_decay = val["round_decay"].asDouble(); initial_r = val["initial_r"].asDouble(); max_jump = val["max_jump"].asInt(); board = new Color[rx*ry]; for (int i=0;i<rx*ry;i++) board[i] = Color(0,0,0); } void ProgressivePhotonMapping::run() { std::cout<<"Start SPPM algorithm..."<<std::endl; bg_pic = new Color[rx*ry]; int* samples_count = new int[rx*ry]; for (int i=0;i<rx*ry;i++) samples_count[i] = 0; double current_r = initial_r; double energy = 1.0 / log(total_round); double current_e = energy; for (int iter = 0; iter < total_round; iter++) { for (int i=0;i<rx*ry;i++) bg_pic[i] = Color(0,0,0); std::cout<<"Round <"<<iter<<"> ... "<<std::endl; std::cout<<"View Point generate..."<<std::endl; std::cout<<"Current R = "<<current_r<<std::endl; std::cout<<"Current E = "<<current_e<<std::endl; view_pts.clear(); srand(time(0)); #ifdef USE_OPENMP #pragma omp parallel num_threads(8) #pragma omp for #endif for (int i=start_rows;i<rx;i++) { if (i % 100 == 0) cout<<"Render row "<<i<<endl; for (int j=start_cols;j<ry;j++) { double _i = i + (rand()*1.0/RAND_MAX-.5)*1; double _j = j + (rand()*1.0/RAND_MAX-.5)*1; Vector rayO; Vector rayD; camera->getRay(_i,_j,rayO,rayD); if (i == 0 && j == 0) cout<<rayO.description()<<" "<<rayD.description()<<endl; RayTracing(rayO,rayD,Color(1,1,1),i,j,0,current_e,0); if (!i && !j) cout<<view_pts[0].C.description()<<endl; samples_count[i*ry+j] ++; } } std::cout<<"Build KDTree..."<<std::endl; buildKDTree(view_root,view_pts); std::cout<<"Photon Scattering"<<std::endl; double brightness = 0; for (auto &lgt : lights) brightness += lgt->getBrightness(); for (auto &lgt : lights) { int lgt_emit_count = photon_num*lgt->getBrightness()/brightness; #ifdef USE_OPENMP #pragma omp for #endif for (int i=0;i<lgt_emit_count;i++) { if (i%1000 == 0) cout<<"Sample Photon "<<i<<" / "<<lgt_emit_count<<endl; Vector rayO,rayD; lgt->randomlyEmit(rayO,rayD); photonTracing(rayO,rayD,lgt->getColor(lgt->getCenter()),0,current_r,total_brightness/photon_num,0,0); } } std::cout<<"Release KDTree..."<<std::endl; releaseKDTree(view_root); for (int i=0;i<rx;i++) for (int j=0;j<ry;j++) board[i*ry+j] += bg_pic[i*ry+j]*(1.0/samples_count[i*ry+j]); current_r *= round_decay; current_e /= round_decay; } } void ProgressivePhotonMapping::photonTracing(const Vector& rayO, const Vector& rayD, const Color& rayC, int depth, double r, double lambda, double dist,int diff_count) { if (depth > max_jump) return ; Collision obj_coll; const Object* obj; obj = findCollidedObject(rayO,rayD,obj_coll); if (!obj)return; dist += obj_coll.dist; if (obj->getMaterial().refl > feps) { Vector resO,resD; obj_coll.reflection(resO,resD); photonTracing(resO,resD,rayC,depth+1,r,lambda*obj->getMaterial().refl,dist,diff_count); } if (obj->getMaterial().refr > feps) { Color acol; if (!obj_coll.face) acol = (obj->getAbsorb().getColor(0,0)*-obj_coll.dist).exp(); else acol = Color(1,1,1); Vector resO,resD; obj_coll.refraction(resO,resD); photonTracing(resO,resD,rayC*acol,depth+1,r,lambda*obj->getMaterial().refr,dist,diff_count); } if (obj->getMaterial().diff > feps) { vector<const ViewPoint*> vps; queryKDTree(view_root,vps,obj_coll.getSurfaceC(),r); for (auto&w : vps) { if ((w->N ^ rayD)<-feps) { Color res = rayC * w->color * pow(((r-(w->C-obj_coll.getSurfaceC()).len())/r),2) * lambda * (1.0/(diff_count+1)) * w->strength; bg_pic[w->x * ry + w->y] += res; } } Vector resO,resD; obj_coll.diffusion(resO,resD); photonTracing(resO,resD,rayC*obj->getColor(obj_coll.C),depth+1,r,lambda*obj->getMaterial().diff,dist+obj_coll.dist,diff_count+1); } } void ProgressivePhotonMapping::RayTracing(const Vector& rayO,const Vector& rayD,Color rayC,int xx,int yy,int depth,double lambda,double dist) { if (lambda < 1e-9 || depth > max_depth) return ; Collision obj_coll,lgt_coll; const Object* obj = findCollidedObject(rayO,rayD,obj_coll); const Light* lgt = findCollidedLight(rayO,rayD,lgt_coll); if ((!obj && !lgt)) { bg_pic[xx*ry+yy] += bg_color*lambda; }else if ((!obj && lgt) || (obj && lgt && obj_coll.dist > lgt_coll.dist)) { bg_pic[xx*ry+yy] += lgt->getColor(lgt_coll.C)*lambda; }else { if (obj->getMaterial().refr > feps) { Color acol; if (!obj_coll.face) acol = (obj->getAbsorb().getColor(0,0)*-obj_coll.dist).exp(); else acol = Color(1,1,1); Vector resO,resD; obj_coll.refraction(resO,resD); RayTracing(resO,resD,rayC*acol,xx,yy,depth+1,lambda*obj->getMaterial().refr,dist+obj_coll.dist); } if (obj->getMaterial().refl > feps) { Vector resO,resD; obj_coll.reflection(resO,resD); RayTracing(resO,resD,rayC,xx,yy,depth+1,lambda*obj->getMaterial().refl,dist+obj_coll.dist); } if (obj->getMaterial().diff > feps) { #pragma omp critical view_pts.push_back( ViewPoint( obj_coll.getSurfaceC(), obj_coll.N, rayC*obj->getColor(obj_coll.C), lambda*obj->getMaterial().diff, xx,yy ) ); } } } void ProgressivePhotonMapping::buildKDTree(KDTreeNode* &node, std::vector<ViewPoint> &lst, int l,int r,int dim) { if (l == -1 && r == -1) l = 0, r = lst.size(); if (l >= r) return; int mid = (l+r)>>1; switch(dim) { case 0: nth_element(lst.begin()+l,lst.begin()+mid,lst.begin()+r,ViewPointComparer<0>()); case 1: nth_element(lst.begin()+l,lst.begin()+mid,lst.begin()+r,ViewPointComparer<1>()); case 2: nth_element(lst.begin()+l,lst.begin()+mid,lst.begin()+r,ViewPointComparer<2>()); } node = new KDTreeNode(); node->value = lst[mid]; //std::cout<<node->value.C.description()<<std::endl; node->lch = node->rch = NULL; node->split_dim = dim; node->bd_max = node->value.C; node->bd_min = node->value.C; buildKDTree(node->lch,lst,l,mid,(dim+1)%3); if (node->lch) { node->bd_max = each_max(node->bd_max,node->lch->bd_max); node->bd_min = each_min(node->bd_min,node->lch->bd_min); } buildKDTree(node->rch,lst,mid+1,r,(dim+1)%3); if (node->rch) { node->bd_max = each_max(node->bd_max,node->rch->bd_max); node->bd_min = each_min(node->bd_min,node->rch->bd_min); } } void ProgressivePhotonMapping::queryKDTree(KDTreeNode* node, vector<const ViewPoint*> &result, const Vector& pos, double r) { double dx,dy,dz; if (pos.getX() <= node->bd_max.getX() && pos.getX() >= node->bd_min.getX()) dx = 0; else dx = min(abs(pos.getX()-node->bd_max.getX()),abs(pos.getX()-node->bd_min.getX())); if (pos.getY() <= node->bd_max.getY() && pos.getY() >= node->bd_min.getY()) dy = 0; else dy = min(abs(pos.getY()-node->bd_max.getY()),abs(pos.getY()-node->bd_min.getY())); if (pos.getZ() <= node->bd_max.getZ() && pos.getZ() >= node->bd_min.getZ()) dz = 0; else dz = min(abs(pos.getZ()-node->bd_max.getZ()),abs(pos.getZ()-node->bd_min.getZ())); if (dx*dx + dy*dy + dz*dz >r*r) return ; if ((node->value.C-pos).len()<=r) result.push_back(&(node->value)); if (node->lch)queryKDTree(node->lch,result,pos,r); if (node->rch)queryKDTree(node->rch,result,pos,r); } void ProgressivePhotonMapping::releaseKDTree(KDTreeNode* &node) { if (!node)return; releaseKDTree(node->lch); releaseKDTree(node->rch); delete node; }