物流网站给做软件,网站优化工具升上去,网监备案网站更换域名,石家庄限号对于光线追踪进行综合运用。
光线与三角形求交
其它的emit那些#xff0c;现在先不用管#xff0c;后面看看作用是什么。
inline Intersection Triangle::getIntersection(Ray ray)
{Intersection inter;if (dotProduct(ray.direction, normal) 0)//光线从里面打现在先不用管后面看看作用是什么。
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_tmp0) return inter;inter.happenedtrue;inter.coordsray(t_tmp);inter.normalnormal;inter.distancet_tmp;inter.objthis;inter.mthis-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::arrayint, 3 dirIsNeg) const
{Vector3f t1(pMin-ray.origin)*invDir;Vector3f t2(pMax-ray.origin)*invDir;Vector3f tMinVector3f::Min(t1,t2);Vector3f tMaxVector3f::Max(t1,t2);float tInterstd::max(tMin.x,std::max(tMin.y,tMin.z));float tExitstd::min(tMax.x,std::min(tMax.y,tMax.z));if(tExit0tIntertExit){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(x0),int(y0),int(z0)], 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(nodenullptr||!node-bounds.IntersectP(ray,ray.direction_inv,{0,0,0})){return res;}if(node-leftnullptrnode-rightnullptr){return node-object-getIntersection(ray);}Intersection leftgetIntersection(node-left,ray);Intersection rightgetIntersection(node-right,ray);resleft.distanceright.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::vectorVector3f 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_distributionfloat dist(0.f, 1.f); // distribution in range [01]return dist(rng);
}
本次实验给出的与框架匹配的伪代码 解决
首先我们需要判断射线是否打到物体上如果没有打到物体上则返回0,0,0.如果打到物体上判断是否是光源如果是光源则返回光源的rgb的值如果不是光源则考虑打到的是物体则需要使用上面的伪代码来求渲染方程。
如果没直接打到光源则求渲染方程的步骤为
求BRDF如果是直接光照使用sampleLight对光源进行采样求交点射线求光源的pdf。利用蒙特卡洛积分求直接光照。 如果是间接光照因为当前只有一条射线求该射线的随机一个方向随机的方向与场景中的物体求交点求得的交点判断是否为光源不是则计算间接光照。根据俄罗斯赌盘理论进行判断是否求该射线的影响。然后求BRDF然后递归求该交点物体的间接光照利用已有的函数求pdf。最后得到间接光照。一条射线 Vector3f Scene::castRay(const Ray ray, int depth) const
{Intersection interintersect(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 pinter.coords;Vector3f p_ninter.normal.normalized();Material* obj_materialinter.m;Vector3f w0ray.direction;Intersection ans;float pdf_light;sampleLight(ans,pdf_light);//对光照进行采样Vector3f xans.coords;Vector3f ws(x-p).normalized();Vector3f NNans.normal.normalized();Vector3f emitans.emit;float d(x-p).norm();Ray objLight(p,ws);Intersection pasintersect(objLight);//判断物体和光照之间是否有物体遮挡if(abs(pas.distance-ans.distance)0.001){Vector3f evalobj_material-eval(w0,ws,p_n);//求PRDFfloat cosinepdotProduct(NN,-ws);float cosinedotProduct(p_n,ws);diremit*cosinep*cosine*eval/(d*d)/pdf_light;//没有物体遮挡计算直接光照}//其他物体间接光照// 其它物体间接光照的方式就是根据欧罗斯赌盘判断你的这根光线目前是否发射出去。float P_RRget_random_float();if(P_RRRussianRoulette){//小于就发射Vector3f wiobj_material-sample(w0,p_n).normalized();Ray obj_to_obj(p,wi);Intersection obj_to_obj_interintersect(obj_to_obj);if(obj_to_obj_inter.happened!obj_to_obj_inter.m-hasEmission()){Vector3f evalobj_material-eval(w0,wi,p_n);float cosinedotProduct(p_n,wi);float pdf_to_objobj_material-pdf(w0,wi,p_n);//求出射方向的概率indircastRay(obj_to_obj,depth1)*eval*cosine/pdf_to_obj/RussianRoulette;//求间接光照}}return indirdir;//Get x, ws , NN , emit from inter// TO DO Implement Path Tracing Algorithm here
}
一次采样spp1的结果 16次spp16 提高1多线程
这里的思路是多线程去批量处理像素点。提高参考
这里的思路就是简单的对高度进行批量处理相当于同时处理24个像素点。
void Renderer::Render(const Scene scene)
{std::vectorVector3f framebuffer(scene.width * scene.height);std::atomic_int process0;//共享变量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_num24;std::thread th[25];int per_heightscene.height/thread_num;std::functionvoid(uint32_t,uint32_t) func[](uint32_t height,uint32_t height_next){for(uint32_t jheight;jheight_next;j){for(uint32_t i0;iscene.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 k0;kspp;k){framebuffer[static_castint(j*scene.width)i]scene.castRay(Ray(eye_pos,dir),0)/spp;}}process;UpdateProgress(process/(float)scene.height);}};for(int i0;ithread_num;i){th[i]std::thread(func,per_height*i,per_height*(i1));}for(int i0;ithread_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 一致的微表面面积的比例。
业界较为主流的法线分布函数是GGXTrowbridge-Reitz因为具有更好的高光长尾。 α为粗糙度∈[0,1]n为宏观平面的法线h为微观平面法线即半程向量入射光线和出射光线的
当α1时D1/π,光线向半球面均匀发散
当α0时只有当nh时D≠0
代码
float D_GGX(const Vector3f n,const Vector3f h,float rough){float a2roughness*roughness;float _pM_PI*pow(pow(dotProduct(normalize(N),normalize(H)),2.0)*(a2-1.0f)1.0f,2.0f);float _passtd::max(_p,0.0000001f);float D_resa2/_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 _normalnormalize(normal);Vector3f _vnormalize(v);float _passtd::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(rough1.0f)*(rough1.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 r0pow((n1-n2)/(n1n2),2.0f);float _passtd::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 Snells 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(cosalpha0.0f){float F;Vector3f specular;specularcaculate(wo,wi,N,F);Vector3f diffuse1.0f/M_PI;float _kd1-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.747f0.058f, 0.747f0.258f, 0.747f) 15.6f * Vector3f(0.740f0.287f,0.740f0.160f,0.740f) 18.4f *Vector3f(0.737f0.642f,0.737f0.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.747f0.058f, 0.747f0.258f, 0.747f) 15.6f * Vector3f(0.740f0.287f,0.740f0.160f,0.740f) 18.4f *Vector3f(0.737f0.642f,0.737f0.159f,0.737f)));light-Kd Vector3f(0.65f);Material* whitem new Material(Microfacet,Vector3f(0.0f));whitem-KsVector3f(0.45f,0.45f,0.45f);whitem-KdVector3f(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_caststd::chrono::hours(stop - start).count() hours\n;std::cout : std::chrono::duration_caststd::chrono::minutes(stop - start).count() minutes\n;std::cout : std::chrono::duration_caststd::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 ScaleVector3f(1.0f,1.0f,1.0f),Vector3f moveVector3f(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_limitsfloat::infinity(),std::numeric_limitsfloat::infinity(),std::numeric_limitsfloat::infinity()};Vector3f max_vert Vector3f{-std::numeric_limitsfloat::infinity(),-std::numeric_limitsfloat::infinity(),-std::numeric_limitsfloat::infinity()};for (int i 0; i mesh.Vertices.size(); i 3) {std::arrayVector3f, 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);vertvert*Scalemove;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::vectorObject* ptrs;for (auto tri : triangles){ptrs.push_back(tri);area tri.area;}bvh new BVHAccel(ptrs);}
