Lambertian working

CMakeLists.txt

@@ -7,6 +7,7 @@ endif()
 #set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_EXPORT_COMPILE_COMMANDS True)
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -pthread -fopenmp")
+#set(CMAKE_CXX_FLAGS_RELEASE "-D_GLIBCXX_DEBUG -g3 -fopenmp -march=native -funroll-loops -std=c++20")
 set(CMAKE_CXX_FLAGS_RELEASE "-DNDEBUG -O3 -fopenmp -march=native -funroll-loops -std=c++20")
 FetchContent_Declare(
     objloader
@@ -47,7 +48,7 @@ endif()
 
 target_link_libraries(chnuscht embree)
 target_link_libraries(chnuscht tbb)
-target_link_libraries(chnuscht "gomp")
-target_link_libraries(chnuscht "sfml-system")
-target_link_libraries(chnuscht "sfml-graphics")
-target_link_libraries(chnuscht "sfml-window")
+#target_link_libraries(chnuscht "gomp")
+#target_link_libraries(chnuscht "sfml-system")
+#target_link_libraries(chnuscht "sfml-graphics")
+#target_link_libraries(chnuscht "sfml-window")

camera.cpp

@@ -7,6 +7,7 @@ struct hit{
     Vector3f normal;
     Vector3f indir;
     Vector3f outdir;
+    const material& mat;
 };
 std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sampler _rng)const{
     using namespace Eigen;
@@ -25,8 +26,9 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
 		for(int i = 0;i < n;i++){
             
             Eigen::Vector3d accum(0,0,0);
-            
-            for(size_t smp = 0;smp < 128;smp++){
+            std::vector<hit> hits;
+            hits.reserve(6);
+            for(size_t smp = 0;smp < 512;smp++){
                 Eigen::Vector3f di = lookray + (left * ((pert_dis(rng.m_rng) + float(i - n2)) / n2 / 1.5)) + realup * ((pert_dis(rng.m_rng) + float(j - n2)) / n2 / 1.5);
                 di.normalize();
                 struct RTCRayHit rayhit;
@@ -38,8 +40,7 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                 rayhit.ray.dir_z = di.z();
                 Vector3f hitpos = pos;
                 Vector3f transport(1,1,1);
-                std::vector<hit> hits;
-                hits.reserve(6);
+                hits.clear();
                 for(size_t bounces = 0;bounces < 6;bounces++){
                     struct RTCIntersectContext context;
                     rtcInitIntersectContext(&context);
@@ -53,7 +54,7 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                     raycasts++;
                     if(rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID){
                         if(sc.emitters.contains(rayhit.hit.geomID)){
-                            Eigen::Vector3d transport(0,1,1);
+                            Eigen::Vector3d transport(sc.added_objects[sc.geom_id_to_object_index_map.find(rayhit.hit.geomID)->second]->mat.tex->eval(0,0).cast<double>());
                             for(auto it = hits.rbegin();it != hits.rend();it++){
                                 transport.array() *= it->outdir.dot(it->normal);
                                 for(const auto& pl : sc.pointlight){
@@ -62,7 +63,7 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                                     Vector3f sh_ray(pl.pos - it->pos);
                                     if(sh_ray.dot(it->normal) < 0)continue;
                                     float sh_ray_length = sh_ray.norm();
-                                    sh_ray.normalize();
+                                    sh_ray /= sh_ray_length;
                                     RTCRay rsh_ray;
                                     rsh_ray.org_x = it->pos.x();
                                     rsh_ray.org_y = it->pos.y();
@@ -97,7 +98,8 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                             Vector3f(rayhit.ray.org_x, rayhit.ray.org_y, rayhit.ray.org_z),
                             hitnormal,
                             Vector3f(rayhit.ray.dir_x, rayhit.ray.dir_y, rayhit.ray.dir_z),
-                            newdir
+                            newdir,
+                            (sc.added_objects[sc.geom_id_to_object_index_map.find(rayhit.hit.geomID)->second]->mat)
                         });
                         transport *= newdir.dot(hitnormal);
                         rayhit.ray.dir_x = newdir.x();
@@ -137,10 +139,41 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                         goto inner;
                     }
                 }
+                {
+                Eigen::Vector3d transport(0,0,0);
+                for(auto it = hits.rbegin();it != hits.rend();it++){
+                    transport.array() *= it->outdir.dot(it->normal);
+                    for(const auto& pl : sc.pointlight){
+                        struct RTCIntersectContext sh_context;
+                        rtcInitIntersectContext(&sh_context);
+                        Vector3f sh_ray(pl.pos - it->pos);
+                        if(sh_ray.dot(it->normal) < 0)continue;
+                        float sh_ray_length = sh_ray.norm();
+                        sh_ray.normalize();
+                        RTCRay rsh_ray;
+                        rsh_ray.org_x = it->pos.x();
+                        rsh_ray.org_y = it->pos.y();
+                        rsh_ray.org_z = it->pos.z();
+                        rsh_ray.dir_x = sh_ray.x();
+                        rsh_ray.dir_y = sh_ray.y();
+                        rsh_ray.dir_z = sh_ray.z();
+                        rsh_ray.tnear = 0.001f;
+                        rsh_ray.tfar = sh_ray_length;
+                        rsh_ray.mask = -1;
+                        rsh_ray.flags = 0;
+                        raycasts++;
+                        rtcOccluded1(sc.m_scene, &sh_context, &rsh_ray);
+                        if(rsh_ray.tfar > 0){
+                            transport += pl.color.cast<double>() * std::abs(sh_ray.dot(it->normal)) * 1 / (sh_ray_length * sh_ray_length);
+                        }
+                    }
+                }
+                accum += transport;
+                }
                 inner:
                 (void)0;
             }
-            img[j * n + i] = (accum.cast<float>() / 32.0f).array().pow(0.5f).matrix();
+            img[j * n + i] = (accum.cast<float>() / 512.0f).array().pow(0.5f).matrix();
         }
     }
     }

main.cpp

@@ -37,7 +37,8 @@ int main(){
            0,1.5,0,0,
            0,0,1,0,
            0,0,0,1;
-    scene.add_object(ajax_mesh, trf, false);
+    unsigned aji = scene.add_object(ajax_mesh, trf, false);
+    scene.added_objects[scene.geom_id_to_object_index_map.find(aji)->second]->mat.m_bsdf = std::make_unique<microfacet_bsdf>(0.1);
     trf << 40,0,0,0,
            0,1,0,-1.5,
            0,0,4,0,
@@ -49,7 +50,7 @@ int main(){
     scene.build();
     camera cam(Vector3f(0, 0.6, 2.4),Vector3f(0,0,0), Vector3f(0,1,0));
     sampler simp;
-    size_t n = 2000;
+    size_t n = 500;
     auto t1 = nanoTime();
     auto x = cam.get_image(scene, n, simp);
     auto t2 = nanoTime();

material.hpp

@@ -0,0 +1,56 @@
+#ifndef TEXTURE_HPP
+#define TEXTURE_HPP
+#include <Eigen/Core>
+#include <memory>
+#include "warping.hpp"
+using Color = Eigen::Vector3f;
+struct texture{
+    virtual Eigen::Vector3f eval(float u, float v)const{
+        return Eigen::Vector3f(1,1,1);
+    }
+};
+struct uni_texture : texture{
+    Color c;
+    uni_texture(const Color& _c) : c(_c){}
+    virtual Eigen::Vector3f eval(float u, float v)const override{
+        return c;
+    }
+};
+struct bsdf{
+    virtual Color eval(const Eigen::Vector3f& in, const Eigen::Vector3f& out, const Eigen::Vector3f& normal)const{
+        return Eigen::Vector3f(1,1,1);
+    }
+};
+struct lambertian_bsdf : bsdf{
+    virtual Color eval(const Eigen::Vector3f& in, const Eigen::Vector3f& out, const Eigen::Vector3f& normal)const override{
+        float x = in.dot(-normal);
+        return Eigen::Vector3f(x,x,x);
+    }
+};
+struct microfacet_bsdf : bsdf{
+    const float alpha;
+    microfacet_bsdf() : alpha(0.1f){}
+    microfacet_bsdf(float _alpha) : alpha(_alpha){}
+    virtual Color eval(const Eigen::Vector3f& in, const Eigen::Vector3f& out, const Eigen::Vector3f& normal)const override{
+        float x = in.dot(-normal);
+        return Eigen::Vector3f(x,x,x);
+        Eigen::Vector3f mirror_ref = in + normal * (in.dot(normal));
+        //float x = beckmann_eval(mirror_ref, alpha);
+        float xd = 3.0f * std::exp(-mirror_ref.dot(out));
+        return Color(x,x,x);
+    }
+};
+struct material{
+    std::unique_ptr<texture> tex;
+    std::unique_ptr<bsdf> m_bsdf;
+    material() : material(Color(1,1,1)){
+
+    }
+    material(const Color& c, std::unique_ptr<bsdf>&& bsdf) : tex(std::make_unique<uni_texture>(c)), m_bsdf(std::move(bsdf)){
+        
+    }
+    material(const Color& c) : tex(std::make_unique<uni_texture>(c)), m_bsdf(std::make_unique<lambertian_bsdf>()){
+        
+    }
+};
+#endif

scene.cpp

@@ -33,7 +33,6 @@ unsigned int scene::add_object(size_t mesh_index, const Eigen::Matrix4f& transfo
         yex.second = 1e30;
         yex.first = -1e30;
     }
-    //std::cout << yex.second - yex.first << "\n";
     if(zex.second - zex.first < 0.0001f){
         zex.second = 1e30;
         zex.first = -1e30;
@@ -82,7 +81,7 @@ unsigned int scene::add_object(size_t mesh_index, const Eigen::Matrix4f& transfo
     if(emitter){
         emitters.emplace(geomID);
     }
-    this->added_objects.push_back(std::make_unique<mesh_object>(mesh_index, std::make_unique<uni_texture>(Color(1,0.5f,1))));
+    this->added_objects.push_back(std::make_unique<mesh_object>(mesh_index, material()));
     this->geom_id_to_object_index_map[geomID] = this->added_objects.size() - 1;
     return geomID;
 }
@@ -104,7 +103,7 @@ unsigned int scene::add_sphere(const Eigen::Vector3f& pos, float rad, bool emitt
     if(emitter){
         emitters.emplace(geomID);
     }
-    this->added_objects.push_back(std::make_unique<sphere_object>(std::make_unique<uni_texture>(Color(1,0.5f,1))));
+    this->added_objects.push_back(std::make_unique<sphere_object>(material()));
     this->geom_id_to_object_index_map[geomID] = this->added_objects.size() - 1;
     return geomID;
 }

scene.hpp

@@ -7,16 +7,18 @@
 #include <Eigen/Dense>
 #include <Eigen/Geometry>
 #include <memory>
-#include "texture.hpp"
-struct scene_object {};
+#include "material.hpp"
+
+struct scene_object {
+    material mat;
+    scene_object(material&& _mat) : mat(std::move(_mat)){}
+};
 struct mesh_object : scene_object{
     size_t mesh_index;
-    std::unique_ptr<texture> tex;
-    mesh_object(size_t mi, std::unique_ptr<texture>&& _tex) : mesh_index(mi), tex(std::move(_tex)){}
+    mesh_object(size_t mi, material&& _mat) : mesh_index(mi), scene_object(std::move(_mat)){}
 };
 struct sphere_object : scene_object{
-    std::unique_ptr<texture> tex;
-    sphere_object(std::unique_ptr<texture>&& _tex) : tex(std::move(_tex)){}
+    sphere_object(material&& _mat) : scene_object(std::move(_mat)){}
 };
 struct pointlight{
     Eigen::Vector3f pos;

texture.hpp

@@ -1,17 +0,0 @@
-#ifndef TEXTURE_HPP
-#define TEXTURE_HPP
-#include <Eigen/Core>
-using Color = Eigen::Vector3f;
-struct texture{
-    virtual Eigen::Vector3f eval(float u, float v)const{
-        return Eigen::Vector3f(1,1,1);
-    }
-};
-struct uni_texture : texture{
-    Color c;
-    uni_texture(const Color& _c) : c(_c){}
-    virtual Eigen::Vector3f eval(float u, float v)const override{
-        return c;
-    }
-};
-#endif

warping.hpp

@@ -51,4 +51,37 @@ inline Eigen::Vector3f cosine_hemisphere(const Eigen::Vector3f& normal, const Ei
     Eigen::Vector3f ret(x, y, z);
 	return Eigen::Vector3f(trf * ret);
 }
+inline float beckmann_eval(const Eigen::Vector3f& m, float alpha) {
+    using std::abs;
+    using std::sqrt;
+    using std::exp;
+    using std::pow;
+
+	if (m.z() > 0 && abs(m.squaredNorm() - 1) < 0.001f) {
+		float costheta = m.z();
+		float r = sqrt(m.x() * m.x() + m.y() * m.y());
+		float tantheta = r / m.z();
+		float upper = exp((-tantheta * tantheta) / (alpha * alpha));
+		float lower = M_PI * alpha * alpha * pow(costheta, 3);
+		return upper / lower;
+	}
+	return 0.0f;
+}
+
+inline Eigen::Vector3f beckmann(const Eigen::Vector2f& sample, float alpha) {
+	using std::cos;
+	using std::sin;
+	using std::sqrt;
+    using std::atan;
+    using std::log;
+
+	float theta = atan(sqrt(-alpha * alpha * log(1 - sample.y())));
+	float pha = 2 * M_PI * sample.x();
+	float z = cos(theta);
+	float r = sqrt(1 - z * z);
+	float y = r * cos(pha);
+	float x = r * sin(pha);
+	return Eigen::Vector3f(x, y, z);
+}
+
 #endif