真实感图像渲染系列：渐进式光子映射算法(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

#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;
}