这是真实感图像渲染系列的第五篇文章。提前声明，这篇文章讲述的算法之和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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 |
#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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 |
#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; } |