真实感图像渲染系列：渐进式光子映射算法(Progressive Photon Mapping)

这是真实感图像渲染系列的第五篇文章。提前声明，这篇文章讲述的算法之和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)，蓝色就变成白色了。这个问题可能导致下图所示的结果，不是我们想要的。

解决思路很巧妙，核心是调整每一轮的光强，我是用的方式是将每一轮光强设置为调和级数。即最终图片光强为1 \times round_1 + \frac{1}{2} \times round_2 + \frac{1}{3} \times round_3 + \dots + \frac{1}{n} \times round_n. 实际上这是一个巧妙地想法，因为后期的round期望对整体亮度不产生改变，而对局部亮度改变，恰好调和级数之和不收敛，但是单纯取出奇数或偶数项的话，他是收敛的假装他就突然收敛了。意味着倘若一个位置在之后的每一次都有很多光子碰撞，则该位置倾向于变量，而倘若之后光子碰撞频率较低的话，该位置亮度倾向于收敛。

实际上，究竟多少光强是一件很灵活的事，例如我发现高亮部分不明显，那么我可以定义每一轮初始光子光强随轮数增强，增强的系数等同于半径衰减系数等等。即R' = R \times \alpha而E' = E\times \frac{1}{\alpha}.

KD树的建立

我工程中使用的KD树是基于维度轮换分割的，写起来非常方便。这里单独讲一讲KD树编写中的一些技巧：

  • KD树建树期间，按照指定维度分割出较小的一半和较大的一半。使用std::sort函数排序固然可以，但是存在更加高效的算法：std::nth_elements(begin,begin+k,end)可以提取序列第k大，将小于第k大数的放在前面，其余放在后面。将k设为n/2，则KD树的分割部分直接通过调用nth_elements完成。
  • KD树通过对子树维护点集包围盒，实现快速判断一个询问点是否在包围盒中。实际上，很多人即使得到了包围盒以及一个点的坐标，仍然无法得到一个高效的算法判断点距离包围盒的距离。这里介绍一个比较快捷的方式：
    • 令包围盒用maxX,minX,maxY,minY,maxZ,minZ定义，询问点(X,Y,Z)
    • 取dx = max(0,|X-maxX|,|X-minX|),dy = max(0,|Y-maxY|,|Y-minY|),dz = max(0,|Z-maxZ|,|Z-minZ|)
    • 答案dist满足：dist = \sqrt{dx^2+dy^2+dz^2}

渲染效果图

源代码

progressive_photon_mapping.h

progressive_photon_mapping.h
C++
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

progressive_photon_mapping.cpp
C++
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;
}

 

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据。