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
@@ -25,6 +26,14 @@ FetchContent_GetProperties(lodepng)
 if(NOT lodepng_POPULATED)
     FetchContent_Populate(lodepng)
 endif()
+FetchContent_Declare(
+    xoshiro
+    GIT_REPOSITORY https://github.com/manuel5975p/xoshiro-cpp.git
+)
+FetchContent_GetProperties(xoshiro)
+if(NOT xoshiro_POPULATED)
+    FetchContent_Populate(xoshiro)
+endif()
 
 add_executable(chnuscht 
     main.cpp
@@ -34,7 +43,9 @@ add_executable(chnuscht
 )
 target_include_directories(chnuscht PUBLIC "${objloader_SOURCE_DIR}/include")
 target_include_directories(chnuscht PUBLIC "${lodepng_SOURCE_DIR}")
+target_include_directories(chnuscht PUBLIC "${xoshiro_SOURCE_DIR}")
 option(EMBREE_TUTORIALS OFF FORCE)
+option(EMBREE_ISPC_SUPPORT OFF FORCE)
 FetchContent_Declare(
     embree
     GIT_REPOSITORY https://github.com/embree/embree/
@@ -47,7 +58,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,8 @@ struct hit{
     Vector3f normal;
     Vector3f indir;
     Vector3f outdir;
+    const material& mat;
+    Vector2f uv;
 };
 std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sampler _rng)const{
     using namespace Eigen;
@@ -25,8 +27,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 < 256;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 +41,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,16 +55,21 @@ 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(Vector2f(0,0)).cast<double>());
                             for(auto it = hits.rbegin();it != hits.rend();it++){
-                                transport.array() *= it->outdir.dot(it->normal);
+                                transport.array() *= (it->mat.m_bsdf->eval(
+                                    -it->outdir,
+                                    -it->indir,
+                                    it->normal
+                                ).array().cast<double>()) * it->mat.tex->eval(it->uv).array().cast<double>();
                                 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();
+                                    sh_ray /= sh_ray_length;
                                     RTCRay rsh_ray;
                                     rsh_ray.org_x = it->pos.x();
                                     rsh_ray.org_y = it->pos.y();
@@ -77,7 +84,7 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                                     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.0 / double(sh_ray_length * sh_ray_length));
+                                        transport.array() += pl.color.cast<double>().array() * it->mat.m_bsdf->eval(-sh_ray.normalized(), it->indir, it->normal).cast<double>().array() * it->mat.tex->eval(it->uv).array().cast<double>() * (1.0 / double(sh_ray_length * sh_ray_length));
                                     }
                                 }
                             }
@@ -93,11 +100,31 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                         if(newdir.dot(hitnormal) <= 0){
                             newdir *= -1.0f;
                         }
+                        Vector2f uvcoordinate(0,0);
+                        {
+                            if(auto v = dynamic_cast<mesh_object*>(sc.added_objects[sc.geom_id_to_object_index_map.find(rayhit.hit.geomID)->second].get())){
+                                std::array<size_t, 3> triangle_indices;
+                                triangle_indices[0] = sc.added_meshes[v->mesh_index].Indices[rayhit.hit.primID * 3 + 0];
+                                triangle_indices[1] = sc.added_meshes[v->mesh_index].Indices[rayhit.hit.primID * 3 + 1];
+                                triangle_indices[2] = sc.added_meshes[v->mesh_index].Indices[rayhit.hit.primID * 3 + 2];
+                                const objl::Vertex& v0 = sc.added_meshes[v->mesh_index].Vertices[triangle_indices[0]];
+                                const objl::Vertex& v1 = sc.added_meshes[v->mesh_index].Vertices[triangle_indices[1]];
+                                const objl::Vertex& v2 = sc.added_meshes[v->mesh_index].Vertices[triangle_indices[2]];
+                                float hu = rayhit.hit.u;
+                                float hv = rayhit.hit.v;
+                                uvcoordinate = 
+                                hu  * Eigen::Map<const Eigen::Vector2f>(reinterpret_cast<const float*>(&v1.TextureCoordinate)) + 
+                                hv  * Eigen::Map<const Eigen::Vector2f>(reinterpret_cast<const float*>(&v2.TextureCoordinate)) + 
+                     (1 - hu  - hv) * Eigen::Map<const Eigen::Vector2f>(reinterpret_cast<const float*>(&v0.TextureCoordinate));
+                            }
+                        }
                         hits.push_back(hit{
                             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),
+                            uvcoordinate
                         });
                         transport *= newdir.dot(hitnormal);
                         rayhit.ray.dir_x = newdir.x();
@@ -107,7 +134,11 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                     else{
                         Eigen::Vector3d transport(0,0,0);
                         for(auto it = hits.rbegin();it != hits.rend();it++){
-                            transport.array() *= it->outdir.dot(it->normal);
+                            transport.array() *= it->mat.m_bsdf->eval(
+                                -it->outdir,
+                                -it->indir,
+                                it->normal
+                            ).array().cast<double>() * it->mat.tex->eval(it->uv).array().cast<double>();
                             for(const auto& pl : sc.pointlight){
                                 struct RTCIntersectContext sh_context;
                                 rtcInitIntersectContext(&sh_context);
@@ -129,7 +160,7 @@ std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sample
                                 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);
+                                    transport.array() += pl.color.cast<double>().array() * it->mat.m_bsdf->eval(-sh_ray.normalized(), it->indir, it->normal).cast<double>().array() * it->mat.tex->eval(it->uv).array().cast<double>() * (1.0 / double(sh_ray_length * sh_ray_length));
                                 }
                             }
                         }
@@ -137,10 +168,45 @@ 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->mat.m_bsdf->eval(
+                        -it->outdir,
+                        -it->indir,
+                        it->normal
+                    ).array().cast<double>()) * it->mat.tex->eval(it->uv).array().cast<double>();
+                    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.array() += pl.color.cast<double>().array() * it->mat.m_bsdf->eval(-sh_ray.normalized(), it->indir, it->normal).cast<double>().array() * it->mat.tex->eval(it->uv).array().cast<double>() * (1.0 / double(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>() / 128.0f).array().pow(0.5f).matrix();
         }
     }
     }

main.cpp

@@ -22,6 +22,8 @@ RTCDevice initializeDevice(){
     return device;
 }
 int main(){
+    using std::cos;
+    using std::sin;
     using namespace Eigen;
     RTCDevice dev = initializeDevice();
 
@@ -33,23 +35,37 @@ int main(){
     objl::Mesh ajax_mesh = ajax.LoadedMeshes[0];
     objl::Mesh quad_mesh = quad.LoadedMeshes[0];
     Matrix4f trf;
-    trf << 1,0,0,0,
-           0,1.5,0,0,
-           0,0,1,0,
-           0,0,0,1;
-    scene.add_object(ajax_mesh, trf, false);
+    std::uniform_real_distribution<float> angles(0, 2 * M_PI);
+    std::uniform_real_distribution<float> o1(0, 1);
+    xoshiro_256 gangles(42);
+    for(float ax = -10;ax < 10;ax += 5){
+        for(float ay = -10;ay < 10;ay += 5){
+            float angle = angles(gangles);
+            trf << cos(angle),0,-sin(angle),ax,
+                   0,1.6,0,0,
+                   sin(angle),0,cos(angle),ay,
+                   0,0,0,1;
+            scene.add_object(ajax_mesh, trf, false);
+            scene.added_objects.back()->mat.tex = std::make_unique<uni_texture>(Color(o1(gangles), o1(gangles), o1(gangles)));
+            //scene.added_objects.back()->mat.m_bsdf = std::make_unique<microfacet_bsdf>(0.8);
+        }
+    }
+    
     trf << 40,0,0,0,
-           0,1,0,-1.5,
-           0,0,4,0,
+           0,1,0,-1.6,
+           0,0,40,0,
            0,0,0,1;
     scene.add_object(quad_mesh, trf, false);
+    scene.added_objects.back()->mat.tex = std::make_unique<quad_checkerboard_texture>(Color(1,0,0), Color(0,1,0), Color(0,0,1), Color(0,0,0));
 
     //scene.add_sphere(Vector3f(-1.5,0,0),0.5, true);
-    scene.pointlight.push_back(pointlight{Eigen::Vector3f(2,0,0), Eigen::Vector3f(1,1,0)});
+    scene.pointlight.push_back(pointlight{Eigen::Vector3f(2,0,0), Eigen::Vector3f(1,0,0)});
+    scene.pointlight.push_back(pointlight{Eigen::Vector3f(2,0,1), Eigen::Vector3f(0,1,0)});
+    scene.pointlight.push_back(pointlight{Eigen::Vector3f(2,0,2), Eigen::Vector3f(0,0,1)});
     scene.build();
-    camera cam(Vector3f(0, 0.6, 2.4),Vector3f(0,0,0), Vector3f(0,1,0));
+    camera cam(Vector3f(0, 0.6 * 5, 2.4 * 5),Vector3f(0,0,0), Vector3f(0,1,0));
     sampler simp;
-    size_t n = 2000;
+    size_t n = 1000;
     auto t1 = nanoTime();
     auto x = cam.get_image(scene, n, simp);
     auto t2 = nanoTime();

material.hpp

@@ -0,0 +1,91 @@
+#ifndef TEXTURE_HPP
+#define TEXTURE_HPP
+#include <Eigen/Core>
+#include <memory>
+#include "warping.hpp"
+using Color = Eigen::Vector3f;
+struct texture{
+    virtual Eigen::Vector3f eval(const Eigen::Vector2f& uv)const{
+        return Eigen::Vector3f(1,1,1);
+    }
+};
+struct uni_texture : texture{
+    Color c;
+    uni_texture(const Color& _c) : c(_c){}
+    virtual Eigen::Vector3f eval(const Eigen::Vector2f& uv)const override{
+        return c;
+    }
+};
+struct checkerboard_texture : texture{
+    Color c1;
+    Color c2;
+    checkerboard_texture(const Color& _c1, const Color& _c2) : c1(_c1), c2(_c2){}
+    virtual Eigen::Vector3f eval(const Eigen::Vector2f& uv)const override{
+        int i = uv(0) * 100.f;
+        int j = uv(1) * 100.f;
+        return ((i ^ j) & 1) ? c1 : c2;
+    }
+};
+
+struct quad_checkerboard_texture : texture{
+    Color c1;
+    Color c2;
+    Color c3;
+    Color c4;
+    quad_checkerboard_texture(const Color& _c1, const Color& _c2, const Color& _c3, const Color& _c4) : c1(_c1), c2(_c2), c3(_c3), c4(_c4){}
+    virtual Eigen::Vector3f eval(const Eigen::Vector2f& uv)const override{
+        int i = uv(0) * 100.f;
+        int j = uv(1) * 100.f;
+        int id = ((i ^ j) & 3);
+        switch(id){
+            case 0:
+            return c1;
+            case 1:
+            return c2;
+            case 2:
+            return c3;
+            case 3:
+            default:
+            return c4;
+        }
+    }
+};
+struct bsdf{
+    virtual Color eval(const Eigen::Vector3f& in, const Eigen::Vector3f& out, const Eigen::Vector3f& normal)const{
+        throw 5;
+        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 x = 1.0f * std::exp(mirror_ref.dot(out) - 1);
+        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,21 @@
 #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)){}
+    virtual void kueder(){}
+};
 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)){}
+    virtual void kueder(){}
 };
 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)){}
+    virtual void kueder(){}
 };
 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