diff --git a/CMakeLists.txt b/CMakeLists.txt index f0e8617..38d7e48 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ add_subdirectory(task02) add_subdirectory(task03) add_subdirectory(task04) add_subdirectory(task05) -# add_subdirectory(task06) +add_subdirectory(task06) # add_subdirectory(task07) # add_subdirectory(task08) # add_subdirectory(task09) diff --git a/README.md b/README.md index 926ac6c..51f45b3 100644 --- a/README.md +++ b/README.md @@ -50,8 +50,8 @@ Topics: |(3)
Apr. 22| **Parametric representation**
Bézier curve, polynominal | [task02](task02) | [[7] ](http://nobuyuki-umetani.com/acg2024s/parametric_curve.pdf) [[8]](http://nobuyuki-umetani.com/acg2024s/polynominal.pdf) | |(5)
May 7| **Coordinate transformation**
Affine, homography transformation | [task03](task03) | [[9]](http://nobuyuki-umetani.com/acg2024s/transformation.pdf) [[10]](http://nobuyuki-umetani.com/acg2024s/transformation_homogeneous_2d.pdf) [[11]](http://nobuyuki-umetani.com/acg2024s/transformation_homogeneous_3d.pdf) | |(4)
May 13| **Graphics pipeline**
depth buffer method, shading, shadow, anti aliasing | [task04](task04) | [[12]](http://nobuyuki-umetani.com/acg2024s/rasterization_3d.pdf) [[13]](http://nobuyuki-umetani.com/acg2024s/graphics_pipeline.pdf) | -|(6)
May 20| **Ray Casting 1**
spatial data structure | [task05](task05) | [[14]](http://nobuyuki-umetani.com/acg2024s/shading.pdf) [[15]](http://nobuyuki-umetani.com/acg2024s/rasterization_subpixel.pdf) [[16]](http://nobuyuki-umetani.com/acg2024s/implicit_modeling.pdf) | -|(7)
May 27| **Ray Casting 2**
Rendering equation, Monte Carlo integration | task06 | | +|(6)
May 20| **Ray Casting 1**
spatial data structure | [task05](task05) | [[14]](http://nobuyuki-umetani.com/acg2024s/shading.pdf) [[16]](http://nobuyuki-umetani.com/acg2024s/implicit_modeling.pdf) | +|(7)
May 27| **Ray Casting 2**
Rendering equation, Monte Carlo integration | [task06](task06) | [[15]](http://nobuyuki-umetani.com/acg2024s/rasterization_subpixel.pdf) [[17]](http://nobuyuki-umetani.com/acg2024s/ray_casting.pdf) [[18]](http://nobuyuki-umetani.com/acg2024s/monte_carlo_integration.pdf) [[19]](http://nobuyuki-umetani.com/acg2024s/ray_triangle_collision.pdf) | |(8)
June 3| **Character animation**
Linear blend skinning | task07 | | |(9)
June 10| Guest lecture by Dr. Rex West | | | |(10)
June 17| **Optimization**
Inverse kinematic | task08 | | @@ -83,8 +83,8 @@ Look at the following document. | [task03](task03) | **Perspectively-correct texture mapping**
rasterization of triangle, barycentric coordinate | | | [task04](task04) | **Vertex shader practice**
Rendering pipeline, mirror reflection, OpenGL | | | [task05](task05) | **Fragment shader practice**
Ray marching method, CSG modeling, implicit modeling | | -| task06 | **Monte Carlo integration**
Multiple importance sampling, path tracing, rendering equation | | -| task07 | TBD | | +| [task06](task06) | **Monte Carlo integration1**
Ambient occlusion, importance sampling, BVH | | +| task07 | **Monte Carlo integration2**
Multiple importance sampling | | | task08 | **Skeletal Character Animation**
Linear blend skinning, articulated rigid body | | | task09 | TBD | | | task10 | TBD | | @@ -100,22 +100,45 @@ Look at the following document. ## Slides - [[1] Introduction](http://nobuyuki-umetani.com/acg2024s/introduction.pdf) + - [[2] C++ programming](http://nobuyuki-umetani.com/acg2024s/cpp.pdf) + - [[3] Git+GitHub](http://nobuyuki-umetani.com/acg2024s/git.pdf) + - [[4] Digital Image](http://nobuyuki-umetani.com/acg2024s/digital_image.pdf) + - [[5] Rasterization in 2D](http://nobuyuki-umetani.com/acg2024s/rasterization_2d.pdf) + - [[6] Barycentric Coordinates](http://nobuyuki-umetani.com/acg2024s/barycentric_coordinates.pdf) + - [[7] Parametric Curve](http://nobuyuki-umetani.com/acg2024s/parametric_curve.pdf) + - [[8] Polynominal Root finding](http://nobuyuki-umetani.com/acg2024s/polynominal.pdf) + - [[9] Coordinate Transformation](http://nobuyuki-umetani.com/acg2024s/transformation.pdf) + - [[10] 2D Homogeneous Transformation](http://nobuyuki-umetani.com/acg2024s/transformation_homogeneous_2d.pdf) + - [[11] 3D Homogeneous Transformation](http://nobuyuki-umetani.com/acg2024s/transformation_homogeneous_3d.pdf) + - [[12] 3D Rasterization](http://nobuyuki-umetani.com/acg2024s/rasterization_3d.pdf) + - [[13] Graphics Pipeline](http://nobuyuki-umetani.com/acg2024s/graphics_pipeline.pdf) + - [[14] Shading](http://nobuyuki-umetani.com/acg2024s/shading.pdf) + - [[15] Subpixel Effect](http://nobuyuki-umetani.com/acg2024s/rasterization_subpixel.pdf) + - [[16] Implicit Modeling](http://nobuyuki-umetani.com/acg2024s/implicit_modeling.pdf) +- [[17]Ray Casting](http://nobuyuki-umetani.com/acg2024s/ray_casting.pdf) + +- [[18]Monte Carlo Integration](http://nobuyuki-umetani.com/acg2024s/monte_carlo_integration.pdf) + +- [[19]Ray Triangle Collision](http://nobuyuki-umetani.com/acg2024s/ray_triangle_collision.pdf) + + + ## Reading Material diff --git a/task06/CMakeLists.txt b/task06/CMakeLists.txt new file mode 100644 index 0000000..16afabc --- /dev/null +++ b/task06/CMakeLists.txt @@ -0,0 +1,35 @@ +# specify the version of cmake (intentionally using old version for those who cannot update CMake) +cmake_minimum_required(VERSION 3.10) + +############################# +# set C++ detail +enable_language(CXX) # we are using C++ +set(CMAKE_CXX_STANDARD 17) # we are using C++ 17 +set(CMAKE_CXX_STANDARD_REQUIRED TRUE) # we are using STL library + +############################# +# set project name + +project(task07) + +############################# +# define macro +add_definitions(-DPROJECT_SOURCE_DIR="${PROJECT_SOURCE_DIR}") + +############################# +# specifying libraries to use + +######################## +# include, build, and link + +include_directories( + ${PROJECT_SOURCE_DIR}/../external + ${PROJECT_SOURCE_DIR}/../external/eigen +) + +add_executable(${PROJECT_NAME} + main.cpp +) + +target_link_libraries(${PROJECT_NAME} +) \ No newline at end of file diff --git a/task06/README.md b/task06/README.md new file mode 100644 index 0000000..be389e8 --- /dev/null +++ b/task06/README.md @@ -0,0 +1,77 @@ +# Task06: Ambient Occlusion and BVH (Monte Carlo Integration, Importance Sampling) + +![preview](preview.png) + +**Deadline: May 30th (Thu) at 15:00pm** + +---- + +## Before Doing Assignment + +If you have not done the [task01](../task01), [task02](../task02) do it first to set up the C++ development environment. + +Follow [this document](../doc/submit.md) to submit the assignment, In a nutshell, before doing the assignment, +- make sure you synchronized the `main ` branch of your local repository to that of remote repository. +- make sure you created branch `task06` from `main` branch. +- make sure you are currently in the `task06` branch (use `git branch -a` command). + +Now you are ready to go! + +--- + +## Problem 1 (compilation practice) + +1. Build the code using cmake +2. Run the code **with release mode** + +The program will output `normal_map.png` that update the image below + +![problem1](normal_map.png) + +## Problem 2 (efficient search using BVH) + +The code in Problem 1 is very slow. It is because the computation of intersection betweeen a ray and a triangle mesh is **brute force**. + +Comment out the brute force intersection computation (from `line#146` to `line#158`) and implement code to evaluate BVH to accelerate the computation around `line #120`. + +The program output the computation time. Fill the table below to compare the timing before/after the acceleraion. Please make sure that you build the code **with release mode** + + + +| brute force | BVH | +| ----------- | ------ | +| ??? ms | ??? ms | + + + + +## Problem 3 (ambient occlusion) + +Now you have the code for fast ray-mesh intersection. Using that code, let's compute the ambient occlusion. + +- First, Comment out the `continue;` in the `main()` at `line #21` in the original code . This will enable the ambient occlusion computation. The problem will output `ao.png`. + +- Then fix the bug in `line #223`in the original code to correctly compute the ambient occlusion. The result should looks like `preview.png` at the begining of this document. + +- Finally, modify the name of the outut image as `ao_uniform.png` + + + +## Problem 4 (importance sampling) + +The computation of ambient occlusion is a bit noisy (i.e., the variance is high). Let's implement the importance sampling to reduce the variance. Write some code to compute the direction and PDF of **cosine-weighted sampling of hemisphere **around `line #53`. Modify the name of the output image `ao.png` as `ao_cosweight.png`. + +| Uniform sample | Cosine weighted sample | +| --------------------------- | ----------------------------- | +| ![problem2](ao_uniform.png) | ![problem3](ao_cosweight.png) | + + + +Observe that the variance (noise) is reduced using the importance sampling + + + + +## After Doing the Assignment + +After modify the code, push the code and submit a pull request. Make sure your pull request only contains the files you edited. Good luck! diff --git a/task06/ao_cosweight.png b/task06/ao_cosweight.png new file mode 100644 index 0000000..554e79d Binary files /dev/null and b/task06/ao_cosweight.png differ diff --git a/task06/ao_uniform.png b/task06/ao_uniform.png new file mode 100644 index 0000000..6680934 Binary files /dev/null and b/task06/ao_uniform.png differ diff --git a/task06/main.cpp b/task06/main.cpp new file mode 100644 index 0000000..3c50373 --- /dev/null +++ b/task06/main.cpp @@ -0,0 +1,255 @@ +#include +#include +#include +#include +#include +// +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include "stb_image_write.h" +#include "Eigen/Core" +#include "Eigen/Geometry" +// +#include "util.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +/** + * Rotation matrix such that z-axis will be direction of `nrm` + * @param nrm transfromed z-axis + * @return 3x3 Rotation matrix + */ +auto local_to_world_vector_transformation( + const Eigen::Vector3f &nrm) -> Eigen::Matrix3f { + auto basis_x = Eigen::Vector3f(1.f, 0.f, 0.f); + const auto basis_y = nrm.cross(basis_x).normalized(); + basis_x = basis_y.cross(nrm); + Eigen::Matrix3f loc2world; + loc2world << basis_x, basis_y, nrm; // initialization from 3 column vectors + return loc2world; +} + +/** + * sample a point on a unit hemisphere + * @param nrm up direction of the hemisphere + * @return sampled direction and its PDF + */ +auto sample_hemisphere( + const Eigen::Vector3f &nrm) -> std::pair { + const auto unirand = Eigen::Vector2f::Random() * 0.5f + Eigen::Vector2f(0.5, 0.5); + const float phi = 2.f * float(M_PI) * unirand.y(); + + // the code to uniformly sample hemisphere (z-up) + const float z = unirand.x(); + const float r = std::sqrt(1.f - z * z); + auto dir_loc = Eigen::Vector3f( // direction in normal coordinate + r * std::cos(phi), + r * std::sin(phi), + z); + float pdf = 0.5f / float(M_PI); + + // For Problem 4, write some code below to sample hemisphere with cosign weight + // (i.e., the sampling frequency is higher at the top) + + + // end of Problem 4. Do not modify the two lines below + const auto dir_out = local_to_world_vector_transformation(nrm) * dir_loc; // rotate the sample (zup -> nrm) + return {dir_out, pdf}; +} + +/** + * ray triangle intersection + * @param ray_org origin of the ray + * @param ray_dir direction of the ray + * @param i_tri index of the triangle + * @param tri2vtx list of triangle index + * @param vtx2xyz list of vertex coordinates + * @return std::nullopt if there is no intersection, otherwise returns a pair of position and normal + */ +auto ray_triangle_intersection( + const Eigen::Vector3f &ray_org, + const Eigen::Vector3f &ray_dir, + unsigned int i_tri, + const Eigen::MatrixX3i &tri2vtx, + const Eigen::MatrixX3f &vtx2xyz) +-> std::optional> { + auto tet_volume = []( + const Eigen::Vector3f &p0, + const Eigen::Vector3f &p1, + const Eigen::Vector3f &p2, + const Eigen::Vector3f &p3) { return (p1 - p0).cross(p2 - p0).dot(p3 - p0) / 6.f; }; + auto p0 = vtx2xyz.row(tri2vtx(i_tri, 0)); + auto p1 = vtx2xyz.row(tri2vtx(i_tri, 1)); + auto p2 = vtx2xyz.row(tri2vtx(i_tri, 2)); + float v0 = tet_volume(p1, p2, ray_org + ray_dir, ray_org); + float v1 = tet_volume(p2, p0, ray_org + ray_dir, ray_org); + float v2 = tet_volume(p0, p1, ray_org + ray_dir, ray_org); + if (v0 < 0.f || v1 < 0.f || v2 < 0.f) { return std::nullopt; } + float sum_inv = 1.f / (v0 + v1 + v2); + const Eigen::Vector3f q0 = (v0 * p0 + v1 * p1 + v2 * p2) * sum_inv; + const Eigen::Vector3f n0 = (p1 - p0).cross(p2 - p0).normalized(); + return std::make_pair(q0, n0); +} + +/** + * search ray-triangle intersection in the triangles under the specified branch of BVH + * @param[in,out] is_hit flag true if there is collision + * @param[in,out] hit_depth update the minimum depth of the intersection location + * @param[in, out] hit_pos update the position of the intersection location + * @param[in, out] hit_normal update normal at the intersection location + * @param[in] i_bvhnode index of BVH node of branch to search intersection + * @param[in] ray_org ray origin + * @param[in] ray_dir ray direction + * @param[in] tri2vtx triangle index + * @param[in] vtx2xyz list of coordinates + * @param[in] bvhnodes list of BVH nodes + */ +void search_collision_in_bvh( + bool &is_hit, + float &hit_depth, + Eigen::Vector3f &hit_pos, + Eigen::Vector3f &hit_normal, + unsigned int i_bvhnode, + const Eigen::Vector3f &ray_org, + const Eigen::Vector3f &ray_dir, + const Eigen::MatrixX3i &tri2vtx, + const Eigen::MatrixX3f &vtx2xyz, + const std::vector &bvhnodes) { + // For problem 2, implement some code here to evaluate BVH + // hint: use following function + // bvhnodes[i_bvhnode].intersect_bv(ray_org, ray_dir) + + if (bvhnodes[i_bvhnode].is_leaf()) { // this is leaf node + const unsigned int i_tri = bvhnodes[i_bvhnode].i_node_left; + // do something + } else { // this is branch node + unsigned int i_node_right = bvhnodes[i_bvhnode].i_node_right; + unsigned int i_node_left =bvhnodes[i_bvhnode].i_node_left; + // do something (hint recursion) + } +} + +auto find_intersection_between_ray_and_triangle_mesh( + const Eigen::Vector3f &ray_org, + const Eigen::Vector3f &ray_dir, + const Eigen::MatrixX3i &tri2vtx, + const Eigen::MatrixX3f &vtx2xyz, + const std::vector &bvhnodes) +-> std::optional> { + bool is_hit = false; + float hit_depth = 1000.; + Eigen::Vector3f hit_pos; + Eigen::Vector3f hit_normal; + + // for Problem 2,3,4, comment out from here + for (unsigned int i_tri = 0; i_tri < tri2vtx.rows(); ++i_tri) { + const auto res = ray_triangle_intersection(ray_org, ray_dir, i_tri, tri2vtx, vtx2xyz); + if (!res) { continue; } + const auto& [q0,n0] = res.value(); + const float depth = (q0 - ray_org).dot(ray_dir); + if (hit_depth > depth) { + is_hit = true; + hit_depth = depth; + hit_pos = q0; + hit_normal = n0; + } + } + // comment out end + + // do not edit from here + search_collision_in_bvh( + is_hit, hit_depth, hit_pos, hit_normal, + 0, // root node index + ray_org, ray_dir, tri2vtx, vtx2xyz, bvhnodes); + // + if (!is_hit) { return std::nullopt; } + return std::make_pair(hit_pos, hit_normal); +} + +auto get_ray_from_camera( + unsigned int width, unsigned int height, + unsigned int iw, unsigned int ih) -> std::pair { + auto cam_ray_src = Eigen::Vector3f(0.5, 0.5, 2.0); // focus point + float ndc_x = ((float(iw) + 0.5f) * 2.f / float(width) - 1.f); // normalized device x-coordinate [-1, +1] + float ndc_y = (1.f - (float(ih) + 0.5f) * 2.f / float(height)); // normalized device y-coordinate [-1, +1] + float sensor_size = 0.25; + Eigen::Vector3f position_on_sensor = Eigen::Vector3f(ndc_x * sensor_size, ndc_y * sensor_size, -1.0) + cam_ray_src; + Eigen::Vector3f cam_ray_dir = (position_on_sensor - cam_ray_src).normalized(); + return {cam_ray_src, cam_ray_dir}; +} + +int main() { + Eigen::MatrixX3f vtx2xyz; + Eigen::MatrixX3i tri2vtx; + std::vector bvhnodes; + acg::load_scene(vtx2xyz, tri2vtx, bvhnodes); + + const unsigned int img_width = 100; + const unsigned int img_height = 100; + // + std::vector img_data_nrm(img_height * img_width * 3, 0.f); + std::vector img_data_ao(img_height * img_width, 0.f); + // + std::chrono::system_clock::time_point start = std::chrono::system_clock::now(); // record starting time + for (unsigned int iw = 0; iw < img_width; ++iw) { + for (unsigned int ih = 0; ih < img_height; ++ih) { + const auto[cam_ray_src, cam_ray_dir] = get_ray_from_camera(img_width, img_height, iw, ih); + const auto& res = find_intersection_between_ray_and_triangle_mesh( + cam_ray_src, cam_ray_dir, tri2vtx, vtx2xyz, bvhnodes); + // draw normal map + if (!res) { + img_data_nrm[(ih * img_width + iw) * 3 + 0] = 0.f; + img_data_nrm[(ih * img_width + iw) * 3 + 1] = 0.f; + img_data_nrm[(ih * img_width + iw) * 3 + 2] = 0.f; + } else { + const auto&[pos, nrm] = res.value(); // position and normal of the first hit point + img_data_nrm[(ih * img_width + iw) * 3 + 0] = nrm.x() * 0.5f + 0.5f; + img_data_nrm[(ih * img_width + iw) * 3 + 1] = nrm.y() * 0.5f + 0.5f; + img_data_nrm[(ih * img_width + iw) * 3 + 2] = nrm.z() * 0.5f + 0.5f; + } + continue; // comment out here for Problem 3,4 + // + if (res) { // ambient occlusion computation + const unsigned int num_sample_ao = 100; + float sum = 0; + for (unsigned int i_sample = 0; i_sample < num_sample_ao; ++i_sample) { + const auto&[pos, nrm] = res.value(); // position and normal of the first hit point + Eigen::Vector3f pos0 = pos + nrm * 0.001f; // offset the position in the direction of normal + const auto[dir, pdf] = sample_hemisphere(nrm); // direction of the sampled light position and its PDF + const auto res1 = find_intersection_between_ray_and_triangle_mesh( + pos0, dir, tri2vtx, vtx2xyz, bvhnodes); + if (!res1) { // if the ray doe not hit anything + sum += 1.f; // Problem 3: This is a bug. write some correct code (hint: use `dir.dot(nrm)`, `pdf`, `M_PI`). + } + } + img_data_ao[ih * img_width + iw] = sum / float(num_sample_ao); // do not change + } + } + } + std::chrono::system_clock::time_point end = std::chrono::system_clock::now(); // record end time + auto elapsed = std::chrono::duration_cast(end - start).count(); + std::cout << "total computation time: " << elapsed << "ms" << std::endl; + + { + std::vector img_data_uchar(img_width * img_height * 3, 0); + for (unsigned int i = 0; i < img_width * img_height; ++i) { + img_data_uchar[i * 3 + 0] = img_data_nrm[i * 3 + 0] * 255.f; + img_data_uchar[i * 3 + 1] = img_data_nrm[i * 3 + 1] * 255.f; + img_data_uchar[i * 3 + 2] = img_data_nrm[i * 3 + 2] * 255.f; + } + stbi_write_png( + (std::filesystem::path(PROJECT_SOURCE_DIR) / "normal_map.png").string().c_str(), + img_width, img_height, 3, img_data_uchar.data(), 0); + } + { + std::vector img_data_uchar(img_width * img_height, 0); + for (unsigned int i = 0; i < img_width * img_height; ++i) { + img_data_uchar[i] = img_data_ao[i] * 255.f; + } + stbi_write_png( + (std::filesystem::path(PROJECT_SOURCE_DIR) / "ao.png").string().c_str(), + img_width, img_height, 1, img_data_uchar.data(), 0); + } +} + diff --git a/task06/normal_map.png b/task06/normal_map.png new file mode 100644 index 0000000..2bbe647 Binary files /dev/null and b/task06/normal_map.png differ diff --git a/task06/preview.png b/task06/preview.png index c3a0adb..ff1ea26 100644 Binary files a/task06/preview.png and b/task06/preview.png differ diff --git a/task06/util.h b/task06/util.h new file mode 100644 index 0000000..dd00452 --- /dev/null +++ b/task06/util.h @@ -0,0 +1,188 @@ +#ifndef UTIL_H_ +#define UTIL_H_ + +#include + +namespace acg { + +class BvhNode { + public: + unsigned int i_node_left; + unsigned int i_node_right; + Eigen::Vector3f v_min; + Eigen::Vector3f v_max; + public: + /** + * check if this BVH node is leaf or not. If it is leaf, `i_node_left` will be the corresponding triangle index + * @return true if its leaf + */ + [[nodiscard]] bool is_leaf() const { + return i_node_right == UINT_MAX; + } + /** + * intersection of ray against the bounding volume + * @param ray_org ray origin + * @param ray_dir ray direction + * @return true if there is intersection + */ + [[nodiscard]] bool intersect_bv( + const Eigen::Vector3f& ray_org, + const Eigen::Vector3f& ray_dir) const { + float tmin = -INFINITY; + float tmax = INFINITY; + for(int i_dim=0;i_dim<3;++i_dim) { + if (std::fabs(ray_dir[i_dim]) > 1.0e-10f) { + const float t1 = (v_min[i_dim] - ray_org[i_dim]) / ray_dir[i_dim]; + const float t2 = (v_max[i_dim] - ray_org[i_dim]) / ray_dir[i_dim]; + tmin = std::max(tmin, std::min(t1, t2)); + tmax = std::min(tmax, std::max(t1, t2)); + } + else if( ray_org[i_dim] < v_min[i_dim] || ray_org[i_dim] > v_max[i_dim] ) { + return false; + } + } + return tmax >= tmin && tmax >= 0.0; + } +}; + +// don't need to understand +void build_bvh_topology( + unsigned int i_start, + unsigned int i_end, + bool is_right, + std::vector &bvhnodes, + unsigned int num_branch) { + if (i_end == i_start + 1) { + if (is_right) { + bvhnodes[i_start].i_node_left = num_branch + i_start; + bvhnodes[i_start].i_node_right = num_branch + i_end; + } else { + bvhnodes[i_end].i_node_left = num_branch + i_start; + bvhnodes[i_end].i_node_right = num_branch + i_end; + } + return; + } + unsigned int i_split = (i_start + i_end - 1) / 2; + if (is_right) { + bvhnodes[i_start].i_node_left = i_split; + bvhnodes[i_start].i_node_right = i_split + 1; + } else { + bvhnodes[i_end].i_node_left = i_split; + bvhnodes[i_end].i_node_right = i_split + 1; + } + build_bvh_topology(i_start, i_split, false, bvhnodes, num_branch); + build_bvh_topology(i_split + 1, i_end, true, bvhnodes, num_branch); +} + +// set the geometry to the bounding volume +void set_bvh_geometry( + unsigned int i_node, + std::vector &bvhnodes, + Eigen::MatrixX3i &tri2vtx, + Eigen::MatrixX3f &vtx2xyz) +{ + if (bvhnodes[i_node].is_leaf()) { + unsigned int i_tri = bvhnodes[i_node].i_node_left; + auto p0 = vtx2xyz.row(tri2vtx(i_tri, 0)); + auto p1 = vtx2xyz.row(tri2vtx(i_tri, 1)); + auto p2 = vtx2xyz.row(tri2vtx(i_tri, 2)); + bvhnodes[i_node].v_min[0] = std::min({p0[0], p1[0], p2[0]}); + bvhnodes[i_node].v_min[1] = std::min({p0[1], p1[1], p2[1]}); + bvhnodes[i_node].v_min[2] = std::min({p0[2], p1[2], p2[2]}); + bvhnodes[i_node].v_max[0] = std::max({p0[0], p1[0], p2[0]}); + bvhnodes[i_node].v_max[1] = std::max({p0[1], p1[1], p2[1]}); + bvhnodes[i_node].v_max[2] = std::max({p0[2], p1[2], p2[2]}); + } else { + set_bvh_geometry(bvhnodes[i_node].i_node_left, bvhnodes, tri2vtx, vtx2xyz); + set_bvh_geometry(bvhnodes[i_node].i_node_right, bvhnodes, tri2vtx, vtx2xyz); + auto min_left = bvhnodes[bvhnodes[i_node].i_node_left].v_min; + auto max_left = bvhnodes[bvhnodes[i_node].i_node_left].v_max; + auto min_right = bvhnodes[bvhnodes[i_node].i_node_right].v_min; + auto max_right = bvhnodes[bvhnodes[i_node].i_node_right].v_max; + bvhnodes[i_node].v_min[0] = std::min(min_left[0], min_right[0]); + bvhnodes[i_node].v_min[1] = std::min(min_left[1], min_right[1]); + bvhnodes[i_node].v_min[2] = std::min(min_left[2], min_right[2]); + bvhnodes[i_node].v_max[0] = std::max(max_left[0], max_right[0]); + bvhnodes[i_node].v_max[1] = std::max(max_left[1], max_right[1]); + bvhnodes[i_node].v_max[2] = std::max(max_left[2], max_right[2]); + } +} + +// don't need to understand this... +uint16_t int_coord_from_morton(uint16_t i_quad) { + return (i_quad & 0x0001) + + ((i_quad & 0x0004) >> 1) + + ((i_quad & 0x0010) >> 2) + + ((i_quad & 0x0040) >> 3) + + ((i_quad & 0x0100) >> 4) + + ((i_quad & 0x0400) >> 5) + + ((i_quad & 0x1000) >> 6) + + ((i_quad & 0x4000) >> 7); +} + +void load_scene( + Eigen::MatrixX3f &vtx2xyz, + Eigen::MatrixX3i &tri2vtx, + std::vector &bvhnodes) +{ + unsigned int num_lev = 7; + auto num_div = static_cast(pow(2, num_lev)); + std::cout << "number of elements on an edge of square mesh: " << num_div << std::endl; + const unsigned int size = num_div + 1; + const unsigned int num_vtx = size * size; + vtx2xyz.resize(num_vtx, 3); + float theta = -3.14 / 4.0; + for (unsigned int i_h = 0; i_h < size; ++i_h) { + float y = float(i_h) / float(num_div); + for (unsigned int i_w = 0; i_w < size; ++i_w) { + float x = float(i_w) / float(num_div); + float r = std::sqrt((x-0.5f)*(x-0.5f)+(y-0.5f)*(y-0.5f)); + float z = 0.5f*exp(-10.f*r*r)*sin(30.f*r)-0.5f; + vtx2xyz(i_h * size + i_w, 0) = x; + vtx2xyz(i_h * size + i_w, 1) = std::cos(theta)*y - std::sin(theta)*z + 0.4f; + vtx2xyz(i_h * size + i_w, 2) = std::sin(theta)*y + std::cos(theta)*z; + } + } + const unsigned int num_quad = num_div * num_div; + const unsigned int num_tri = num_quad * 2; + tri2vtx.resize(num_tri, 3); + for (unsigned int i_quad = 0; i_quad < num_quad; ++i_quad) { + unsigned int i_w = int_coord_from_morton(i_quad); + unsigned int i_h = int_coord_from_morton(i_quad>>1); + unsigned int i0_tri = (i_h * num_div + i_w) * 2 + 0; + tri2vtx(i0_tri, 0) = i_h * size + i_w; + tri2vtx(i0_tri, 1) = i_h * size + (i_w + 1); + tri2vtx(i0_tri, 2) = (i_h + 1) * size + (i_w + 1); + unsigned int i1_tri = (i_h * num_div + i_w) * 2 + 1; + tri2vtx(i1_tri, 0) = i_h * size + i_w; + tri2vtx(i1_tri, 1) = (i_h + 1) * size + (i_w + 1); + tri2vtx(i1_tri, 2) = (i_h + 1) * size + i_w; + } + bvhnodes.resize(num_quad * 2 - 1 + num_quad * 2); + build_bvh_topology(0, num_quad - 1, true, bvhnodes, num_quad - 1); + for (unsigned int i_quad = 0; i_quad < num_quad; ++i_quad) { + unsigned int i_node = i_quad + num_quad - 1; + bvhnodes[i_node].i_node_left = num_quad * 2 - 1 + i_quad * 2 + 0; + bvhnodes[i_node].i_node_right = num_quad * 2 - 1 + i_quad * 2 + 1; + } + for (unsigned int i_tri = 0; i_tri < num_tri; ++i_tri) { + unsigned int i_node = num_quad * 2 - 1 + i_tri; + bvhnodes[i_node].i_node_left = i_tri; + bvhnodes[i_node].i_node_right = -1; + } + set_bvh_geometry(0, bvhnodes, tri2vtx, vtx2xyz); + + /* + for(int i_node=0;i_node