月度归档： 2018年7月

这是真实感图像渲染系列的第七篇文章。

如果将渲染修改为并行，那么渲染效率将提升约8倍（我的电脑是八核）。非常非常有诱惑力。而且基于OpenMP框架的并行化甚至还非常简单！

平台编译

使用OpenMP最麻烦的地方是框架安装问题，倘若本机的g++版本支持-fopenmp ,那么恭喜你，你已经完成了90%的工作了。倘若本机g++不支持该编译选项，那么需要在网上搜索对应的安装方法。安装方法不同系统都有较大的差别，故不在这里细讲。

OpenMP程序并行化指令都是用#pragma 实现的，也就是说倘若编译器不支持，也不会产生什么问题。

OpenMP简单使用

这是一个最简单的for语句，对应于PPM算法中的光子映射。

for (int i=0;i<total_photons;i++) {
    PhotonTracing(...);
}

该for语句满足如下条件：

• 循环起始位置和终结位置确定
• 不同的循环没有依赖性，

满足这样条件的for语句，理论上，将循环不同的i并行处理是没有问题的。OpenMP使用了最简单的预编译指令帮助我们实现for语句的并行，修改为

#include <omp.h>

#pragma omp parallel for
for (int i=0;i<total_photons;i++) {
    PhotonTracing(...);
}

OpenMP临界区设置

#pragma omp critical 编译器将自动将#pragma 指令后面紧接着的语句块进行并行。

一般来说，并行语句块内使用的变量，都是读写局部变量以及读全局变量，这时，不会存在什么问题，但是总有特殊情况，我们多个并行的循环需要修改同一个全局变量。这是可能导致错误的。如PPM算法中，计算“视点”，需要将结果保存在全局数组view_pts 里面。这时，多个线程并行修改会导致一些问题（诸如SegmentationFault）。这时，使用#pragma omp critical 可以强制后面紧接着的那个语句块通过“加锁”保证同时只有一个线程在运行该代码块。

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

 

这是真实感图像渲染系列的第六篇文章。

动机

一幅图片渲染时间少则一分钟，多则半天，等到渲染完成再看效果未免也太佛系了。况且诸如PPM算法亮度调节什么的一般在前两轮的时候，就已经大概知道了，并不需要渲染完。因此我们希望写一个实时显示渲染效果的类。

实现

显示部分在技术上没有任何难度，就是OpenCV中的imshow而已，如果必要，每次imshow之前，使用imwrite备份图片到文件。

我们需要一边显示一遍渲染，这个就需要多线程了。由于opencv的限制，imshow和waitKey只能在主线程调用，因此我们把图像渲染放在子线程里面去。这里线程的创建我使用了pthreads库。

具体细节就看代码啦。init函数传入一个Color类数组的指针，而PaintBoard能够实时监测该指针指向内容的变化，并且显示出来。（需要用户在窗口按回车键刷新）

 

#ifndef PAINT_BOARD_H
#define PAINT_BOARD_H
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

#include "../core/color.h"
using namespace std;

class PaintBoard {
	private:
		cv::Mat *image;
		Color **board;
		int sizX,sizY;
	public:
		PaintBoard();
		~PaintBoard();
		void init(int sizX,int sizY,Color** _board);
		void update();
		void display();
		void save();
		void save_raw();
};

#endif

#include "paint_board.h"
#include "glog/logging.h"
#include <iostream>
using namespace cv;

PaintBoard::PaintBoard() {
	sizX = -1,sizY = -1;
	image = NULL;
}
PaintBoard::~PaintBoard() {
}

void PaintBoard::init(int sizX,int sizY,Color** _board) {
	this->sizX = sizX;
	this->sizY = sizY;
	image = new Mat(this->sizX,this->sizY,CV_64FC3,Scalar(1.0,1.0,1.0));
	board = _board;
}

void PaintBoard::update() {
	assert(image);
	for (int i=0;i<sizX;i++) {
		for (int j=0;j<sizY;j++) {
			//Color& c = (*board)[(sizX-i-1)*sizY+j];
			Color& c = (*board)[i*sizY+j];
			Vec3f v(c.getR(),c.getG(),c.getB());
			image->at<Vec3d>(i,j) = v;
		}
	}
}
void PaintBoard::display() {
	imshow("picture", *image);
	waitKey(0);
}
void PaintBoard::save_raw() {
    FILE* fout = fopen("result.ppm","w");
    fprintf(fout,"P3 %d %d\n",sizY,sizX);
    fprintf(fout,"255\n");
    for (int i=0;i<sizX;i++)
    {
        for (int j=0;j<sizY;j++)
        {
			Color c = (*board)[i*sizY+j]*255;
            fprintf(fout,"% 3d % 3d % 3d",int(c.getR()),int(c.getG()),int(c.getB()));
        }
        fprintf(fout,"\n");
    }
    fclose(fout);
}
void PaintBoard::save() {
    vector<int>compression_params;
    compression_params.push_back(IMWRITE_PNG_COMPRESSION);
    compression_params.push_back(9);
    imwrite("result.bmp", (*image)*255);
}

#include <iostream>
#include "display/paint_board.h"
#include "json/json.h"
#include "render/raytracing.h"
#include "render/progressive_photon_mapping.h"
#include <unistd.h>
#include <fstream>
#include <pthread.h>

bool show_graph_flag = true;

void* show_graph(void* params) {
	Render *PPM = (Render*)params;
	PPM->run();
	show_graph_flag = false;
	return NULL;
}


int main(int argc, char** args)
{
	std::ifstream ifs("configures/config.json");
	Json::CharReaderBuilder reader;
	Json::Value root;
	JSONCPP_STRING errs;
	Json::parseFromStream(reader, ifs, &root, &errs);

	Render *PPM = new ProgressivePhotonMapping();
	//Render *PPM = new RayTracing();
	PaintBoard PB;
	PPM->accept(root);
	PPM->registerPaintBoard(&PB);

	pthread_t watcher;
	pthread_create(&watcher,NULL,show_graph,PPM);

	while (show_graph_flag) {
		PB.update();
		PB.save();
		PB.display();
		usleep(100000);
	}

	show_graph_flag = false;
	pthread_join(watcher,NULL);

	PB.update();
	PB.save();
}

 

这是真实感图像渲染系列的第五篇文章。提前声明，这篇文章讲述的算法之和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;
}

 

这是真实感图像渲染系列的第四篇文章。

算法概述

光线追踪算法可以算是图像渲染中的Hello World了，不过这个Hello World可没有那么好写。通过前两篇文章，我们已经知道了渲染算法包含的核心要素，如颜色，物体等表示。

算法本身网上已经有很多详细的教程，多说无益，本质上RayTracing就是最直接的通过视线反推光路，并忽视一切不容易追踪的光路的算法。

光线追踪本身可以处理的光路需要满足两个限制：

1. 只经过了一次漫反射面
2. 漫反射面是作为摄像机到光源光路中最后一个界面存在

我们可以通过一个玻璃球在地上是否会聚焦光线形成光斑来判断一个算法是不是普通的RayTracing。这里，漫反射面（地面）并没有作为光路的最后一个面（中间还插入了玻璃球的投射面）。光线追踪算法会在玻璃球下产生一个黑色的阴影。

算法步骤

下面，我将以最简练的方式概述RayTracing的算法。

对于屏幕上的每一个像素，我们可以通过摄像机(camera)找到一条视线光路，RayTracing算法目的就是得到这条光路的颜色。当然，倘若光路直接射到了光源内部（即与光源产生碰撞），那么颜色就是光源颜色。否则，倘若光路射到了无穷远处，那么我们将提前设置好的背景光颜色作为实现颜色。

多数情况，视线既不会与光源碰撞，也不会到达无穷远处。视线沿着这个光路行进，在物体表面可能由于“折射”，“反射”等而分叉，我们递归得处理折射，反射后视线的颜色。这里，递归总深度，光线强度皆可作为停止递归的信号。

最终，部分视线将到达一个漫反射平面，这时，光线亮度设置为能够直接照亮该位置的光源亮度。对于复杂光源，诸如面光源，我们需要在光源处随机取点，计算出光源中多少光可以不受遮挡到达反射面。

基于hash的边缘超采样

通过基本的光线追踪算法，我们可以得到如下的图片

可以发现，球体的边沿存在较严重的锯齿效应。通过基于hash的边缘超采样后，效果如下

锯齿处的颜色变味了两边颜色的混合，从而在肉眼看来消除了锯齿。

在RayTracing中，实现的函数调用可以看做是一棵树结构，每一个节点都会碰撞一个物体。在每一次碰撞后，我们可以维护光路hash值。通过如下方法迭代出来的hash值，只要中途光路的折射顺序有任何细微差别，都会导致hash值的变化。我们对于摄像机每一个像素的hash值进行比较，倘若一个像素四周hash值都和他相同，则其不在边缘上。否则，我们对该像素(x,y)进行超采样。其颜色为(x,y),(x+1/3,y+1/3),(x-1/3,y+1/3),…,(x-1/3,y-1/3)这九个像素位置对应颜色的平均值。

在实际运行中，超采样的运行时间相对会比较长。毕竟稍微复杂一点的场景，会有很多边缘点，而每个边缘点的重采样时间是原来的9倍。

源代码

ray_tracing.h

#ifndef RAYTRACING_H
#define RAYTRACING_H

#include <vector>
#include "../core/light.h"
#include "../core/object.h"
#include "../core/color.h"
#include "scene.h"
#include "json/json.h"
#include "render.h"

class RayTracing : private Scene,public Render{
	private:
		unsigned **hash_table;
		int max_depth;
		int shade_quality;
		int spec_power;
		int start_rows;
		int bazier_quality;
	public:
		//INIT
		RayTracing();
		~RayTracing();
		virtual void accept(const Json::Value& val);


		//RayTrace
		Color calcReflection(const Object& obj, const Collision& obj_coll, int depth, unsigned& hash);
		Color calcRefraction(const Object& obj, const Collision& obj_coll, int depth, unsigned& hash);
		Color calcDiffusion(const Object& obj, const Collision& obj_coll);

		Color rayTrace(const Vector& rayO, const Vector& rayD, int depth, unsigned& hash);

		void run();
};

#endif

ray_tracing.cpp

#include "raytracing.h"
#include <string>
#include <iostream>

#ifdef USE_OPENMP
#include <omp.h>
#pragma message("The openmp is enable!")
#endif

RayTracing::RayTracing() {
}

RayTracing::~RayTracing() {
}

void RayTracing::accept(const Json::Value& val) {
	Render::accept(val);
	Scene::accept(val,rx,ry);

	max_depth = val["max_depth"].asInt();
	shade_quality = val["shade_quality"].asInt();
	spec_power = val["spec_power"].asInt();
	start_rows = val["start_rows"].asInt();
	bazier_quality = val["bazier_quality"].asInt();
	board = new Color[rx*ry];
	for (int i=0;i<rx*ry;i++)
		board[i] = Color(1,1,1);
}


Color RayTracing::calcReflection(const Object& obj, const Collision& obj_coll, int depth, unsigned&hash) {
	if (depth > max_depth) {
		hash = hash * 17 + 233;
		return bg_color;
	}
	hash = hash * 17 + obj.getHash();
	Vector resO,resD;
	obj_coll.reflection(resO,resD);
	return (rayTrace(resO,resD,depth,hash) * obj.getMaterial().refl).adjust();
}

Color RayTracing::calcRefraction(const Object& obj, const Collision& obj_coll, int depth, unsigned& hash) {
	if (depth > max_depth) {
		hash = hash * 19 + 233;
		return bg_color;
	}
	hash = hash * 19 + obj.getHash();
	Vector rayD;
	Vector rayO;
	Color acol;
	obj_coll.refraction(rayO,rayD);
	Color rcol = rayTrace(rayO,rayD,depth,hash);
	if (obj_coll.face)
		rcol = rcol*obj.getMaterial().refr;
	if (!obj_coll.face)
		acol = (obj.getAbsorb().getColor(0,0)*-obj_coll.dist).exp();
	else
		acol = Color(1,1,1);

	return (rcol*acol).adjust();
}

Color RayTracing::calcDiffusion(const Object& obj, const Collision& obj_coll) {
	Color color = obj.getColor(obj_coll.C);
	Color ret = color * bg_color *(1-obj.getMaterial().diff-obj.getMaterial().spec);
	for (auto &lgt : lights) {
		double shade = lgt->getShade(obj_coll.getSurfaceC(),objects,shade_quality);
		ret += color * lgt->getColor(lgt->getCenter()) * shade * obj.getMaterial().diff;
		double spec_ratio = (lgt->getCenter() - obj_coll.C).unit() ^ obj_coll.N.unit();
		if (spec_ratio > 0)
			ret += color * lgt->getColor(lgt->getCenter()) * shade * pow(spec_ratio,spec_power) * obj.getMaterial().spec;
	}
	return ret.adjust();
}

Color RayTracing::rayTrace(const Vector& rayO, const Vector& rayD, int depth, unsigned& hash) {
	Collision obj_coll,lgt_coll;
	const Object* obj = findCollidedObject(rayO,rayD,obj_coll);
	const Light* lgt = findCollidedLight(rayO,rayD,lgt_coll);
	if (!obj && !lgt) {
		hash = hash * 17 + 235;
		return bg_color;
	}else if (!obj || (obj && lgt && obj_coll.dist > lgt_coll.dist)) {
		hash = hash * 17 + lgt->getHash();
		return lgt->getColor(lgt->getCenter());
	}else {
		Color ret;
		if (obj->getMaterial().refl>feps)
			ret += calcReflection(*obj,obj_coll,depth+1,hash);
		if (obj->getMaterial().diff>feps || obj->getMaterial().spec>feps)
			ret += calcDiffusion(*obj,obj_coll);
		if (obj->getMaterial().refr>feps)
			ret += calcRefraction(*obj,obj_coll,depth+1,hash);
		return ret.adjust();
	}
}

void RayTracing::run() {
	hash_table = new unsigned*[rx];
	for (int i=0;i<rx;i++)
		hash_table[i] = new unsigned[ry];
#ifdef USE_OPENMP
#pragma omp parallel
#pragma omp for schedule(dynamic,3)
#endif
	for (int _i=start_rows;_i<rx;_i++){
		int i = _i;
#ifdef USE_OPENMP
		if (!i)
			std::cout<<omp_get_num_threads()<<std::endl;
#endif
		printf("Render row #%d\n",i);
		for (int j=0;j<ry;j++) {
			Vector rayO,rayD;
			camera->getRay(i,j,rayO,rayD);
			Color cc = rayTrace(rayO,rayD,0,hash_table[i][j]);
			board[i*ry+j] = cc;
		}
	}
#ifdef USE_OPENMP
#pragma omp for schedule(dynamic,5)
#endif
	for (int i=0;i<rx;i++) {
		printf("Resample row #%d\n",i);
		for (int j=0;j<ry;j++) {
			bool flag = false;
			if (i!=0 && hash_table[i][j] != hash_table[i-1][j])
				flag = true;
			if (i!=rx-1 && hash_table[i][j] != hash_table[i+1][j])
				flag = true;
			if (j!=0 && hash_table[i][j] != hash_table[i][j-1])
				flag = true;
			if (j!=ry-1 && hash_table[i][j] != hash_table[i][j+1])
				flag = true;
			if (flag) {
				Color res;
				for (int k1=-1; k1<=1; k1++) {
					for (int k2=-1; k2<=1; k2++) {
						unsigned hash = 0;
						Vector rayO,rayD;
						camera->getRay(i+k1/3.0,j+k2/3.0,rayO,rayD);
						res += rayTrace(rayO,rayD,0,hash);
					}
				}
				res = res*(1.0/9);
				board[i*ry + j] = res;
			}
		}
	}
	for (int i=0;i<rx;i++) {
		delete[] hash_table[i];
	}
	delete[] hash_table;
}

 

这是真实感图像渲染系列的第三篇文章。

场景初始化的必要性

不同于普通的编程作业，图像渲染的场景结构非常复杂，比如同样是物体，一个球体需要读入中心坐标以及半径，而平面则只需要一个定位坐标以及法线方向就能够确定。如何初始化一个场景，为每一个参数赋予正确的初始值是非常重要的。

在参考过往代码发现场景的初始化无非两种方式：一种是将场景作为一个类写在源代码中；另一种是将场景写在文本文件中，对于所有需要初始化的类写一个接口来读取文件。在我的代码中，使用了第二种初始化方式。就编译效率上来说，第二种方式可以防止每次修改模型都重新编译。这在后期可以优化出很多时间。

场景初始化方式

一方面，我使用json格式储存场景，json格式的读取在C++中确实不算简单，需要安装jsoncpp库。jsoncpp库提供了一个数据类型Json::Value ，这个数据类型支持整数下表访问（用于访问列表），以及字符串下表访问（用于访问字典），以及各种转换函数如val.asDouble() 等。最终场景的json文件附在本文末尾。

jsoncpp的简单使用可以参考 通过jsoncpp库的CharReaderBuilder解析字符串

我使用伪·visitor设计模式初始化场景，对于每一个需要初始化的类申明一个成员函数void accept(const Json::Value& val); 接受一个Json::Value 用于初始化自身。

例如，对于vector来说，accept函数如下。

void Vector :: accept(const Json::Value& val) {
	if (!val.isMember("x") || !val.isMember("y") || !val.isMember("z"))
		std::cout<<"The vector not found..."<<std::endl;
	x = val["x"].asDouble();
	y = val["y"].asDouble();
	z = val["z"].asDouble();
}

而一个包含vector的sphere类，可以拥有如下accept函数

void Sphere::accept(const Json::Value& val) {
	Object::accept(val);
	O.accept(val["center"]);
	radius = val["radius"].asDouble();
}

可见，至少在安装了库之后，json的读取是足够简单的。

为场景添加功能

场景本身会被后续不同的渲染算法所使用，为了增加代码复用性，我们可以把渲染算法公共的一些函数放在场景中，例如求一个光线最近碰撞到的物体的函数findCollidedObject 等。最终场景类如下

class Scene {
	protected:
		Camera *camera;
		std::vector<Light*> lights;
		std::vector<Object*> objects;
		Color bg_color;
	public:
		Scene();
		~Scene();
		void accept(const Json::Value& val,int _rx,int _ry);
		const Object* findCollidedObject(const Vector &rayO,const Vector& rayD,Collision& coll);
		const Light* findCollidedLight(const Vector& rayO,const Vector& rayD, Collision& coll);
};

后续的光线追踪等代码，都会直接继承scene来使用scene的这些接口。

源代码

场景json文件

{
 "rx": 512,
 "ry": 1024,
 "_rx": 5,
 "_ry": 10,
 "photon_num": 20000,
 "total_round": 80,
 "total_brightness": 1310.0,
 "round_decay": 0.95,
 "initial_r": 5,
 "r_decay": 0.8,
 "camera": {
  "type": "default",
  "position": {
   "x": 20,
   "y": -10,
   "z": 0
  },
  "dx": {
   "x": -10,
   "y": 0,
   "z": 0
  },
  "dy": {
   "x": 0,
   "y": 20,
   "z": 0
  },
  "origin": {
   "x": 40,
   "y": 0,
   "z": -50
  },
  "fdepth": 0.02
 },
 "lights": [
  {
   "type": "point_light",
   "name": "light1",
   "brightness": 40,
   "position": {
    "x": 12,
    "y": -5,
    "z": 27
   },
   "texture": {
    "type": "pure",
    "r": 1,
    "g": 1,
    "b": 1
   },
   "material": {
    "refl": 0,
    "diff": 0,
    "spec": 0,
    "refr": 0.0,
    "refr_k": 1.0
   }
  },
  {
   "type": "point_light",
   "name": "light1",
   "brightness": 50,
   "position": {
    "x": 37,
    "y": 0.6,
    "z": -45
   },
   "texture": {
    "type": "pure",
    "r": 1,
    "g": 1,
    "b": 1
   },
   "material": {
    "refl": 0,
    "diff": 0,
    "spec": 0,
    "refr": 0.0,
    "refr_k": 1.0
   }
  }
 ],
 "objects": [
  {
   "type": "sphere",
   "name": "marble_ball",
   "center": {
    "x": 3,
    "y": 0,
    "z": 30
   },
   "radius": 1.7,
   "texture": {
    "type": "picture",
    "filename": "materials/marble.bmp",
    "rx": 2,
    "ry": 3.14
   },
   "absorb": {
    "type": "pure",
    "r": 0.5,
    "g": 0.7,
    "b": 0.5
   },
   "material": {
    "refl": 0.3,
    "diff": 0.8,
    "spec": 0.6,
    "refr": 0.0,
    "refr_k": 1.3
   }
  },
  {
   "type": "sphere",
   "name": "marble_ball",
   "center": {
    "x": 3,
    "y": -6,
    "z": 29
   },
   "radius": 1.7,
   "texture": {
    "type": "picture",
    "filename": "materials/marble2.jpg",
    "rx": 2,
    "ry": 3.14
   },
   "absorb": {
    "type": "pure",
    "r": 0.5,
    "g": 0.7,
    "b": 0.5
   },
   "material": {
    "refl": 0.3,
    "diff": 0.8,
    "spec": 0.6,
    "refr": 0.0,
    "refr_k": 1.3
   }
  },
  {
   "type": "sphere",
   "name": "ball",
   "center": {
    "x": 3,
    "y": 0,
    "z": 22
   },
   "radius": 2,
   "texture": {
    "type": "pure",
    "r": 0.7,
    "g": 0,
    "b": 0.7
   },
   "absorb": {
    "type": "pure",
    "r": 0.4,
    "g": 0.1,
    "b": 0.01
   },
   "material": {
    "refl": 0.1,
    "diff": 0.0,
    "spec": 0.6,
    "refr": 0.9,
    "refr_k": 1.3
   }
  },
  {
   "type": "sphere",
   "name": "marble_ball",
   "center": {
    "x": 35,
    "y": 1,
    "z": -40
   },
   "radius": 0.2,
   "texture": {
    "type": "picture",
    "filename": "materials/marble2.jpg",
    "rx": 2,
    "ry": 3.14
   },
   "absorb": {
    "type": "pure",
    "r": 0.5,
    "g": 0.7,
    "b": 0.5
   },
   "material": {
    "refl": 0,
    "diff": 1.0,
    "spec": 0.6,
    "refr": 0,
    "refr_k": 1.3
   }
  },
  {
   "type": "bazier_curves",
   "name": "cup",
   "position": {
    "x": 0,
    "y": 0,
    "z": 30
   },
   "ctrl_pts": [
    [
     [
      0.0,
      1.1,
      1.2,
      2.0
     ],
     [
      0.04,
      0.04,
      -0.0,
      -0.0
     ]
    ],
    [
     [
      2.0,
      2.2,
      2.0,
      1.9
     ],
     [
      -0.0,
      -0.0,
      0.1,
      0.15
     ]
    ],
    [
     [
      1.9,
      1.73,
      0.4,
      0.2
     ],
     [
      0.15,
      0.23,
      -0.2,
      0.4
     ]
    ],
    [
     [
      0.2,
      -0.1,
      3.8,
      3.8
     ],
     [
      0.4,
      1.3,
      1.4,
      2.0
     ]
    ],
    [
     [
      3.8,
      3.8,
      3.6,
      3.6
     ],
     [
      2.0,
      2.2,
      2.2,
      2.0
     ]
    ],
    [
     [
      3.6,
      3.6,
      0.25,
      0.0
     ],
     [
      2.0,
      1.9,
      1.3,
      1.3
     ]
    ]
   ],
   "texture": {
    "type": "pure",
    "r": 0.7,
    "g": 0.7,
    "b": 0.3
   },
   "absorb": {
    "type": "pure",
    "r": 0.1,
    "g": 0.1,
    "b": 0.5
   },
   "material": {
    "refl": 0.0,
    "diff": 0.0,
    "spec": 0.0,
    "refr": 0.9,
    "refr_k": 1.3
   }
  },
  {
   "type": "plane",
   "name": "floor",
   "position": {
    "x": 0,
    "y": -10,
    "z": 25
   },
   "dx": {
    "x": 0,
    "y": 0,
    "z": 30
   },
   "dy": {
    "x": 0,
    "y": 20,
    "z": 0
   },
   "texture": {
    "type": "picture",
    "filename": "materials/floor.jpg",
    "rx": 8,
    "ry": 8
   },
   "material": {
    "refl": 0.4,
    "diff": 0.6,
    "spec": 0.0,
    "refr": 0.0,
    "refr_k": 1.0
   },
   "border": false
  }
 ],
 "bg_color": {
  "type": "pure",
  "r": 0.0,
  "g": 0.0,
  "b": 0.0
 },
 "max_depth": 20,
 "max_jump": 20,
 "start_rows": 0,
 "start_cols": 0
}

 

这是真实感图像渲染系列的第二篇文章。

颜色

颜色，向量与物体是各种渲染算法的基础。

颜色定义为三个取值在0-1之间的浮点数，分别表示RGB分量。由于分量的值恰好为0-1之间的浮点数，颜色不仅可以理解为每种颜色的强度，还可以理解为每种颜色光线的穿透率。换句话来说，这样定义的颜色同时支持乘法和加法。

class Color
{
	private:
		double r,g,b;
	public:
		Color();
		Color(double r,double g,double b);
		Color operator += (const Color&a);
		void accept(const Json::Value& val);
		double getR()const {return r;}
		double getG()const {return g;}
		double getB()const {return b;}
		friend Color operator +(const Color& a,const Color &b);
		friend Color operator *(const Color& a,double k);
		friend Color operator *(const Color& a,const Color& b);
		Color exp()const;
		std::string description()const;
		Color adjust()const;
};

特别的，对于透明物体的折射，光线穿过物体，存在一定的吸收，该吸收率可以看做是距离的指数函数。即，如在红光吸收率0.1的物体中穿过2单位距离，则射出的红光应该为原来射入红光强度的.

向量

向量的定义，也是三个0-1之间的浮点数，表示xyz坐标分量。向量同时具有表征方向和表征位置的功能，坐标本身在某种意义上也是坐标原点O加上某个表示方向的向量。

一个光线，我们就可以用两个向量和分别表示他的起始位置以及方向，是一个单位向量，某点在光线上当且仅当，且。

class Vector
{
	double x,y,z;
	public:
		Vector();
		Vector(double x,double y,double z);
		Vector(const Eigen::Vector3d vec);
		double getX() const {return x;}
		double getY() const {return y;}
		double getZ() const {return z;}
		virtual void accept(const Json::Value &sl);
		friend Vector operator * (const Vector &a,const Vector &b);
		friend double operator ^ (const Vector &a,const Vector &b);
		friend Vector operator + (const Vector &a,const Vector &b);
		friend Vector operator - (const Vector &a,const Vector &b);
		friend Vector operator - (const Vector &a);
		friend Vector operator * (double k,const Vector &a);
		friend Vector operator * (const Vector &a,double k);
		friend bool operator == (const Vector &a,const Vector &b);
		friend Vector operator /(const Vector& a,double k);
		friend Vector each_min(const Vector& v1,const Vector& v2);
		friend Vector each_max(const Vector& v1,const Vector& v2);
		Vector unit()const;
		Vector reverse()const;
		double len()const;
		double sqrlen()const;
		bool isUnit()const;
		std::string description()const;
		Eigen::Vector3d eigen()const;
		static Vector randomVectorOnSphere();
};

 

物体

由于C++的多态机制，物体本身是一个抽象类，只需要支持：

• 判断和某个光线的交点
• 判断某个位置对应的贴图颜色
• 提取出诸如漫反射率等表面材质的性质
• 询问物体hash值（用于抗锯齿等）

struct Material {
	double refl;//reflection ratio
	double diff;//diffusion ratio
	double spec;//high light diffusion
	double refr;//refraction ratio
	double refr_k;
	void accept(const Json::Value& val);
};

class Object {
	protected:
		unsigned hash;
		Material material;
		std::string name;
		Texture* texture;
		Texture* absorb;
	public:
		Object();
		virtual ~Object();
		virtual void accept(const Json::Value& val);
		unsigned getHash()const {return hash;}

		const Material& getMaterial()const {return material;}

		virtual bool collideWith(const Vector& rayO, const Vector& rayD,Collision& coll) = 0;
		virtual Color getColor(const Vector &pos)const = 0;
		std::string getName()const{return name;}
		const Texture& getAbsorb()const;
};

碰撞

碰撞是渲染里面非常核心也很复杂的一个模块，我们使用一个专门的类记录碰撞的信息，实际上，对于碰撞，可以通过碰撞位置，碰撞法线方向和如何光线方向完全表述。有了这三个量，我们可以简单的计算出这个碰撞的反射、折射方向（其中折射方向需要物体折射率）。这里讲一个小的注意点，就是反射，折射光线的起始点不应该都设置为，反射需要保证光线的起始位置和入社光线位置同侧，而折射恰恰相反，由于精度误差，可能实际在物体表面的某一侧而我们用作为新光线的起点是不妥的。解决方案很简单，对于反射，我们使用，折射则使用作为光线的其实位置。

class Object;

struct Collision {
	Vector C;//The Center of Collision
	Vector N;//The normal of the plane
	Vector I;//The direction of reflaction
	double dist;
	bool face;
	Object* belongs;
	std::string description()const;
	void refraction(Vector& resO, Vector& resD)const;
	void reflection(Vector& resO, Vector& resD)const;
	void diffusion(Vector& resO, Vector& resD)const;
	void diffusion_hl(Vector& resO, Vector& resD)const;
	Vector getSurfaceC()const;
	Vector getBackfaceC()const;
};

 

 

这是真实感图像渲染系列的第一篇文章。

这是图形学课程的一次大作业，基础要求是实现光线追踪算法(Ray Tracing)，包括反射，折射以及漫反射，以及参数曲面求交(Bazier Curves)，在此基础上，我增加了渐进式光子映射算法(Progressive Photon Mapping)以及景深效果(Depth of Focus). 不过由于本人毫无美术细胞，所以场景设计一团糟，只是粗略体现出了每一个算法的现象而已。结果如图

包含两个含贴图的球体，一个看起来像水果盘的贝塞尔曲线，一个透视的玻璃球，以及一个距离摄像头太近的，对焦失败的球。

工程源代码可以在 https://github.com/mhy12345/RayTracing-ProgressivePhotonMaping 找到。

其中，该工程依赖了glog（实际上master版本不依赖，但是需要他进行一些必要的调试信息输出），eigen（用于矩阵运算），opencv（用户实时显示渲染结果以及bmp格式保存）以及jsoncpp（用于解析json文件）。后面三个库都需要提前安装在电脑中。同时，工程需要你的编译器支持-fopenmp  选项，例如，在我的笔记本上，编译器 clang-omp 支持该选项，编译命令即为

cmake -DCMAKE_CXX_COMPILER=/usr/local/bin/clang-omp++ 
-DCMAKE_C_COMPILER=/usr/local/bin/clang-omp ..

根目录下的./serial 和./main 分别是渲染程序的串行版本和并行版本。