gaht eigentlich

.gitignore

@@ -0,0 +1,2 @@
+build/
+.vscode/

CMakeLists.txt

@@ -0,0 +1,53 @@
+cmake_minimum_required (VERSION 3.16)
+include(FetchContent)
+project(Flaschenray)
+if (NOT CMAKE_BUILD_TYPE)
+    set(CMAKE_BUILD_TYPE Release)
+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 "-DNDEBUG -Ofast -fopenmp -march=native -funroll-loops -std=c++20")
+FetchContent_Declare(
+    objloader
+    GIT_REPOSITORY https://github.com/manuel5975p/OBJ-Loader
+)
+FetchContent_GetProperties(objloader)
+if(NOT objloader_POPULATED)
+    FetchContent_Populate(objloader)
+    #add_subdirectory(${objloader_SOURCE_DIR} ${objloader_BINARY_DIR})
+endif()
+FetchContent_Declare(
+    lodepng
+    GIT_REPOSITORY https://github.com/lvandeve/lodepng.git
+)
+FetchContent_GetProperties(lodepng)
+if(NOT lodepng_POPULATED)
+    FetchContent_Populate(lodepng)
+endif()
+
+add_executable(chnuscht 
+    main.cpp
+    camera.cpp
+    scene.cpp
+    "${lodepng_SOURCE_DIR}/lodepng.cpp"
+)
+target_include_directories(chnuscht PUBLIC "${objloader_SOURCE_DIR}/include")
+target_include_directories(chnuscht PUBLIC "${lodepng_SOURCE_DIR}")
+option(EMBREE_TUTORIALS OFF FORCE)
+FetchContent_Declare(
+    embree
+    GIT_REPOSITORY https://github.com/embree/embree/
+)
+FetchContent_GetProperties(embree)
+if(NOT embree_POPULATED)
+    FetchContent_Populate(embree)
+    add_subdirectory(${embree_SOURCE_DIR} ${embree_BINARY_DIR})
+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")

camera.cpp

@@ -0,0 +1,85 @@
+#include "camera.hpp"
+#include <Eigen/Dense>
+#include "warping.hpp"
+std::vector<Eigen::Vector3f> camera::get_image(const scene& sc, size_t n, sampler _rng)const{
+    using namespace Eigen;
+    std::vector<Eigen::Vector3f> img(n * n, Eigen::Vector3f(0,0,0));
+    const Eigen::Vector3f lookray = (look_at - pos).normalized();
+	const Eigen::Vector3f left = lookray.cross(up).normalized();
+	const Eigen::Vector3f realup = left.cross(lookray);
+    const int n2 = n / 2;
+    size_t raycasts = 0;
+    #pragma omp parallel
+    {
+    sampler rng;
+    std::uniform_real_distribution<float> pert_dis(-0.5,0.5);
+    #pragma omp for collapse(2) schedule(guided) reduction(+:raycasts)
+    for(int j = 0;j < n;j++){
+		for(int i = 0;i < n;i++){
+            
+            Eigen::Vector3d accum(0,0,0);
+            
+            for(size_t smp = 0;smp < 4096;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;
+                rayhit.ray.org_x = pos.x();
+                rayhit.ray.org_y = pos.y();
+                rayhit.ray.org_z = pos.z();
+                rayhit.ray.dir_x = di.x();
+                rayhit.ray.dir_y = di.y();
+                rayhit.ray.dir_z = di.z();
+                Vector3f hitpos = pos;
+                Vector3f transport(1,1,1);
+                for(size_t bounces = 0;bounces < 5;bounces++){
+                    struct RTCIntersectContext context;
+                    rtcInitIntersectContext(&context);
+                    rayhit.ray.tnear = 0.00002f;
+                    rayhit.ray.tfar = std::numeric_limits<float>::infinity();
+                    rayhit.ray.mask = -1;
+                    rayhit.ray.flags = 0;
+                    rayhit.hit.geomID = RTC_INVALID_GEOMETRY_ID;
+                    rayhit.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
+                    rtcIntersect1(sc.m_scene, &context, &rayhit);
+                    raycasts++;
+                    if(rayhit.hit.geomID != RTC_INVALID_GEOMETRY_ID){
+                        if(sc.emitters.contains(rayhit.hit.geomID)){
+                            accum += (Eigen::Vector3d(0,1,1).array() * transport.cast<double>().array()).matrix();
+                            //std::cout << "sd" << std::endl;
+                            goto inner;
+                        }
+                        //std::cout << std::to_string(rayhit.ray.tfar) + "\n";
+                        rayhit.ray.org_x += rayhit.ray.dir_x * rayhit.ray.tfar;
+                        rayhit.ray.org_y += rayhit.ray.dir_y * rayhit.ray.tfar;
+                        rayhit.ray.org_z += rayhit.ray.dir_z * rayhit.ray.tfar;
+                        Vector3f hitnormal = Vector3f(rayhit.hit.Ng_x, rayhit.hit.Ng_y, rayhit.hit.Ng_z).normalized();
+                        //std::cout << hitnormal << "\n\n";
+                        //Vector3f newdir = cosine_hemisphere(hitnormal, rng.next2D()).normalized();
+                        Vector3f newdir = uniform_sphere(rng.next2D());
+                        if(newdir.dot(hitnormal) <= 0){
+                            newdir *= -1.0f;
+                        }
+                        transport *= newdir.dot(hitnormal);
+                        //std::cout << std::to_string(newdir.dot(Vector3f(rayhit.hit.Ng_x, rayhit.hit.Ng_y, rayhit.hit.Ng_z).normalized())) + std::string("\n");
+                        rayhit.ray.dir_x = newdir.x();
+                        rayhit.ray.dir_y = newdir.y();
+                        rayhit.ray.dir_z = newdir.z();
+                    }
+                    else{
+                        goto inner;
+                    }
+                }
+                inner:
+                (void)0;
+            }
+            img[j * n + i] = (accum.cast<float>() / 512.0f).array().pow(0.6f).matrix();
+        }
+    }
+    }
+    std::cout << std::to_string(raycasts) + " Raycasts\n";
+    return img;
+}
+
+camera::camera(const Eigen::Vector3f& loc, const Eigen::Vector3f& look, const Eigen::Vector3f& u) : pos(loc), look_at(look), up(u){
+	
+}

camera.hpp

@@ -0,0 +1,16 @@
+#ifndef E_CAMERA_HPP
+#define E_CAMERA_HPP
+#include <Eigen/Core>
+#include <embree3/rtcore.h>
+#include <xoshiro.hpp>
+#include "sampler.hpp"
+#include "scene.hpp"
+struct scene;
+struct camera{
+    Eigen::Vector3f pos;
+    Eigen::Vector3f look_at;
+    Eigen::Vector3f up;
+    std::vector<Eigen::Vector3f> get_image(const scene& sc, size_t n, sampler rng)const;
+    camera(const Eigen::Vector3f& loc, const Eigen::Vector3f& look, const Eigen::Vector3f& u);
+};
+#endif

main.cpp

@@ -0,0 +1,69 @@
+#include <embree3/rtcore.h>
+#include <cstdio>
+#include <limits>
+#include <iostream>
+#include <OBJ_Loader.h>
+#include <bmp.hpp>
+#include "scene.hpp"
+#include "camera.hpp"
+#include <time.hpp>
+#include <lodepng.h>
+#include <algorithm>
+#include <Eigen/Geometry>
+using std::printf;
+void errorFunction(void* userPtr, enum RTCError error, const char* str){
+    printf("error %d: %s\n", error, str);
+}
+RTCDevice initializeDevice(){
+    RTCDevice device = rtcNewDevice("threads=24,user_threads=24");
+    if (!device)
+        printf("error %d: cannot create device\n", rtcGetDeviceError(NULL));
+    rtcSetDeviceErrorFunction(device, errorFunction, NULL);
+    return device;
+}
+int main(){
+    using namespace Eigen;
+    RTCDevice dev = initializeDevice();
+
+    scene scene(dev);
+    objl::Loader ajax;
+    objl::Loader quad;
+    ajax.LoadFile("ajax.obj");
+    quad.LoadFile("quad.obj");
+    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);
+    trf << 40,0,0,0,
+           0,1,0,-1.5,
+           0,0,4,0,
+           0,0,0,1;
+    scene.add_object(quad_mesh, trf, false);
+
+    scene.add_sphere(Vector3f(-1.5,0,0),0.5, true);
+    scene.build();
+    camera cam(Vector3f(0, 0.6 * 3, 2.4 * 3),Vector3f(0,0,0), Vector3f(0,1,0));
+    sampler simp;
+    size_t n = 2000;
+    auto t1 = nanoTime();
+    auto x = cam.get_image(scene, n, simp);
+    auto t2 = nanoTime();
+    std::cout << (t2 - t1) / 1000 / 1000.0 << " miliseconds\n";
+    std::vector<unsigned char> lodepng_image(x.size() * 4);
+    for(size_t i = 0;i < x.size();i++){
+        lodepng_image[i * 4] =     (unsigned char)(std::min(x[x.size() - i - 1].x() * 256,255.f));
+        lodepng_image[i * 4 + 1] = (unsigned char)(std::min(x[x.size() - i - 1].y() * 256,255.f));
+        lodepng_image[i * 4 + 2] = (unsigned char)(std::min(x[x.size() - i - 1].z() * 256,255.f));
+        lodepng_image[i * 4 + 3] = 255;
+    }
+    lodepng::encode("bild.png", lodepng_image, n, n);
+    //drawbmp("bild.bmp", x, n, n);
+    
+    //std::cout << (rayhit.ray.tfar) << "\n";
+    //std::cout << (rayhit.hit.Ng_x) << "\n";
+    return 0;
+}

sampler.hpp

@@ -0,0 +1,22 @@
+#ifndef E_SAMPLER_HPP
+#define E_SAMPLER_HPP
+#include <xoshiro.hpp>
+#include <Eigen/Core>
+template<typename RNG>
+struct sampler_tmplt{
+    RNG m_rng;
+    std::uniform_real_distribution<float> m_dis;
+    sampler_tmplt() : m_rng(), m_dis(0,1){};
+    sampler_tmplt(uint64_t seed) : m_rng(seed), m_dis(0,1){};
+    float next1D(){
+        return m_dis(m_rng);
+    }
+    Eigen::Vector2f next2D(){
+        return Eigen::Vector2f(m_dis(m_rng), m_dis(m_rng));
+    }
+    Eigen::Vector3f next3D(){
+        return Eigen::Vector3f(m_dis(m_rng), m_dis(m_rng), m_dis(m_rng));
+    }
+};
+using sampler = sampler_tmplt<xoshiro_256>;
+#endif

scene.cpp

@@ -0,0 +1,183 @@
+#include "scene.hpp"
+#include <utility>
+#include <algorithm>
+using namespace Eigen;
+scene::scene(RTCDevice dev){
+    m_scene = rtcNewScene(dev);
+    m_dev = dev;
+}
+unsigned int scene::add_object(const objl::Mesh& mesh, const Eigen::Matrix4f& transform, bool emitter){
+    added_meshes.push_back(mesh);
+    return add_object(added_meshes.size() - 1, transform, emitter);
+}
+unsigned int scene::add_object(size_t mesh_index, const Eigen::Matrix4f& transform, bool emitter){
+    using namespace objl;
+    const objl::Mesh& mesh = this->added_meshes[mesh_index];
+    std::pair<float, float> xex = {mesh.Vertices[0].Position.X, mesh.Vertices[0].Position.X};
+    std::pair<float, float> yex = {mesh.Vertices[0].Position.Y, mesh.Vertices[0].Position.Y};
+    std::pair<float, float> zex = {mesh.Vertices[0].Position.Z, mesh.Vertices[0].Position.Z};
+    for(size_t i = 1;i < mesh.Vertices.size();i++){
+        xex.first = std::min(mesh.Vertices[i].Position.X, xex.first);
+        yex.first = std::min(mesh.Vertices[i].Position.Y, yex.first);
+        zex.first = std::min(mesh.Vertices[i].Position.Z, zex.first);
+
+        xex.second = std::max(mesh.Vertices[i].Position.X, xex.second);
+        yex.second = std::max(mesh.Vertices[i].Position.Y, yex.second);
+        zex.second = std::max(mesh.Vertices[i].Position.Z, zex.second);
+    }
+    if(xex.second - xex.first < 0.0001f){
+        xex.second = 1e30;
+        xex.first = -1e30;
+    }
+    if(yex.second - yex.first < 0.0001f){
+        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;
+    }
+    std::vector<objl::Vertex> rectified_vertices(mesh.Vertices.size());
+    for(size_t i = 0;i < mesh.Vertices.size();i++){
+        Vector4f rec_pos;
+        rec_pos.x() = (mesh.Vertices[i].Position.X - xex.first) * (2.0f / (xex.second - xex.first)) - 1.0f;
+        rec_pos.y() = (mesh.Vertices[i].Position.Y - yex.first) * (2.0f / (yex.second - yex.first)) - 1.0f;
+        rec_pos.z() = (mesh.Vertices[i].Position.Z - zex.first) * (2.0f / (zex.second - zex.first)) - 1.0f;
+        rec_pos.w() = 1;
+        rec_pos = transform * rec_pos;
+        rectified_vertices[i].Position.X = rec_pos.x();
+        rectified_vertices[i].Position.Y = rec_pos.y();
+        rectified_vertices[i].Position.Z = rec_pos.z();
+    }
+    RTCGeometry m_geom = rtcNewGeometry(m_dev, RTC_GEOMETRY_TYPE_TRIANGLE);
+
+            
+            float* vertices = (float*) rtcSetNewGeometryBuffer(m_geom,
+                                                     RTC_BUFFER_TYPE_VERTEX,
+                                                     0,
+                                                     RTC_FORMAT_FLOAT3,
+                                                     3 * sizeof(float),
+                                                     rectified_vertices.size());
+            unsigned* indices = (unsigned*) rtcSetNewGeometryBuffer(m_geom,
+                                                          RTC_BUFFER_TYPE_INDEX,
+                                                          0,
+                                                          RTC_FORMAT_UINT3,
+                                                          3 * sizeof(unsigned),
+                                                          mesh.Indices.size() / 3);
+    if (vertices && indices){
+        for(size_t i = 0;i < mesh.Indices.size();i++){
+            indices[i] = mesh.Indices[i];
+        }
+        for(size_t i = 0;i < rectified_vertices.size();i++){
+            vertices[3 * i] =     rectified_vertices[i].Position.X;
+            vertices[3 * i + 1] = rectified_vertices[i].Position.Y;
+            vertices[3 * i + 2] = rectified_vertices[i].Position.Z;
+        }
+    }
+    rtcUpdateGeometryBuffer(m_geom, RTC_BUFFER_TYPE_VERTEX, 0);
+    rtcUpdateGeometryBuffer(m_geom, RTC_BUFFER_TYPE_INDEX, 0);
+    rtcCommitGeometry(m_geom);
+    unsigned int geomID = rtcAttachGeometry(m_scene, m_geom);
+    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->geom_id_to_object_index_map[geomID] = this->added_objects.size() - 1;
+    return geomID;
+}
+unsigned int scene::add_sphere(const Eigen::Vector3f& pos, float rad, bool emitter){
+    RTCGeometry m_geom = rtcNewGeometry(m_dev, RTC_GEOMETRY_TYPE_SPHERE_POINT);
+    float* vertices = (float*) rtcSetNewGeometryBuffer(m_geom,
+                                                     RTC_BUFFER_TYPE_VERTEX,
+                                                     0,
+                                                     RTC_FORMAT_FLOAT4,
+                                                     4 * sizeof(float),
+                                                     1);
+    vertices[0] = pos.x();
+    vertices[1] = pos.y();
+    vertices[2] = pos.z();
+    vertices[3] = rad;
+    rtcUpdateGeometryBuffer(m_geom, RTC_BUFFER_TYPE_VERTEX, 0);
+    rtcCommitGeometry(m_geom);
+    unsigned int geomID = rtcAttachGeometry(m_scene, m_geom);
+    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->geom_id_to_object_index_map[geomID] = this->added_objects.size() - 1;
+    return geomID;
+}
+void scene::build(){
+    rtcSetSceneBuildQuality(m_scene, RTC_BUILD_QUALITY_HIGH);
+    rtcCommitScene(m_scene);
+}
+/*scene get_scene(RTCDevice dev){
+    RTCScene rtcsc = rtcNewScene(dev);
+    objl::Loader loader;
+    loader.LoadFile("ajax.obj");
+    objl::Mesh& ajaxmesh = loader.LoadedMeshes.at(0);
+    //std::cout << ajaxmesh.Vertices.size() << "\n";
+    if(ajaxmesh.Indices.size() % 3 != 0){
+        std::terminate();
+    }
+    RTCGeometry m_geom = rtcNewGeometry(dev, RTC_GEOMETRY_TYPE_TRIANGLE);
+
+            
+            float* vertices = (float*) rtcSetNewGeometryBuffer(m_geom,
+                                                     RTC_BUFFER_TYPE_VERTEX,
+                                                     0,
+                                                     RTC_FORMAT_FLOAT3,
+                                                     3 * sizeof(float),
+                                                     ajaxmesh.Vertices.size());
+            unsigned* indices = (unsigned*) rtcSetNewGeometryBuffer(m_geom,
+                                                          RTC_BUFFER_TYPE_INDEX,
+                                                          0,
+                                                          RTC_FORMAT_UINT3,
+                                                          3 * sizeof(unsigned),
+                                                          ajaxmesh.Indices.size() / 3);
+    if (vertices && indices){
+        for(size_t i = 0;i < ajaxmesh.Indices.size();i++){
+            indices[i] = ajaxmesh.Indices[i];
+        }
+        for(size_t i = 0;i < ajaxmesh.Vertices.size();i++){
+            vertices[3 * i] =     ajaxmesh.Vertices[i].Position.X;
+            vertices[3 * i + 1] = ajaxmesh.Vertices[i].Position.Y;
+            vertices[3 * i + 2] = ajaxmesh.Vertices[i].Position.Z;
+        }
+    }
+    rtcUpdateGeometryBuffer(m_geom, RTC_BUFFER_TYPE_VERTEX, 0);
+    rtcUpdateGeometryBuffer(m_geom, RTC_BUFFER_TYPE_INDEX, 0);
+    rtcCommitGeometry(m_geom);
+    unsigned int geomID = rtcAttachGeometry(rtcsc, m_geom);
+    RTCGeometry m_geom2 = rtcNewGeometry(dev, RTC_GEOMETRY_TYPE_TRIANGLE);
+    vertices = (float*) rtcSetNewGeometryBuffer(m_geom2,
+                                                     RTC_BUFFER_TYPE_VERTEX,
+                                                     0,
+                                                     RTC_FORMAT_FLOAT3,
+                                                     3 * sizeof(float),
+                                                     3);
+    indices = (unsigned*) rtcSetNewGeometryBuffer(m_geom2,
+                                                          RTC_BUFFER_TYPE_INDEX,
+                                                          0,
+                                                          RTC_FORMAT_UINT3,
+                                                          3 * sizeof(unsigned),
+                                                          1);
+    
+    vertices[0] = -11; vertices[1] = 10; vertices[2] = 5;
+    vertices[3] = -11; vertices[4] = 15; vertices[5] = 5;
+    vertices[6] = -20; vertices[7] = 15; vertices[8] = -5;
+    indices[0] = 0;
+    indices[1] = 1;
+    indices[2] = 2;
+    rtcUpdateGeometryBuffer(m_geom2, RTC_BUFFER_TYPE_VERTEX, 0);
+    rtcUpdateGeometryBuffer(m_geom2, RTC_BUFFER_TYPE_INDEX, 0);
+    rtcCommitGeometry(m_geom2);
+    unsigned int geid2 = rtcAttachGeometry(rtcsc, m_geom2);
+    rtcSetSceneBuildQuality(rtcsc, RTC_BUILD_QUALITY_HIGH);
+    rtcCommitScene(rtcsc);
+    scene sc(rtcsc);
+    //sc.emitting.insert({geomID, true});
+    sc.emitting.insert({geid2, true});
+    return sc;
+}*/

scene.hpp

@@ -0,0 +1,40 @@
+#ifndef SCENE_HPP
+#define SCENE_HPP
+#include <embree3/rtcore.h>
+#include <unordered_set>
+#include <unordered_map>
+#include <OBJ_Loader.h>
+#include <Eigen/Dense>
+#include <Eigen/Geometry>
+#include <memory>
+#include "texture.hpp"
+struct scene_object {};
+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)){}
+};
+struct sphere_object : scene_object{
+    std::unique_ptr<texture> tex;
+    sphere_object(std::unique_ptr<texture>&& _tex) : tex(std::move(_tex)){}
+};
+struct pointlight{
+    Eigen::Vector3f pos;
+    Eigen::Vector3f color;
+};
+struct scene{
+    RTCScene m_scene;
+    RTCDevice m_dev;
+    std::unordered_set<unsigned int> emitters;
+    std::unordered_map<unsigned int, size_t> geom_id_to_object_index_map;
+    std::vector<objl::Mesh> added_meshes;
+    std::vector<pointlight> pointlight;
+    std::vector<std::unique_ptr<scene_object>> added_objects;
+    scene(RTCDevice sc);
+    unsigned int add_object(const objl::Mesh& mesh, const Eigen::Matrix4f& transform, bool emitter = false);
+    unsigned int add_object(size_t mesh_index, const Eigen::Matrix4f& transform, bool emitter = false);
+    unsigned int add_sphere(const Eigen::Vector3f& pos, float rad, bool emitter = false);
+    void build();
+};
+scene get_scene(RTCDevice dev);
+#endif

texture.hpp

@@ -0,0 +1,17 @@
+#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

@@ -0,0 +1,54 @@
+#ifndef WARPING_HPP
+#define WARPING_HPP
+#include <Eigen/Dense>
+#include <cmath>
+#include <random>
+inline Eigen::Vector3f uniform_sphere(const Eigen::Vector2f& sample){
+    using std::acos;
+    using std::cos;
+    using std::sin;
+    std::uniform_real_distribution<float> dis(0, 1);
+    float theta = 2 * M_PI * sample(0);
+    float phi = acos(2 * sample(1) - 1);
+    Eigen::Vector3f ret;
+    ret(0) = sin(phi) * cos(theta);
+    ret(1) = sin(phi) * sin(theta);
+    ret(2) = cos(phi);
+    return ret;
+}
+inline Eigen::Vector3f cosine_hemisphere(const Eigen::Vector3f& normal, const Eigen::Vector2f& sample) {
+    using std::asin;
+    using std::cos;
+    using std::sin;
+    using std::sqrt;
+	//constexpr float r = 1;
+	//float theta = asin(sqrt(sample.y()));
+	//float phi = 2 * M_PI * sample.x();
+	//float z = r * cos(theta);
+	//float _r = sqrt(r * r - z * z);
+	//float x = _r * cos(phi);
+	//float y = _r * sin(phi);
+    float r = sqrt(sample.x());
+    float phi = 2 * M_PI * sample.y();
+    float x = r * cos(phi);
+    float y = r * sin(phi);
+    float z = sqrt(1 - r * r);
+    Eigen::Vector3f notnormal = normal;
+    if(normal.y() == 0){
+        std::swap(notnormal.x(), notnormal.z());
+        notnormal.z() *= -1;
+    }
+    else{
+        std::swap(notnormal.x(), notnormal.y());
+        notnormal.y() *= -1;
+    }
+    Eigen::Vector3f notnormal2 = normal.cross(notnormal);
+    notnormal = notnormal2.cross(normal);
+    Eigen::Matrix3f trf;
+    trf.col(0) = notnormal;
+    trf.col(1) = notnormal2;
+    trf.col(2) = normal;
+    Eigen::Vector3f ret(x, y, z);
+	return Eigen::Vector3f(trf * ret);
+}
+#endif