games101作业7光线追踪含多线程和微表面提高

对于光线追踪进行综合运用。

光线与三角形求交

其它的emit那些，现在先不用管，后面看看作用是什么。

inline Intersection Triangle::getIntersection(Ray ray)
{Intersection inter;if (dotProduct(ray.direction, normal) > 0)//光线从里面打，无交点return inter;double u, v, t_tmp = 0;Vector3f pvec = crossProduct(ray.direction, e2);double det = dotProduct(e1, pvec);if (fabs(det) < EPSILON)return inter;double det_inv = 1. / det;Vector3f tvec = ray.origin - v0;u = dotProduct(tvec, pvec) * det_inv;if (u < 0 || u > 1)return inter;Vector3f qvec = crossProduct(tvec, e1);v = dotProduct(ray.direction, qvec) * det_inv;if (v < 0 || u + v > 1)return inter;t_tmp = dotProduct(e2, qvec) * det_inv;if(t_tmp<0) return inter;inter.happened=true;inter.coords=ray(t_tmp);inter.normal=normal;inter.distance=t_tmp;inter.obj=this;inter.m=this->m;// TODO find ray triangle intersection// bool happened;// Vector3f coords;// Vector3f tcoords;// Vector3f normal;// Vector3f emit;// double distance;// Object* obj;// Material* m;return inter;
}

判断包围盒与光线是否相交

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg) const
{Vector3f t1=(pMin-ray.origin)*invDir;Vector3f t2=(pMax-ray.origin)*invDir;Vector3f tMin=Vector3f::Min(t1,t2);Vector3f tMax=Vector3f::Max(t1,t2);float tInter=std::max(tMin.x,std::max(tMin.y,tMin.z));float tExit=std::min(tMax.x,std::min(tMax.y,tMax.z));if(tExit>=0&&tInter<=tExit){return true;}// invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division// dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic// TODO test if ray bound intersectsreturn false;
}

BVH的查找

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{Intersection res;if(node==nullptr||!node->bounds.IntersectP(ray,ray.direction_inv,{0,0,0})){return res;}if(node->left==nullptr&&node->right==nullptr){return node->object->getIntersection(ray);}Intersection left=getIntersection(node->left,ray);Intersection right=getIntersection(node->right,ray);res=left.distance<=right.distance?left:right;return res;// TODO Traverse the BVH to find intersection}

实现光线追踪

参考

改动

首先是调用castRay的Renderer::Render(const Scene& scene)

它对同一个pixel采样了16次，求了平均。(这里我不是很懂，理论上不应该是对单个像素的分小像素来取样吗？）

void Renderer::Render(const Scene& scene)
{std::vector<Vector3f> framebuffer(scene.width * scene.height);float scale = tan(deg2rad(scene.fov * 0.5));float imageAspectRatio = scene.width / (float)scene.height;Vector3f eye_pos(278, 273, -800);int m = 0;// change the spp value to change sample ammountint spp = 16;std::cout << "SPP: " << spp << "\n";for (uint32_t j = 0; j < scene.height; ++j) {for (uint32_t i = 0; i < scene.width; ++i) {// generate primary ray directionfloat x = (2 * (i + 0.5) / (float)scene.width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;Vector3f dir = normalize(Vector3f(-x, y, 1));for (int k = 0; k < spp; k++){framebuffer[m] += scene.castRay(Ray(eye_pos, dir), 0) / spp;  }m++;}UpdateProgress(j / (float)scene.height);}UpdateProgress(1.f);// save framebuffer to fileFILE* fp = fopen("binary.ppm", "wb");(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);for (auto i = 0; i < scene.height * scene.width; ++i) {static unsigned char color[3];color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));fwrite(color, 1, 3, fp);}fclose(fp);    
}

Object.hpp

增加了获得面积的方式

virtual float getArea()=0;
virtual void Sample(Intersection &pos, float &pdf)=0;
virtual bool hasEmit()=0;

get_random_float()的改动（针对windows）

inline float get_random_float()
{static std::random_device dev;static std::mt19937 rng(dev());static std::uniform_real_distribution<float> dist(0.f, 1.f); // distribution in range [0，1]return dist(rng);
}

本次实验给出的与框架匹配的伪代码

解决

首先我们需要判断射线是否打到物体上，如果没有打到物体上则返回0,0,0.如果打到物体上判断是否是光源，如果是光源，则返回光源的rgb的值，如果不是光源则考虑打到的是物体，则需要使用上面的伪代码来求渲染方程。

如果没直接打到光源，则求渲染方程的步骤为：

求BRDF
如果是直接光照，使用sampleLight对光源进行采样求交点射线，求光源的pdf。利用蒙特卡洛积分求直接光照。

如果是间接光照，因为当前只有一条射线，求该射线的随机一个方向，随机的方向与场景中的物体求交点，求得的交点判断是否为光源，不是则计算间接光照。根据俄罗斯赌盘理论进行判断是否求该射线的影响。然后求BRDF，然后递归求该交点物体的间接光照，利用已有的函数求pdf。最后得到间接光照。（一条射线）

Vector3f Scene::castRay(const Ray &ray, int depth) const
{Intersection inter=intersect(ray);Vector3f dir={0.0,0.0,0.0};Vector3f indir=(0.0,0.0,0.0);if(!inter.happened){return dir;}if(inter.m->hasEmission()){return inter.m->getEmission();//如果与之相交的物体材质属于光照，则返回光照}//获取物体的信息// bool happened;// Vector3f coords;// Vector3f tcoords;// Vector3f normal;// Vector3f emit;// double distance;// Object* obj;// Material* m;Vector3f p=inter.coords;Vector3f p_n=inter.normal.normalized();Material* obj_material=inter.m;Vector3f w0=ray.direction;Intersection ans;float pdf_light;sampleLight(ans,pdf_light);//对光照进行采样Vector3f x=ans.coords;Vector3f ws=(x-p).normalized();Vector3f NN=ans.normal.normalized();Vector3f emit=ans.emit;float d=(x-p).norm();Ray objLight(p,ws);Intersection pas=intersect(objLight);//判断物体和光照之间是否有物体遮挡if(abs(pas.distance-ans.distance)<0.001){Vector3f eval=obj_material->eval(w0,ws,p_n);//求PRDFfloat cosinep=dotProduct(NN,-ws);float cosine=dotProduct(p_n,ws);dir+=emit*cosinep*cosine*eval/(d*d)/pdf_light;//没有物体遮挡，计算直接光照}//其他物体间接光照// 其它物体间接光照的方式就是根据欧罗斯赌盘，判断你的这根光线目前是否发射出去。float P_RR=get_random_float();if(P_RR<RussianRoulette){//小于就发射Vector3f wi=obj_material->sample(w0,p_n).normalized();Ray obj_to_obj(p,wi);Intersection obj_to_obj_inter=intersect(obj_to_obj);if(obj_to_obj_inter.happened&&!obj_to_obj_inter.m->hasEmission()){Vector3f eval=obj_material->eval(w0,wi,p_n);float cosine=dotProduct(p_n,wi);float pdf_to_obj=obj_material->pdf(w0,wi,p_n);//求出射方向的概率indir+=castRay(obj_to_obj,depth+1)*eval*cosine/pdf_to_obj/RussianRoulette;//求间接光照}}return indir+dir;//Get x, ws , NN , emit from inter// TO DO Implement Path Tracing Algorithm here
}

一次采样（spp==1）的结果

16次（spp==16）

提高1多线程

这里的思路是多线程去批量处理像素点。提高参考

这里的思路就是简单的对高度进行批量处理，相当于同时处理24个像素点。

void Renderer::Render(const Scene& scene)
{std::vector<Vector3f> framebuffer(scene.width * scene.height);std::atomic_int process=0;//共享变量float scale = tan(deg2rad(scene.fov * 0.5));float imageAspectRatio = scene.width / (float)scene.height;Vector3f eye_pos(278, 273, -800);int m = 0;// change the spp value to change sample ammountint spp = 16;std::cout << "SPP: " << spp << "\n";float thread_num=24;std::thread th[25];int per_height=scene.height/thread_num;std::function<void(uint32_t,uint32_t)> func=[&](uint32_t height,uint32_t height_next){for(uint32_t j=height;j<height_next;j++){for(uint32_t i=0;i<scene.width;i++){float x = (2 * (i + 0.5) / (float)scene.width - 1) *imageAspectRatio * scale;float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;Vector3f dir = normalize(Vector3f(-x, y, 1));for(int k=0;k<spp;k++){framebuffer[static_cast<int>(j*scene.width)+i]+=scene.castRay(Ray(eye_pos,dir),0)/spp;}}process++;UpdateProgress(process/(float)scene.height);}};for(int i=0;i<thread_num;i++){th[i]=std::thread(func,per_height*i,per_height*(i+1));}for(int i=0;i<thread_num;i++){th[i].join();}UpdateProgress(1.f);// save framebuffer to fileFILE* fp = fopen("binary.ppm", "wb");(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);for (auto i = 0; i < scene.height * scene.width; ++i) {static unsigned char color[3];color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));fwrite(color, 1, 3, fp);}fclose(fp);    
}

修改cmake文件，链接thread库

cmake_minimum_required(VERSION 3.10)
project(RayTracing)set(CMAKE_CXX_STANDARD 17)add_executable(RayTracing main.cpp Object.hpp Vector.cpp Vector.hpp Sphere.hpp global.hpp Triangle.hpp Scene.cppScene.hpp Light.hpp AreaLight.hpp BVH.cpp BVH.hpp Bounds3.hpp Ray.hpp Material.hpp Intersection.hppRenderer.cpp Renderer.hpp)
find_package(Threads)
target_link_libraries (${PROJECT_NAME} ${CMAKE_THREAD_LIBS_INIT})

同时还需要改变自己虚拟机的运行内核

速度提升的还是比较明显的

提高2 Microfacet

参考1 参考2 计算BRDF

微表面模型（microfacet model） 是一种基于物理的局部光照模型，它假设物体的表面是凹凸不平的，宏观的表面由许多微小的平面（即微表面）构成，光线在每个微平面上发生理想镜面反射或折射。

计算法线分布函数 D

微平面的法线分布函数D(m)描述了微观表面上的表面法线m的统计分布。（法线分布函数：某个点上，法线指向指定方向上的概率。）它表示在单位面积上，法线向量与半程向量 hhh 一致的微表面面积的比例。

业界较为主流的法线分布函数是GGX（Trowbridge-Reitz），因为具有更好的高光长尾。

α为粗糙度∈[0,1]
n为宏观平面的法线
h为微观平面法线，即半程向量（入射光线和出射光线的）

当α=1时，D=1/π,光线向半球面均匀发散

当α=0时，只有当n=h时，D≠=0

代码

float D_GGX(const Vector3f& n,const Vector3f& h,float rough){float a2=roughness*roughness;float _p=M_PI*pow(pow(dotProduct(normalize(N),normalize(H)),2.0)*(a2-1.0f)+1.0f,2.0f);float _pas=std::max(_p,0.0000001f);float D_res=a2/_pas;return D_res;
}

几何函数

目前较为常用的是其中最为简单的形式，分离遮蔽阴影（Separable Masking and Shadowing Function）

该形式将几何项G分为两个独立的部分：光线方向（light）和视线方向（view），并对两者用相同的分布函数来描述。

其中UE4的方案是Schlick-GGX，即基于Schlick近似，将k映射为k=α/2，去匹配GGX Smith方程：

G(n,v,l,k)返回遮蔽阴影函数的计算结果。

N为宏观表面法线
V为光线反射方向（可理解为视口的方向）
L为光线进入的反方向（可理解为光源的方向）
roughness为粗糙度。

实际上就是分别使用入射和出射方向求Gdirection然后相乘

float GeometrySchlickGGX(const Vector3f& normal,const Vector3f& v,const float k){Vector3f _normal=normalize(normal);Vector3f _v=normalize(v);float _pas=std::max(dotProduct(_normal,_v),0.0f);float _mpas=_pas*(1.0f-k)+k;return _pas/_mpas;
}
float GeometrySmith(const Vector3f& n,const Vector3f& wi,const Vector3f wo,float rough){float k=(rough+1.0f)*(rough+1.0f)/8.0f;return GeometrySchlickGGX(n,wi,k)*GeometrySchlickGGX(n,wo,k);
}

Fresnel方程 F（菲涅尔项）

float Fresnel(const float& n1,const float& n2,const Vector3f& wi,const Vector3f n,float& F){float r0=pow((n1-n2)/(n1+n2),2.0f);float _pas=std::max(dotProduct(normalize(wi),normalize(n)),0.0f);return r0+(1-r0)*pow((1-_pas),5.0f);
}

但是在工程中，实现的是完整版的

void fresnel(const Vector3f &I, const Vector3f &N, const float &ior, float &kr) const{float cosi = clamp(-1, 1, dotProduct(I, N));float etai = 1, etat = ior;if (cosi > 0) {  std::swap(etai, etat); }// Compute sini using Snell's lawfloat sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));// Total internal reflectionif (sint >= 1) {kr = 1;}else {float cost = sqrtf(std::max(0.f, 1 - sint * sint));cosi = fabsf(cosi);float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));kr = (Rs * Rs + Rp * Rp) / 2;}// As a consequence of the conservation of energy, transmittance is given by:// kt = 1 - kr;}

计算微表面材质

为了能量守恒，漫反射归到折射光照中，高光为ks也就是镜面反射。

Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){switch(m_type){case DIFFUSE:{// calculate the contribution of diffuse   modelfloat cosalpha = dotProduct(N, wo);if (cosalpha > 0.0f) {Vector3f diffuse = Kd / M_PI;return diffuse;}elsereturn Vector3f(0.0f);break;}case Microfacet://微表面材质{float cosalpha = dotProduct(N, wo);if(cosalpha>0.0f){float F;Vector3f specular;specular=caculate(wo,wi,N,F);Vector3f diffuse=1.0f/M_PI;float _kd=1-F;return specular*Ks+_kd*Kd*diffuse;}else{return Vector3f(0.0f);}break;}}
}

其它在Material中用到diffuse材质的部分，也需要加上微表面材质

float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){switch(m_type){case DIFFUSE:case Microfacet:{// uniform sample probability 1 / (2 * PI)if (dotProduct(wo, N) > 0.0f)return 0.5f / M_PI;elsereturn 0.0f;break;}}
}

Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){switch(m_type){case DIFFUSE:case Microfacet:{// uniform sample on the hemispherefloat x_1 = get_random_float(), x_2 = get_random_float();float z = std::fabs(1.0f - 2.0f * x_1);float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);return toWorld(localRay, N);break;}}
}

Material* red = new Material(Microfacet, Vector3f(0.0f));red->Kd = Vector3f(0.63f, 0.065f, 0.05f);Material* green = new Material(Microfacet, Vector3f(0.0f));green->Kd = Vector3f(0.14f, 0.45f, 0.091f);Material* white = new Material(Microfacet, Vector3f(0.0f));white->Kd = Vector3f(0.725f, 0.71f, 0.68f);Material* light = new Material(Microfacet, (8.0f * Vector3f(0.747f+0.058f, 0.747f+0.258f, 0.747f) + 15.6f * Vector3f(0.740f+0.287f,0.740f+0.160f,0.740f) + 18.4f *Vector3f(0.737f+0.642f,0.737f+0.159f,0.737f)));light->Kd = Vector3f(0.65f);

换成兔子模型参考

int main(int argc, char** argv)
{// Change the definition here to change resolutionScene scene(784, 784);Material* red = new Material(DIFFUSE, Vector3f(0.0f));red->Kd = Vector3f(0.63f, 0.065f, 0.05f);Material* green = new Material(DIFFUSE, Vector3f(0.0f));green->Kd = Vector3f(0.14f, 0.45f, 0.091f);Material* white = new Material(DIFFUSE, Vector3f(0.0f));white->Kd = Vector3f(0.725f, 0.71f, 0.68f);Material* light = new Material(DIFFUSE, (8.0f * Vector3f(0.747f+0.058f, 0.747f+0.258f, 0.747f) + 15.6f * Vector3f(0.740f+0.287f,0.740f+0.160f,0.740f) + 18.4f *Vector3f(0.737f+0.642f,0.737f+0.159f,0.737f)));light->Kd = Vector3f(0.65f);Material* whitem = new Material(Microfacet,Vector3f(0.0f));whitem->Ks=Vector3f(0.45f,0.45f,0.45f);whitem->Kd=Vector3f(0.3f,0.3f,0.25f);MeshTriangle bunny("../models/bunny/bunny.obj",whitem);MeshTriangle floor("../models/cornellbox/floor.obj", white);// MeshTriangle shortbox("../models/cornellbox/shortbox.obj", white);// MeshTriangle tallbox("../models/cornellbox/tallbox.obj", white);MeshTriangle left("../models/cornellbox/left.obj", red);MeshTriangle right("../models/cornellbox/right.obj", green);MeshTriangle light_("../models/cornellbox/light.obj", light);scene.Add(&floor);// scene.Add(&shortbox);// scene.Add(&tallbox);scene.Add(&left);scene.Add(&right);scene.Add(&light_);scene.Add(&bunny);scene.buildBVH();Renderer r;auto start = std::chrono::system_clock::now();r.Render(scene);auto stop = std::chrono::system_clock::now();std::cout << "Render complete: \n";std::cout << "Time taken: " << std::chrono::duration_cast<std::chrono::hours>(stop - start).count() << " hours\n";std::cout << "          : " << std::chrono::duration_cast<std::chrono::minutes>(stop - start).count() << " minutes\n";std::cout << "          : " << std::chrono::duration_cast<std::chrono::seconds>(stop - start).count() << " seconds\n";return 0;
}

但是这样是看不到兔子的

需要对兔子进行位置和大小变换

只需要将兔子每个顶点进行放缩和位移

MeshTriangle bunny("../models/bunny/bunny.obj",whitem,Vector3f(2000.0f),Vector3f(300.0f,0.0f,300.0f));

MeshTriangle(const std::string& filename, Material *mt = new Material(),Vector3f Scale=Vector3f(1.0f,1.0f,1.0f),Vector3f move=Vector3f(0.0f,0.0f,0.0f)){objl::Loader loader;loader.LoadFile(filename);area = 0;m = mt;assert(loader.LoadedMeshes.size() == 1);auto mesh = loader.LoadedMeshes[0];Vector3f min_vert = Vector3f{std::numeric_limits<float>::infinity(),std::numeric_limits<float>::infinity(),std::numeric_limits<float>::infinity()};Vector3f max_vert = Vector3f{-std::numeric_limits<float>::infinity(),-std::numeric_limits<float>::infinity(),-std::numeric_limits<float>::infinity()};for (int i = 0; i < mesh.Vertices.size(); i += 3) {std::array<Vector3f, 3> face_vertices;for (int j = 0; j < 3; j++) {auto vert = Vector3f(mesh.Vertices[i + j].Position.X,mesh.Vertices[i + j].Position.Y,mesh.Vertices[i + j].Position.Z);vert=vert*Scale+move;face_vertices[j] = vert;min_vert = Vector3f(std::min(min_vert.x, vert.x),std::min(min_vert.y, vert.y),std::min(min_vert.z, vert.z));max_vert = Vector3f(std::max(max_vert.x, vert.x),std::max(max_vert.y, vert.y),std::max(max_vert.z, vert.z));}triangles.emplace_back(face_vertices[0], face_vertices[1],face_vertices[2], mt);}bounding_box = Bounds3(min_vert, max_vert);std::vector<Object*> ptrs;for (auto& tri : triangles){ptrs.push_back(&tri);area += tri.area;}bvh = new BVHAccel(ptrs);}