diff --git a/.gitignore b/.gitignore index 0c2705c8..d5dc2bcb 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,5 @@ target NJUCG.mdj cmake-build* .idea -./doxygen-doc \ No newline at end of file +./doxygen-doc +scenes/**/*.hdr \ No newline at end of file diff --git a/.gitmodules b/.gitmodules index a1e9f05f..cce648f8 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,3 +14,6 @@ [submodule "tools/ScenePreviewer"] path = tools/ScenePreviewer url = https://github.com/r1ckhu/Moer_ScenePreviewer.git +[submodule "ext/autodiff"] + path = ext/autodiff + url = https://github.com/autodiff/autodiff.git \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 75dd3014..3940cfbe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,8 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON) if (MSVC) add_compile_options("/MP") add_compile_options("$<$:/utf-8>") - add_compile_options("$<$:/utf-8>") + add_compile_options("$<$:/utf-8>") + add_compile_options("/bigobj") # set unicode as character set # add_definitions(-DUNICODE -D_UNICODE) endif (MSVC) @@ -18,6 +19,7 @@ endif (FM_SPEED_DEFAULT) option(EMBREE_USE_TBB "Enable TBB in embree." OFF) option(USE_SAMPLED_SPECTRUM "" OFF) +option(ENABLE_GPISMEDIUM "enable gpis medium feature" ON) set(MESH_LOADER "Tinyobjloader" CACHE STRING "Select mesh loader") set_property(CACHE MESH_LOADER PROPERTY STRINGS Assimp Tinyobjloader) diff --git a/README.md b/README.md index eb926469..71114765 100644 --- a/README.md +++ b/README.md @@ -22,6 +22,8 @@ Moer is developed by NJUMeta (www.njumeta.com) from Nanjing University. ![hair](https://yjpp.oss-cn-hangzhou.aliyuncs.com/uPic/C6966B38E3AD16F8AC9CC9D86C73921E.jpg) ![teapot](https://yjpp.oss-cn-hangzhou.aliyuncs.com/uPic/tea-pot.png) ![disney](https://yjpp.oss-cn-hangzhou.aliyuncs.com/uPic/disneybsdf.png) +![micrograin-teapot](https://github.com/Cchen-77/images/blob/main/micrograin.jpg) +![gpis](https://github.com/Cchen-77/images/blob/main/gpis.jpg) ## Checklist ### Accelerator @@ -63,6 +65,7 @@ Moer is developed by NJUMeta (www.njumeta.com) from Nanjing University. ### Media - [X] homogeneous - [X] hetergeneous +- [X] gaussian process implicit surface - [ ] others ### BxDF @@ -73,6 +76,7 @@ Moer is developed by NJUMeta (www.njumeta.com) from Nanjing University. - [X] disney - [X] plastic - [X] hair +- [X] micrograin for porous layer - [ ] others ### Shape diff --git a/ext/CMakeLists.txt b/ext/CMakeLists.txt index 51235051..2425b97e 100644 --- a/ext/CMakeLists.txt +++ b/ext/CMakeLists.txt @@ -180,4 +180,16 @@ if (NOT TARGET stb) target_link_libraries(${TARGET_NAME} PUBLIC stb) target_include_directories(stb PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/stb) set_ide_folder(stb ${3DPARTY_FOLDER}) +endif() + +if(NOT TARGET autodiff) + message(STATUS "============start config autodiff============") + set(AUTODIFF_BUILD_TESTS OFF) + set(AUTODIFF_BUILD_PYTHON OFF) + set(AUTODIFF_BUILD_EXAMPLES OFF) + set(AUTODIFF_BUILD_DOCS OFF) + add_subdirectory(autodiff) + target_link_libraries(${TARGET_NAME} PUBLIC autodiff) + target_include_directories(${TARGET_NAME} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/autodiff) + set_ide_folder(autodiff ${3DPARTY_FOLDER}) endif() \ No newline at end of file diff --git a/ext/autodiff b/ext/autodiff new file mode 160000 index 00000000..cb2d6e8b --- /dev/null +++ b/ext/autodiff @@ -0,0 +1 @@ +Subproject commit cb2d6e8b9fa98b7b8cf3ffaefa231943c56361ae diff --git a/scenes/gpis-1/scene.json b/scenes/gpis-1/scene.json new file mode 100644 index 00000000..0b3e0eab --- /dev/null +++ b/scenes/gpis-1/scene.json @@ -0,0 +1,142 @@ +{ + "renderer" : { + "spp" : 1, + "output_file" : "gpis" + }, + + "mediums" : [ + { + "name" : "gpis", + "type" : "gpis", + "num_sample_points": 8, + "marching_step_size": 10000, + "gaussian_process":{ + "type":"std", + "mean":{ + "type":"procedural", + "func":"knob", + "transform":{ + "position":[0,0,0], + "scale":[1,1,1], + "rotation":[0,30,0] + } + }, + "covariance":{ + "type":"squared_exponential", + "lengthScale":0.1, + "sigma":0.01 + } + }, + "phase":{ + "type":"lambert", + "albedo":[0.9, 0.3, 0.35] + } + } + ], + + "materials" : [ + { + "name" : "leftWall", + "type" : "lambert", + "albedo" : [0.63, 0.0065, 0.05] + }, + { + "name" : "rightWall", + "type" : "lambert", + "albedo" : [0.14, 0.45, 0.091] + }, + { + "name" : "floor", + "type" : "lambert", + "albedo" : [0.725, 0.71, 0.68] + }, + { + "name" : "ceiling", + "type" : "lambert", + "albedo" : [0.725, 0.71, 0.68] + }, + { + "name" : "backWall", + "type" : "lambert", + "albedo" : [0.725, 0.71, 0.68] + }, + { + "name" : "light", + "type" : "lambert", + "albedo" : [0, 0, 0] + }, + { + "name" : "interface", + "type" : "null", + "in_medium" : "gpis" + } + ], + + "entities" : [ + { + "type" : "quad", + "base" : [0, 0, 0], + "edge0" : [2, 0, 0], + "edge1" : [0, 0, 2], + "material" : "floor" + }, + { + "type" : "quad", + "base" : [0, 2, 0], + "edge0" : [0, 0, -2], + "edge1" : [-2, 0, 0], + "material" : "ceiling" + }, + { + "type" : "quad", + "base" : [0, 1, -1], + "edge0" : [-2, 0, 0], + "edge1" : [0, 2, 0], + "material" : "backWall" + }, + { + "type" : "quad", + "base" : [-1, 1, 0], + "edge0" : [0, 0, 2], + "edge1" : [0, 2, 0], + "material" : "leftWall" + }, + { + "type" : "quad", + "base" : [1, 1, 0], + "edge0" : [0, 0, -2], + "edge1" : [0, 2, 0], + "material" : "rightWall" + }, + { + "type" : "quad", + "transform": { + "position": [-0.005, 1.98, -0.03], + "scale": [0.47, 0.1786, 0.38], + "rotation": [0, 180, 180] + }, + "emission" : [30, 28, 20], + "material" : "light" + }, + { + "type" : "cube", + "material" : "interface", + "transform" : { + "scale" : [3, 3, 3], + "position" : [0, 1, 0] + } + } + ], + + "camera" : { + "resolution" : [512, 512], + "type" : "pinhole", + "fov" : 19.5, + + "transform" : { + "position" : [0, 1, 7], + "up" : [0, 1, 0], + "look_at" : [0, 1, 0] + } + } +} \ No newline at end of file diff --git a/scenes/gpis-2/scene.json b/scenes/gpis-2/scene.json new file mode 100644 index 00000000..e217bc4c --- /dev/null +++ b/scenes/gpis-2/scene.json @@ -0,0 +1,163 @@ +{ + "renderer" : { + "spp" : 32, + "output_file" : "gpis" + }, + + "mediums" : [ + { + "name" : "gpis-hard", + "type" : "gpis", + "marching_num_sample_points":64, + "marching_step_size": 10000, + "gaussian_process":{ + "type":"std", + "mean":{ + "type":"procedural", + "func":"knob", + "transform":{ + "position":[-0.8,0.37,0], + "scale":[0.5,0.5,0.5], + "rotation":[0,0,0] + } + }, + "covariance":{ + "type":"squared_exponential", + "lengthScale":0.1, + "sigma":0.01 + }, + "global_condition":{ + "enable":false, + "points":[ + [0,0,0], + [1,1,1] + ], + "derivative_types":[ + "none", + "first" + ], + "derivative_dirs":[ + [0,0,0], + [1,0,0] + ], + "values":[ + 0, + 1 + ] + } + }, + "phase":{ + "type":"lambert", + "albedo":[0.9, 0.3, 0.35] + } + }, + { + "name" : "gpis-soft", + "type" : "gpis", + "marching_num_sample_points":32, + "marching_step_size": 10000, + "gaussian_process":{ + "type":"std", + "mean":{ + "type":"procedural", + "func":"knob", + "transform":{ + "position":[0.8,0.37,0], + "scale":[0.5,0.5,0.5], + "rotation":[0,0,0] + } + }, + "covariance":{ + "type":"squared_exponential", + "lengthScale":0.1, + "sigma":0.05 + } + }, + "phase":{ + "type":"lambert", + "albedo":[0.9, 0.3, 0.35] + } + }, + { + "name" : "homogeneous", + "type" : "homogeneous", + "sigmaT" : [20, 20, 20], + "albedo" : [0.8, 0.8, 0.8] + } + ], + + "materials" : [ + { + "name": "floor", + "albedo": { + "type": "checker", + "on_color": [ + 0.725, + 0.71, + 0.68 + ], + "off_color": [ + 0.325, + 0.31, + 0.25 + ], + "res_u": 20, + "res_v": 20 + }, + "type": "lambert" + }, + { + "name" : "interface-hard", + "type" : "null", + "in_medium":"gpis-hard" + }, + { + "name" : "interface-soft", + "type" : "null", + "in_medium":"gpis-soft" + } + ], + + "entities" : [ + { + "type" : "quad", + "base" : [0, 0, 0], + "edge0" : [6, 0, 0], + "edge1" : [0, 0, 6], + + "material" : "floor" + }, + { + "type" : "infinite_sphere", + "emission" : "textures/envmap.hdr" + }, + { + "type" : "cube", + "material" : "interface-hard", + "transform" : { + "scale" : [2, 2, 2], + "position" : [-0.8, 1.05, 0] + } + }, + { + "type" : "cube", + "material" : "interface-soft", + "transform" : { + "scale" : [2, 2, 2], + "position" : [0.8, 1.05, 0] + } + } + ], + + "camera" : { + "resolution" : [1024, 1024], + "type" : "pinhole", + "fov" : 19.5, + + "transform" : { + "position" : [0, 4, 7], + "up" : [0, 1, 0], + "look_at" : [0, 0, 0] + } + } +} \ No newline at end of file diff --git a/scenes/gpis-2/textures/envmap.hdr b/scenes/gpis-2/textures/envmap.hdr new file mode 100644 index 00000000..e93e1cf5 Binary files /dev/null and b/scenes/gpis-2/textures/envmap.hdr differ diff --git a/scenes/heterogeneous2/models/smoke.nvdb b/scenes/heterogeneous2/models/smoke.nvdb new file mode 100644 index 00000000..342efd59 Binary files /dev/null and b/scenes/heterogeneous2/models/smoke.nvdb differ diff --git a/scenes/heterogeneous2/scene.json b/scenes/heterogeneous2/scene.json new file mode 100644 index 00000000..f8bff837 --- /dev/null +++ b/scenes/heterogeneous2/scene.json @@ -0,0 +1,99 @@ +{ + "renderer" : { + "spp" : 1, + "output_file" : "smoke" + }, + + "mediums" : [ + { + "name" : "gpis", + "type" : "gpis", + "marching_num_sample_points": 8, + "marching_step_size": 0.1, + "gaussian_process":{ + "type":"std", + "mean":{ + "type":"procedural", + "func":"knob", + "transform":{ + "position":[0,1,0], + "scale":[1,1,1], + "rotation":[0,0,0] + } + }, + "covariance":{ + "type":"squared_exponential", + "lengthScale":0.1, + "sigma":0.01 + } + }, + "phase":{ + "type":"lambert", + "albedo":[0.9, 0.3, 0.35] + } + }, + { + "name" : "homogeneous", + "type" : "homogeneous", + "sigmaT" : [20, 20, 20], + "albedo" : [0.8, 0.8, 0.8] + }, + { + "name" : "smoke", + "type" : "heterogeneous", + "filepath" : "models/smoke.nvdb", + "sigma_scale" : 10000, + "transform" : { + "position" : [0, 0.05, 0], + "scale" : 0.03 + } + } + ], + + "materials" : [ + { + "name" : "floor", + "type" : "lambert", + "albedo" : [0.725, 0.71, 0.68] + }, + { + "name" : "interface", + "type" : "null", + "in_medium":"smoke" + } + ], + + "entities" : [ + { + "type" : "quad", + "base" : [0, 0, 0], + "edge0" : [3, 0, 0], + "edge1" : [0, 0, 3], + "material" : "floor" + }, + { + "type" : "infinite_sphere", + "emission" : "textures/envmap.hdr" + }, + { + "type" : "cube", + "material" : "interface", + "transform" : { + "scale" : [2, 2, 2], + "position" : [0, 1.05, 0] + } + } + ], + + "camera" : { + "resolution" : [1024, 1024], + "type" : "pinhole", + "fov" : 19.5, + + "transform" : { + "position" : [2, 7, 7], + "up" : [0, 1, 0], + "look_at" : [0, 0, 0] + } + } +} \ No newline at end of file diff --git a/scenes/micrograin-teapot/scene.json b/scenes/micrograin-teapot/scene.json index 5bc2ac53..795753f8 100644 --- a/scenes/micrograin-teapot/scene.json +++ b/scenes/micrograin-teapot/scene.json @@ -11,8 +11,8 @@ "type":"plastic" }, "micrograinType":"conductor", - "tau0":0.27, - "beta":0.1, + "tau0":0.5, + "beta":0.8, "R0":[0.588,0.294,0.001], "k":[0,0,0] }, @@ -137,7 +137,7 @@ "enable_resume_render": false, "stratified_sampler": false, "scene_bvh": true, - "spp": 1024, + "spp": 512, "spp_step": 16, "checkpoint_interval": "0", "timeout": "0", diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b60d25f0..b17d5f15 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -38,6 +38,10 @@ if (USE_SAMPLED_SPECTRUM) target_compile_definitions(${TARGET_NAME} PRIVATE USING_SAMPLED_SPECTRUM) endif() +if(ENABLE_GPISMEDIUM) + target_compile_definitions(${TARGET_NAME} PRIVATE ENABLE_GPISMEDIUM) +endif() + target_compile_definitions(${TARGET_NAME} PRIVATE CMAKE_DEF_SPECTRUM_SAMPLES=60) set_target_properties(${TARGET_NAME} PROPERTIES DEBUG_POSTFIX "_d") diff --git a/src/CoreLayer/Geometry/Vector.h b/src/CoreLayer/Geometry/Vector.h index af2b88d7..eb3a68fc 100644 --- a/src/CoreLayer/Geometry/Vector.h +++ b/src/CoreLayer/Geometry/Vector.h @@ -245,6 +245,18 @@ struct TVector3 { bool isZero() const { return x==0 && y==0 && z==0; } + + TVector2 xy() const { + return TVector2(x, y); + } + + TVector2 yz() const { + return TVector2(y, z); + } + + TVector2 xz() const { + return TVector2(x, z); + } }; diff --git a/src/FunctionLayer/GaussianProcess/GPFunctions.cpp b/src/FunctionLayer/GaussianProcess/GPFunctions.cpp new file mode 100644 index 00000000..2d99cf55 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GPFunctions.cpp @@ -0,0 +1,81 @@ +#include "GPFunctions.h" + +double CovarianceFunction::dcov_dx(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirX) const { + autodiff::Vector3real2nd px{pointX.x, pointX.y, pointX.z}; + autodiff::Vector3real2nd py{pointY.x, pointY.y, pointY.z}; + Eigen::Array3d vx{ddirX.x, ddirX.y, ddirX.z}; + Eigen::Array3d vy(0.); + auto dfdv = autodiff::derivatives([this](const autodiff::Vector3real2nd &x, const autodiff::Vector3real2nd &y) { return cov(x, y); }, autodiff::along(vx, vy), autodiff::at(px, py)); + return dfdv[1]; +} +double CovarianceFunction::dcov_dy(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirY) const { + autodiff::Vector3real2nd px{pointX.x, pointX.y, pointX.z}; + autodiff::Vector3real2nd py{pointY.x, pointY.y, pointY.z}; + Eigen::Array3d vx(0.); + Eigen::Array3d vy{ddirY.x, ddirY.y, ddirY.z}; + auto dfdv = autodiff::derivatives([this](const autodiff::Vector3real2nd &x, const autodiff::Vector3real2nd &y) { return cov(x, y); }, autodiff::along(vx, vy), autodiff::at(px, py)); + return dfdv[1]; +} +double CovarianceFunction::dcov2_dxdy(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirX, const Vec3d &ddirY) const { + autodiff::Vector3dual2nd px{pointX.x, pointX.y, pointX.z}; + autodiff::Vector3dual2nd py{pointY.x, pointY.y, pointY.z}; + Eigen::Matrix3d hess = autodiff::hessian([&](const autodiff::Vector3dual2nd &px, const autodiff::Vector3dual2nd &py) { return cov(px, py); }, wrt(px, py), at(px, py)).block(3, 0, 3, 3); + Eigen::Array3d vx{ddirX.x, ddirX.y, ddirX.z}; + Eigen::Array3d vy{ddirY.x, ddirY.y, ddirY.z}; + double res = vy.transpose().matrix() * hess * vx.matrix(); + return res; +} + +ProceduralMean::ProceduralMean(const Json &json) { + func = SdfFunctions::funcStringToEnum(json["func"]); + + Json transform = getOptional(json, "transform", Json()); + + Point3d transformPos = getOptional(transform, "position", Point3d(0.0)); + transformMatrix.setTranslate(transformPos.x, transformPos.y, transformPos.z); + Vec3d transformScale = getOptional(transform, "scale", Vec3d(1, 1, 1)); + transformMatrix.setScale(transformScale.x, transformScale.y, transformScale.z); + Vec3d transformRot = getOptional(transform, "rotation", Vec3d(0, 0, 0)); + transformMatrix.setRotateEuler(Angle(transformRot.x, Angle::EAngleType::ANGLE_DEG), + Angle(transformRot.y, Angle::EAngleType::ANGLE_DEG), + Angle(transformRot.z, Angle::EAngleType::ANGLE_DEG), + EulerType::EULER_YZX); + + invTransformMatrix = transformMatrix.getInverse(); + + scale = getOptional(json, "scale", 1.); + offset = getOptional(json, "offset", 0.); +} + +double ProceduralMean::mean(const Point3d &point) const { + auto invTransformedPoint = invTransformMatrix * point; + double sd = SdfFunctions::eval(func, invTransformedPoint); + return scale * sd + offset; +} + +Vec3d ProceduralMean::dmean_dp(const Point3d &point) const { + constexpr double eps = 0.001; + + std::array vals = { + mean(point + Vec3d(eps, 0., 0.)), + mean(point + Vec3d(0., eps, 0.)), + mean(point + Vec3d(0., 0., eps)), + mean(point)}; + + return Vec3d(vals[0] - vals[3], vals[1] - vals[3], vals[2] - vals[3]) / eps; +} + +SquaredExponentialCovariance::SquaredExponentialCovariance(const Json &json) { + sigma = getOptional(json, "sigma", 0.01); + lengthScale = getOptional(json, "lengthScale", 0.1); +} + +autodiff::real2nd SquaredExponentialCovariance::cov(const autodiff::real2nd &dis2) const { + return sqr(sigma) * exp(-(dis2 / (2 * sqr(lengthScale)))); +} +autodiff::dual2nd SquaredExponentialCovariance::cov(const autodiff::dual2nd &dis2) const { + return sqr(sigma) * exp(-(dis2 / (2 * sqr(lengthScale)))); +} +double SquaredExponentialCovariance::cov(double dis2) const { + return sqr(sigma) * fm::exp(-(dis2 / (2 * sqr(lengthScale)))); +} diff --git a/src/FunctionLayer/GaussianProcess/GPFunctions.h b/src/FunctionLayer/GaussianProcess/GPFunctions.h new file mode 100644 index 00000000..7c7e5c04 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GPFunctions.h @@ -0,0 +1,118 @@ +#pragma once +#include "CoreLayer/Adapter/JsonUtil.h" +#include "CoreLayer/Geometry/Geometry.h" +#include "CoreLayer/Geometry/Matrix.h" + +#include "FunctionLayer/SDFFunction/SdfFunctions.h" + +#include "autodiff/forward/dual.hpp" +#include "autodiff/forward/real.hpp" +#include "autodiff/forward/dual/eigen.hpp" +#include "autodiff/forward/real/eigen.hpp" + +enum class DerivativeType { + None, + First, +}; + +class MeanFunction { +public: + MeanFunction() = default; + virtual double operator()(const DerivativeType &derivativeType, const Point3d &point, const Vec3d &derivativeDir = {}) const { + if (derivativeType == DerivativeType::None) { + return mean(point); + } else { + return dot(dmean_dp(point), derivativeDir); + } + } + +protected: + virtual double mean(const Point3d &point) const = 0; + virtual Vec3d dmean_dp(const Point3d &point) const = 0; +}; + +class ProceduralMean : public MeanFunction { +public: + ProceduralMean(const Json &json); + +protected: + virtual double mean(const Point3d &point) const override; + virtual Vec3d dmean_dp(const Point3d &point) const override; + + mutable TransformMatrix3D transformMatrix; + mutable TransformMatrix3D invTransformMatrix; + + SdfFunctions::Function func; + double scale; + double offset; +}; + +class CovarianceFunction { +public: + CovarianceFunction() = default; + virtual double operator()(const DerivativeType &derivativeTypeX, const Point3d &pointX, + const DerivativeType &derivativeTypeY, const Point3d &pointY, + const Vec3d &derivativeDirX = {}, const Vec3d &derivativeDirY = {}) const { + if (derivativeTypeX == DerivativeType::None) { + if (derivativeTypeY == DerivativeType::None) { + return cov(pointX, pointY); + } else { + return dcov_dy(pointX, pointY, derivativeDirY); + } + } else { + if (derivativeTypeY == DerivativeType::None) { + return dcov_dx(pointX, pointY, derivativeDirX); + } else { + return dcov2_dxdy(pointX, pointY, derivativeDirX, derivativeDirY); + } + } + return 0.; + } + +protected: + // common interface + virtual double cov(const Point3d &pointX, const Point3d &pointY) const = 0; + // interfaces for autodiff + virtual autodiff::real2nd cov(const autodiff::Vector3real2nd &pointX, const autodiff::Vector3real2nd &pointY) const = 0; + virtual autodiff::dual2nd cov(const autodiff::Vector3dual2nd &pointX, const autodiff::Vector3dual2nd &pointY) const = 0; + + double dcov_dx(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirX) const; + double dcov_dy(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirY) const; + double dcov2_dxdy(const Point3d &pointX, const Point3d &pointY, const Vec3d &ddirX, const Vec3d &ddirY) const; +}; + +class StationaryCovariance : public CovarianceFunction { +protected: + virtual autodiff::real2nd cov(const autodiff::real2nd &dis2) const = 0; + virtual autodiff::dual2nd cov(const autodiff::dual2nd &dis2) const = 0; + virtual double cov(double dis) const = 0; + + virtual double cov(const Point3d &pointX, const Point3d &pointY) const { + auto d = pointX - pointY; + return cov(dot(d, d)); + } + virtual autodiff::real2nd cov(const autodiff::Vector3real2nd &pointX, const autodiff::Vector3real2nd &pointY) const { + auto d = pointX - pointY; + return cov(d.dot(d)); + } + virtual autodiff::dual2nd cov(const autodiff::Vector3dual2nd &pointX, const autodiff::Vector3dual2nd &pointY) const { + auto d = pointX - pointY; + return cov(d.dot(d)); + } +}; + +class SquaredExponentialCovariance : public StationaryCovariance { +public: + SquaredExponentialCovariance(const Json &json); + +protected: + virtual autodiff::real2nd cov(const autodiff::real2nd &dis2) const override; + virtual autodiff::dual2nd cov(const autodiff::dual2nd &dis2) const override; + virtual double cov(double dis2) const; + +protected: + double sigma; + double lengthScale; +}; + + diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp b/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp new file mode 100644 index 00000000..c0815e40 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp @@ -0,0 +1,321 @@ +#include "GaussianProcess.h" + +#include "GaussianProcessUtils.h" +#include "CoreLayer/Geometry/Frame.h" + +GPRealization::GPRealization(const GaussianProcess *_gp, + const Point3d *_points, + const DerivativeType *_derivativeTypes, + const Vec3d *_derivativeDirs, + const double *_values, + size_t numPoints, const Vec3d &derivativeDir) : gp(_gp) { + for (int i = 0; i < numPoints; ++i) { + Vec3d ddir = _derivativeDirs ? _derivativeDirs[i] : derivativeDir; + points.push_back(_points[i]); + derivativeTypes.push_back(_derivativeTypes[i]); + derivativeDirections.push_back(ddir); + values.push_back(_values[i]); + } +} + +void GPRealization::makeIntersection(size_t p, double offset) { + auto preV = values[p - 1]; + auto curV = values[p]; + + Point3d zeroCrossing = lerp(points[p - 1], points[p], offset); + double gradValue = (curV - preV) / ((points[p] - points[p - 1]).length()); + + points.push_back(zeroCrossing); + derivativeTypes.push_back(DerivativeType::None); + values.push_back(lerp(preV, curV, offset)); + // place holder + derivativeDirections.push_back({}); + + points.push_back(zeroCrossing); + derivativeTypes.push_back(DerivativeType::First); + values.push_back(gradValue); + derivativeDirections.push_back(normalize(points[p] - points[p - 1])); + + justIntersected = true; +} + +Vec3d GPRealization::sampleGradient(Point3d pos, Vec3d rayDir, Sampler &sampler) { + std::array gradPs{pos, pos, pos}; + std::array gradDerivs{DerivativeType::First, DerivativeType::First, DerivativeType::First}; + + Frame frame(rayDir); + + Vec3d sampleGrad = {}; + + // we just apply 'makeIntersect' method + if (justIntersected) { + std::array gradDirs{ + vec_conv(frame.s), + vec_conv(frame.t)}; + auto realization = gp->sampleCond(gradPs.data(), gradDerivs.data(), gradDirs.data(), gradDirs.size(), {}, + points.data(), derivativeTypes.data(), derivativeDirections.data(), values.data(), points.size(), {}, sampler); + // intersection's gradient is already known since we perform linear interpolation between points + sampleGrad = frame.toWorld({realization.values[0], realization.values[1], values[values.size() - 1]}); + + } else { + std::array gradDirs{ + vec_conv(frame.s), + vec_conv(frame.t), + vec_conv(frame.n)}; + auto realization = gp->sampleCond(gradPs.data(), gradDerivs.data(), gradDirs.data(), gradDirs.size(), {}, + points.data(), derivativeTypes.data(), derivativeDirections.data(), values.data(), points.size(), {}, sampler); + sampleGrad = frame.toWorld({realization.values[0], realization.values[1], realization.values[2]}); + } + return lastSampledGrad = sampleGrad; +} + +void GPRealization::applyMemoryModel(Vec3d rayDir, MemoryModel memoryModel) { + std::vector pointsNew; + std::vector derivativeDirectionsNew; + std::vector derivativeTypesNew; + std::vector valuesNew; + + size_t pointSize = points.size(); + switch (memoryModel) { + case MemoryModel::None: + break; + case MemoryModel::GlobalN: + // TODO(Cchen77): GlobalN memory model + break; + case MemoryModel::Renewal: { + size_t p = pointSize - 1; + if (justIntersected) { + p = pointSize - 2; + } + pointsNew.push_back(points[p]); + derivativeDirectionsNew.push_back(derivativeDirections[p]); + derivativeTypesNew.push_back(derivativeTypes[p]); + valuesNew.push_back(values[p]); + + points = pointsNew; + derivativeDirections = derivativeDirectionsNew; + derivativeTypes = derivativeTypesNew; + values = valuesNew; + break; + } + case MemoryModel::RenewalPlus: { + size_t p = pointSize - 1; + if (justIntersected) { + p = pointSize - 2; + } + pointsNew.push_back(points[p]); + derivativeDirectionsNew.push_back({}); + derivativeTypesNew.push_back(DerivativeType::None); + valuesNew.push_back(values[p]); + + pointsNew.push_back(points[p]); + derivativeDirectionsNew.push_back(rayDir); + derivativeTypesNew.push_back(DerivativeType::First); + valuesNew.push_back(dot(lastSampledGrad, rayDir)); + + points = pointsNew; + derivativeDirections = derivativeDirectionsNew; + derivativeTypes = derivativeTypesNew; + values = valuesNew; + break; + } + default: + break; + } +} + +GaussianProcess::GaussianProcess(std::shared_ptr _mean, std::shared_ptr _cov, const GPRealization &_globalCondition) : meanFunction(_mean), covFunction(_cov) { + initGlobalCondition(_globalCondition); +} + +void GaussianProcess::initGlobalCondition(const GPRealization &_globalCondition) { + globalCondition = _globalCondition; + + // solver picking strategy from https://github.com/daseyb/gpis-light-transport/blob/main/src/core/math/GaussianProcess.cpp + bool succesfullSolve = false; + auto kCC = covPriorSym(EXPAND_GPREALIZATION(globalCondition)); + if (kCC.rows() <= 64) { + Eigen::LDLT solver(kCC.triangularView()); + if (solver.info() == Eigen::ComputationInfo::Success && solver.isPositive()) { + globalCondtionSolver = solver; + succesfullSolve = true; + } + } + if (!succesfullSolve) { + Eigen::BDCSVD solver; + solver.compute(kCC.triangularView()); + if (solver.info() != Eigen::ComputationInfo::Success) { + globalCondtionSolver = solver; + succesfullSolve = true; + } + } + if (!succesfullSolve) { + std::cerr << "Fail to find a suitable solver for global condition!\n"; + } +} + +GPRealization GaussianProcess::sample(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir, Sampler &sampler) const { + std::vector values; + auto [_mean, _cov] = meanAndCov(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + MultiVariableNormalDistribution mvn(_mean, _cov); + auto v = mvn.sample(sampler); + for (int i = 0; i < numPoints; ++i) { + values.push_back(v(i)); + } + + return GPRealization(this, points, derivativeTypes, derivativeDirs, values.data(), numPoints, derivativeDir); +} + +GPRealization GaussianProcess::sampleCond(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir, + const Point3d *pointsCond, const DerivativeType *derivativeTypesCond, const Vec3d *derivativeDirsCond, const double *valuesCond, size_t numPointsCond, const Vec3d &derivativeDirCond, + Sampler &sampler) const { + if (numPointsCond == 0) { + return sample(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir, sampler); + } + + Eigen::MatrixXd kCC = covSym(pointsCond, derivativeTypesCond, derivativeDirsCond, numPointsCond, derivativeDirCond); + Eigen::MatrixXd kCx = cov(pointsCond, derivativeTypesCond, derivativeDirsCond, numPointsCond, derivativeDirCond, + points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + // we should solve kCC*x = kCx => xT kCC = kxC => xT = kxC*kCC^-1 + // solving strategy from https://github.com/daseyb/gpis-light-transport/blob/main/src/core/math/GaussianProcess.cpp + Eigen::MatrixXd solved; + bool succesfullSolve = false; + if (kCC.rows() <= 64) { + Eigen::LDLT solver(kCC.triangularView()); + if (solver.info() == Eigen::ComputationInfo::Success && solver.isPositive()) { + solved = solver.solve(kCx).transpose(); + if (solver.info() == Eigen::ComputationInfo::Success) { + succesfullSolve = true; + } else { + std::cerr << "Conditioning solving failed (LDLT)!\n"; + } + } + } + if (!succesfullSolve) { + Eigen::BDCSVD solver; + solver.compute(kCC.triangularView()); + + if (solver.info() != Eigen::ComputationInfo::Success) { + std::cerr << "Conditioning decomposition failed (BDCSVD)!\n"; + } + + solved = solver.solve(kCx).transpose(); + if (solver.info() != Eigen::ComputationInfo::Success) { + std::cerr << "Conditioning solving failed (BDCSVD)!\n"; + } + } + + auto [meanxx, kxx] = meanAndCov(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + + Eigen::VectorXd _mean = meanxx + + (solved * (Eigen::Map(valuesCond, numPointsCond) - mean(pointsCond, derivativeTypesCond, derivativeDirsCond, numPointsCond, derivativeDirCond))); + Eigen::MatrixXd _cov = kxx - + (solved * kCx); + MultiVariableNormalDistribution mvn(_mean, _cov); + auto v = mvn.sample(sampler); + std::vector values; + for (int i = 0; i < numPoints; ++i) { + values.push_back(v(i)); + } + return GPRealization(this, points, derivativeTypes, derivativeDirs, values.data(), numPoints, derivativeDir); +} + +Eigen::VectorXd GaussianProcess::meanPrior(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { + Eigen::VectorXd _mean(numPoints); + for (int i = 0; i < numPoints; ++i) { + Vec3d ddir_i = derivativeDirs ? derivativeDirs[i] : derivativeDir; + _mean(i) = (*meanFunction)(derivativeTypes[i], points[i], ddir_i); + } + return _mean; +} + +Eigen::MatrixXd GaussianProcess::covPrior(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, size_t numPointsX, const Vec3d &derivativeDirX, + const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, size_t numPointsY, const Vec3d &derivativeDirY) const { + Eigen::MatrixXd _cov(numPointsX, numPointsY); + for (int i = 0; i < numPointsX; ++i) { + Vec3d ddir_i = derivativeDirsX ? derivativeDirsX[i] : derivativeDirX; + for (int j = 0; j < numPointsY; ++j) { + Vec3d ddir_j = derivativeDirsY ? derivativeDirsY[j] : derivativeDirY; + _cov(i, j) = (*covFunction)(derivativeTypesX[i], pointsX[i], derivativeTypesY[j], pointsY[j], ddir_i, ddir_j); + } + } + return _cov; +} + +Eigen::MatrixXd GaussianProcess::covPriorSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { + Eigen::MatrixXd _cov(numPoints, numPoints); + for (int i = 0; i < numPoints; ++i) { + Vec3d ddir_i = derivativeDirs ? derivativeDirs[i] : derivativeDir; + for (int j = 0; j <= i; ++j) { + Vec3d ddir_j = derivativeDirs ? derivativeDirs[j] : derivativeDir; + double cov_i_j = (*covFunction)(derivativeTypes[i], points[i], derivativeTypes[j], points[j], ddir_i, ddir_j); + // so our kenrel should satisify cov(x,y) = cov(y,x) + _cov(i, j) = _cov(j, i) = cov_i_j; + } + } + return _cov; +} + +Eigen::VectorXd GaussianProcess::mean(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { + if (globalCondition.size() == 0) { + return meanPrior(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + } + auto mean_prior = meanPrior(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + + auto kCx = covPrior(EXPAND_GPREALIZATION(globalCondition), points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + auto solved = std::visit([&kCx](auto &solver) -> Eigen::MatrixXd { return solver.solve(kCx).transpose(); }, globalCondtionSolver); + return mean_prior + + (solved * (Eigen::Map(globalCondition.values.data(), globalCondition.size()) - meanPrior(EXPAND_GPREALIZATION(globalCondition)))); +} + +Eigen::MatrixXd GaussianProcess::cov(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, size_t numPointsX, const Vec3d &derivativeDirX, + const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, size_t numPointsY, const Vec3d &derivativeDirY) const { + if (globalCondition.size() == 0) { + return covPrior(pointsX, derivativeTypesX, derivativeDirsX, numPointsX, derivativeDirX, + pointsY, derivativeTypesY, derivativeDirsY, numPointsY, derivativeDirY); + } + auto cov_prior = covPrior(pointsX, derivativeTypesX, derivativeDirsX, numPointsX, derivativeDirX, + pointsY, derivativeTypesY, derivativeDirsY, numPointsY, derivativeDirY); + + auto kCx = covPrior(EXPAND_GPREALIZATION(globalCondition), pointsX, derivativeTypesX, derivativeDirsX, numPointsX, derivativeDirX); + auto kCy = covPrior(EXPAND_GPREALIZATION(globalCondition), pointsY, derivativeTypesY, derivativeDirsY, numPointsY, derivativeDirY); + auto solved = std::visit([&kCx](auto &solver) -> Eigen::MatrixXd { return solver.solve(kCx).transpose(); }, globalCondtionSolver); + return cov_prior - solved * kCy; +} + +Eigen::MatrixXd GaussianProcess::covSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { + if (globalCondition.size() == 0) { + return covPriorSym(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + } + auto cov_prior = covPriorSym(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + auto kCx = covPrior(EXPAND_GPREALIZATION(globalCondition), points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + auto solved = std::visit([&kCx](auto &solver) -> Eigen::MatrixXd { return solver.solve(kCx).transpose(); }, globalCondtionSolver); + return cov_prior - solved * kCx; +} + +std::tuple GaussianProcess::meanAndCov(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { + Eigen::VectorXd _mean(numPoints); + Eigen::MatrixXd _cov(numPoints, numPoints); + for (int i = 0; i < numPoints; ++i) { + Vec3d ddir_i = derivativeDirs ? derivativeDirs[i] : derivativeDir; + _mean(i) = (*meanFunction)(derivativeTypes[i], points[i], ddir_i); + for (int j = 0; j <= i; ++j) { + Vec3d ddir_j = derivativeDirs ? derivativeDirs[j] : derivativeDir; + double cov_i_j = (*covFunction)(derivativeTypes[i], points[i], derivativeTypes[j], points[j], ddir_i, ddir_j); + // so our kenrel should satisify cov(x,y) = cov(y,x) + _cov(i, j) = _cov(j, i) = cov_i_j; + } + } + if (globalCondition.size() == 0) { + return {_mean, _cov}; + } + auto kCx = covPrior(EXPAND_GPREALIZATION(globalCondition), points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + auto solved = std::visit([&kCx](auto &solver) -> Eigen::MatrixXd { return solver.solve(kCx).transpose(); }, globalCondtionSolver); + _mean += solved * (Eigen::Map(globalCondition.values.data(), globalCondition.size()) - meanPrior(EXPAND_GPREALIZATION(globalCondition))); + _cov -= solved * kCx; + return {_mean, _cov}; +} + +double GaussianProcess::goodStepSize(Point3d p, Vec3d rd, double desiredCov) const { + return 0.; +} diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcess.h b/src/FunctionLayer/GaussianProcess/GaussianProcess.h new file mode 100644 index 00000000..3befab43 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcess.h @@ -0,0 +1,83 @@ +#pragma once +#include "CoreLayer/Geometry/Geometry.h" +#include "CoreLayer/Adapter/JsonUtil.h" +#include "FunctionLayer/Sampler/Sampler.h" +#include "Eigen/Dense" + +#include "GPFunctions.h" + +#include +enum class MemoryModel { + None, + GlobalN, + Renewal, + RenewalPlus, +}; + +struct GaussianProcess; +struct GPRealization { + GPRealization() : gp(nullptr){}; + GPRealization(const GaussianProcess *_gp, const Point3d *_points, const DerivativeType *_derivativeTypes, const Vec3d *_derivativeDirs, const double *_values, size_t numPoints, const Vec3d &derivativeDir); + + std::vector points; + std::vector derivativeTypes; + std::vector derivativeDirections; + std::vector values; + + // there is a zero-crossing between values[p-1] ~ values[p] + virtual void makeIntersection(size_t p, double offset); + virtual Vec3d sampleGradient(Point3d pos, Vec3d rayDir, Sampler &sampler); + virtual void applyMemoryModel(Vec3d rayDir, MemoryModel memoryModel = MemoryModel::None); + + size_t size() const { return points.size(); } + bool isEmpty() const { return points.empty(); } + + bool justIntersected = false; + Vec3d lastSampledGrad; + const GaussianProcess *gp; + + void reset() { + points.clear(); + derivativeTypes.clear(); + derivativeDirections.clear(); + values.clear(); + + justIntersected = false; + lastSampledGrad = {}; + gp = nullptr; + } +}; +struct GaussianProcess { + GaussianProcess(std::shared_ptr _mean, std::shared_ptr _cov, const GPRealization &_globalCondition); + virtual GPRealization sample(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir, Sampler &sampler) const; + virtual GPRealization sampleCond(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir, + const Point3d *pointsCond, const DerivativeType *derivativeTypesCond, const Vec3d *derivativeDirsCond, const double *valuesCond, size_t numPointsCond, const Vec3d &derivativeDirCond, + Sampler &sampler) const; + std::shared_ptr meanFunction; + std::shared_ptr covFunction; + + virtual double goodStepSize(Point3d p, Vec3d rd, double desiredCov) const; + +protected: + GPRealization globalCondition; + std::variant, Eigen::BDCSVD> globalCondtionSolver; + void initGlobalCondition(const GPRealization &_globalCondition); + + // for mean and cov before apply global condition + Eigen::VectorXd meanPrior(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const; + Eigen::MatrixXd covPrior(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, size_t numPointsX, const Vec3d &derivativeDirX, + const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, size_t numPointsY, const Vec3d &derivativeDirY) const; + Eigen::MatrixXd covPriorSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPointsconst, const Vec3d &derivativeDir) const; + + // for mean and cov after apply global condition + Eigen::VectorXd mean(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const; + Eigen::MatrixXd cov(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, size_t numPointsX, const Vec3d &derivativeDirX, + const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, size_t numPointsY, const Vec3d &derivativeDirY) const; + Eigen::MatrixXd covSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPointsconst, const Vec3d &derivativeDir) const; + std::tuple meanAndCov(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const; +}; + +#define EXPAND_GPREALIZATION(real) \ + real.points.data(), real.derivativeTypes.data(), real.derivativeDirections.data(), real.size(), {} +#define EXPAND_GPREALIZATION_WITH_VALUE(real) \ + real.points.data(), real.derivativeTypes.data(), real.derivativeDirections.data(), real.values.data(), real.size(), {} \ No newline at end of file diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.cpp b/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.cpp new file mode 100644 index 00000000..b42b2a58 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.cpp @@ -0,0 +1,51 @@ +#include "GaussianProcessFactory.h" + +std::shared_ptr GaussianProcessFactory::LoadGaussianProcessFromJson(const Json &json) { + std::string gpType = json.at("type"); + if (gpType == "std") { + auto meanFunction = MeanFunctionFactory::LoadMeanFunctionFromJson(json["mean"]); + auto covarianceFunction = CovarianceFunctionFactory::LoadCovarianceFunctionFromJson(json["covariance"]); + GPRealization globalCondition; + if (json.find("global_condition") != json.end() && getOptional(json["global_condition"], "enable", false)) { + auto gcJson = json["global_condition"]; + auto gcPoints = gcJson["points"]; + auto gcDerivativeTypes = gcJson["derivative_types"]; + auto gcDerivativeDirs = gcJson["derivative_dirs"]; + auto gcValues = gcJson["values"]; + for (int i = 0; i < gcPoints.size(); ++i) { + globalCondition.points.push_back(gcPoints[i]); + if (gcDerivativeTypes[i] == "none") { + globalCondition.derivativeTypes.push_back(DerivativeType::None); + + } else if (gcDerivativeTypes[i] == "first") { + globalCondition.derivativeTypes.push_back(DerivativeType::First); + } + globalCondition.derivativeDirections.push_back(gcDerivativeDirs[i]); + globalCondition.values.push_back(gcValues[i]); + } + } + return std::make_shared(meanFunction, covarianceFunction, globalCondition); + } else if (gpType == "csg") { + // TODO(Cchen77): csg type process + return nullptr; + } + return nullptr; +} + +std::shared_ptr MeanFunctionFactory::LoadMeanFunctionFromJson(const Json &json) { + std::string meanType = json.at("type"); + if (meanType == "procedural") { + return std::make_shared(json); + } else { + return nullptr; + } +} + +std::shared_ptr CovarianceFunctionFactory::LoadCovarianceFunctionFromJson(const Json &json) { + std::string covType = json.at("type"); + if (covType == "squared_exponential") { + return std::make_shared(json); + } else { + return nullptr; + } +} diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.h b/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.h new file mode 100644 index 00000000..b1c2db01 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.h @@ -0,0 +1,12 @@ +#pragma once +#include"GaussianProcess.h" +#include"GPFunctions.h" +namespace GaussianProcessFactory { + std::shared_ptr LoadGaussianProcessFromJson(const Json& json); +}// namespace MediumFactory +namespace MeanFunctionFactory { + std::shared_ptr LoadMeanFunctionFromJson(const Json& json); +} +namespace CovarianceFunctionFactory { + std::shared_ptr LoadCovarianceFunctionFromJson(const Json &json); +} diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.cpp b/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.cpp new file mode 100644 index 00000000..de24d1e3 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.cpp @@ -0,0 +1,45 @@ +#include"GaussianProcessUtils.h" + +// https://github.com/daseyb/gpis-light-transport/blob/main/src/core/sampling/Gaussian.cpp +MultiVariableNormalDistribution::MultiVariableNormalDistribution(const Eigen::VectorXd &_mean, const Eigen::MatrixXd &_cov) { + + Eigen::LLT chol(_cov.triangularView()); + // We can only use the cholesky decomposition if + // the covariance matrix is symmetric, pos-definite. + // But a covariance matrix might be pos-semi-definite. + // In that case, we'll go to an EigenSolver + if (chol.info() == Eigen::Success) { + // Use cholesky solver + normTransform = chol.matrixL(); + } else { + Eigen::SelfAdjointEigenSolver eigs(_cov); + + if (eigs.info() != Eigen::ComputationInfo::Success) { + std::cerr << "Matrix square root failed!\n"; + } + + normTransform = eigs.eigenvectors() * eigs.eigenvalues().cwiseMax(0).cwiseSqrt().asDiagonal(); + } + + mean = _mean; +} + +// https://github.com/daseyb/gpis-light-transport/blob/main/src/core/sampling/Gaussian.cpp +Eigen::VectorXd MultiVariableNormalDistribution::sample(Sampler &sampler) const { + // Generate a vector of standard normal variates with the same dimension as the mean + Eigen::VectorXd z = Eigen::VectorXd(mean.size()); + + // We're always getting two samples, so make use of that + for (int i = 0; i < mean.size() / 2; i++) { + Vec2d norm_samp = rand_normal_2(sampler); + z(i * 2) = norm_samp.x; + z(i * 2 + 1) = norm_samp.y; + } + if (mean.size() % 2) + { + Vec2d norm_samp = rand_normal_2(sampler); + z(mean.size() - 1) = norm_samp.x; + } + Eigen::VectorXd sample = mean + normTransform * z; + return sample; +} diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.h b/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.h new file mode 100644 index 00000000..8df37b81 --- /dev/null +++ b/src/FunctionLayer/GaussianProcess/GaussianProcessUtils.h @@ -0,0 +1,50 @@ +#pragma once +#include"Eigen/Dense" +#include"Eigen/Sparse" +#include"CoreLayer/Geometry/Geometry.h" +#include"FunctionLayer/Sampler/Sampler.h" +#include"CoreLayer/Math/Common.h" + +template +inline To vec_conv(const From& vd) { + return To{vd.x, vd.y, vd.x}; +} +// Box muller transform +inline Vec2d rand_normal_2(Sampler &sampler) { + double u1 = sampler.sample1D(); + double u2 = sampler.sample1D(); + + double r = fm::sqrt(-2 * log(1. - u1)); + double x = fm::cos(2 * M_PI * u2); + double y = fm::sin(2 * M_PI * u2); + double z1 = r * x; + double z2 = r * y; + + return Vec2d(z1, z2); +} + +inline Eigen::MatrixXd project_to_psd(const Eigen::MatrixXd& in) { + + Eigen::SelfAdjointEigenSolver es(in); + Eigen::VectorXd eigenValues = es.eigenvalues(); + for (int i = 0; i < eigenValues.size(); ++i) { + if (eigenValues(i) < 0) { + eigenValues(i) = 0; + } + } + Eigen::MatrixXd psdMatrix = es.eigenvectors() * eigenValues.asDiagonal() * es.eigenvectors().transpose(); + return psdMatrix; + +} + +struct MultiVariableNormalDistribution { + Eigen::VectorXd mean; + + Eigen::BDCSVD svd; + + Eigen::MatrixXd normTransform; + + MultiVariableNormalDistribution(const Eigen::VectorXd &_mean, const Eigen::MatrixXd &_cov); + + Eigen::VectorXd sample(Sampler &sampler) const; +}; \ No newline at end of file diff --git a/src/FunctionLayer/Integrator/VolPathIntegrator.cpp b/src/FunctionLayer/Integrator/VolPathIntegrator.cpp index b617fc45..308d1f46 100644 --- a/src/FunctionLayer/Integrator/VolPathIntegrator.cpp +++ b/src/FunctionLayer/Integrator/VolPathIntegrator.cpp @@ -1,11 +1,11 @@ /** * @file VolPathIntegrator.cpp * @author Chenxi Zhou - * @brief + * @brief * @version 0.1 * @date 2022-09-22 - * - * @copyright NJUMeta (c) 2022 + * + * @copyright NJUMeta (c) 2022 * www.njumeta.com */ @@ -32,24 +32,36 @@ Spectrum VolPathIntegrator::Li(const Ray &initialRay, std::shared_ptr sce Ray ray = initialRay; - const double eps = 1e-4; + const double eps = 1e-5; int nBounces = 0; bool specularBounce = false; PathIntegratorLocalRecord prevLightSampleRecord; std::shared_ptr medium = nullptr; + // mainly for GPIS medium now,but maybe some other medium need it also. + MediumState mediumState{*sampler}; + auto itsOpt = scene->intersect(ray); while (true) { - MediumSampleRecord mRec; - if (medium && medium->sampleDistance(&mRec, ray, itsOpt.value(), sampler->sample2D())) { + MediumSampleRecord mRec{}; + mRec.mediumState = &mediumState; + if (medium && medium->sampleDistanceSafe(&mRec, ray, itsOpt, sampler->sample2D())) { // Handle medium distance sampling throughput *= mRec.tr * mRec.sigmaS / mRec.pdf; - Intersection mediumScatteringPoint = fulfillScatteringPoint(mRec.scatterPoint, ray.direction, medium); + + Intersection mediumScatteringPoint; + if (!mRec.needAniso) { + mediumScatteringPoint = fulfillScatteringPoint(mRec.scatterPoint, ray.direction, medium); + } else { + // abuse slightly for gpis medium + mediumScatteringPoint = fulfillScatteringPoint(mRec.scatterPoint, mRec.aniso, medium); + mediumScatteringPoint.geometryNormal = mRec.aniso; + } //* ----- Luminaire Sampling ----- for (int i = 0; i < nDirectLightSamples; ++i) { - PathIntegratorLocalRecord sampleLightRecord = sampleDirectLighting(scene, mediumScatteringPoint, ray); + PathIntegratorLocalRecord sampleLightRecord = sampleDirectLighting2(scene, mediumScatteringPoint, ray, &mediumState); PathIntegratorLocalRecord evalScatterRecord = evalScatter(mediumScatteringPoint, ray, sampleLightRecord.wi); if (!sampleLightRecord.f.isBlack()) { double misw = MISWeight(sampleLightRecord.pdf, evalScatterRecord.pdf); @@ -67,7 +79,7 @@ Spectrum VolPathIntegrator::Li(const Ray &initialRay, std::shared_ptr sce throughput *= sampleScatterRecord.f / sampleScatterRecord.pdf; ray = Ray{mediumScatteringPoint.position + sampleScatterRecord.wi * eps, sampleScatterRecord.wi}; itsOpt = scene->intersect(ray); - auto [sampleIts, tr] = intersectIgnoreSurface(scene, ray, medium); + auto [sampleIts, tr] = intersectIgnoreSurface2(scene, ray, medium, &mediumState); auto evalLightRecord = evalEmittance(scene, sampleIts, ray); if (!evalLightRecord.f.isBlack()) { double misw = MISWeight(sampleScatterRecord.pdf, evalLightRecord.pdf); @@ -82,24 +94,26 @@ Spectrum VolPathIntegrator::Li(const Ray &initialRay, std::shared_ptr sce PathIntegratorLocalRecord evalLightRecord = evalEmittance(scene, itsOpt, ray); if (nBounces == 0) { L += throughput * evalLightRecord.f; - if (!itsOpt) break; } } + // there will be case that itsOpt is null when light cross null interfaces and hit nothing, + // for that case we still need to calculate evnlight for it. + if (!itsOpt) break; auto its = itsOpt.value(); its.medium = medium; if (its.material->getBxDF(its)->isNull()) { medium = getTargetMedium(its, ray.direction); + mediumState.reset(); ray = Ray{its.position + eps * ray.direction, ray.direction}; itsOpt = scene->intersect(ray); - if (!itsOpt) break; continue; } //* Direct Illumination for (int i = 0; i < nDirectLightSamples; ++i) { - PathIntegratorLocalRecord sampleLightRecord = sampleDirectLighting(scene, its, ray); + PathIntegratorLocalRecord sampleLightRecord = sampleDirectLighting2(scene, its, ray,&mediumState); PathIntegratorLocalRecord evalScatterRecord = evalScatter(its, ray, sampleLightRecord.wi); if (!sampleLightRecord.f.isBlack()) { @@ -122,7 +136,7 @@ Spectrum VolPathIntegrator::Li(const Ray &initialRay, std::shared_ptr sce ray = Ray{its.position + sampleScatterRecord.wi * eps, sampleScatterRecord.wi}; itsOpt = scene->intersect(ray); - auto [sampleIts, tr] = intersectIgnoreSurface(scene, ray, medium); + auto [sampleIts, tr] = intersectIgnoreSurface2(scene, ray, medium, &mediumState); auto evalLightRecord = evalEmittance(scene, sampleIts, ray); if (!evalLightRecord.f.isBlack()) { @@ -337,9 +351,9 @@ std::shared_ptr VolPathIntegrator::getTargetMedium(const Intersection &i Spectrum VolPathIntegrator::evalTransmittance(std::shared_ptr scene, const Intersection &its, Point3d pointOnLight) const { - //if (!its.material) { - // std::cout << "Stop!\n"; - //} + // if (!its.material) { + // std::cout << "Stop!\n"; + // } float tmax = (pointOnLight - its.position).length(); Ray shadowRay{its.position, normalize(pointOnLight - its.position), 1e-4f, tmax - 1e-4f}; @@ -455,4 +469,126 @@ VolPathIntegrator::intersectIgnoreSurface(std::shared_ptr scene, testRayItsOpt = scene->intersect(marchRay); } return {testRayItsOpt, tr}; -} \ No newline at end of file +} + +/// @brief Iteratively eval the transmittance from intersection point its to direction wi. +/// @param scene Scene description used to operate ray intersection. +/// @param its Current intersection point or scattering point. +/// @param pointOnLight Sampled point on light. +/// @param meidumState Inital meidum state +/// @return The transmittance from pointOnLight to its. Transmittance will be zero if ray hits a non-null surface. +Spectrum VolPathIntegrator::evalTransmittance2(std::shared_ptr scene, const Intersection &its, Point3d pointOnLight, const MediumState *mediumState) const { + MediumState transientMeidumState = *mediumState; + float tmax = (pointOnLight - its.position).length(); + Ray shadowRay{its.position, normalize(pointOnLight - its.position), 1e-4f, tmax - 1e-4f}; + std::shared_ptr medium = its.medium; + Spectrum tr(1.f); + while (true) { + auto itsOpt = scene->intersect(shadowRay); + + if (medium) { + if (!itsOpt) { + tr *= medium->evalTransmittance2(shadowRay.origin, pointOnLight, &transientMeidumState); + break; + } + + if (!itsOpt->material->getBxDF(*itsOpt)->isNull()) { + tr = .0f; + break; + } + + tr *= medium->evalTransmittance2(shadowRay.origin, itsOpt->position, &transientMeidumState); + medium = getTargetMedium(*itsOpt, shadowRay.direction); + transientMeidumState.reset(); + shadowRay.origin = itsOpt->position; + shadowRay.timeMax -= itsOpt->t; + } else { + if (!itsOpt) break; + + if (!itsOpt->material->getBxDF(*itsOpt)->isNull()) { + tr = .0f; + break; + } + medium = getTargetMedium(*itsOpt, shadowRay.direction); + transientMeidumState.reset(); + shadowRay.origin = itsOpt->position; + shadowRay.timeMax -= itsOpt->t; + } + } + + return tr; +} +/// @brief Sample on the distribution of direct lighting and take medium transmittance into account. +/// @param scene Scene description. Multiple shadow ray intersect operations will be performed. +/// @param its Current intersection point which returned pdf dependent on. +/// @param ray Current ray. Should only be applied for time records. +/// @param meidumState Inital meidum state +/// @return Sampled direction on the distribution of direct lighting and corresponding solid angle dependent pdf. An extra flag indicites that whether it sampled on a delta distribution. +PathIntegratorLocalRecord VolPathIntegrator::sampleDirectLighting2(std::shared_ptr scene, const Intersection &its, const Ray &ray, const MediumState *mediumState) { + auto [light, pdfChooseLight] = chooseOneLight(scene, sampler->sample1D()); + auto record = light->sampleDirect(its, sampler->sample2D(), ray.timeMin); + double pdfDirect = record.pdfDirect * pdfChooseLight;// pdfScatter with respect to solid angle + Vec3d dirScatter = record.wi; + Point3d posL = record.dst; + Point3d posS = its.position; + auto transmittance = evalTransmittance2(scene, its, record.dst, mediumState); + // if (!its.material && transmittance.sum() < 2.9f) { + // std::cout << transmittance.sum() << "\n"; + // } + return {dirScatter, transmittance * record.s, pdfDirect, record.isDeltaPos}; +} +/// @brief Intersect in scene but ignore bsdf with isNull()==true. +/// @param scene Scene description where multiple intersect operation will be performed. +/// @param ray Initial ray. +/// @param medium Initial medium. +/// @param meidumState Inital meidum state +/// @return Intersected point and transmittance alone the way. +std::pair, Spectrum> VolPathIntegrator::intersectIgnoreSurface2(std::shared_ptr scene, const Ray &ray, std::shared_ptr medium, const MediumState *mediumState) const { + MediumState transientMeidumState = *mediumState; + + const double eps = 1e-5; + Vec3d dir = ray.direction; + + Spectrum tr(1.0); + Ray marchRay{ray.origin + dir * eps, dir}; + std::shared_ptr currentMedium = medium; + + Point3d lastScatteringPoint = ray.origin; + auto testRayItsOpt = scene->intersect(marchRay); + + // calculate the transmittance of last segment from lastScatteringPoint to testRayItsOpt. + while (true) { + + // corner case: infinite medium or infinite light source. + if (!testRayItsOpt.has_value()) { + if (currentMedium != nullptr) + tr = Spectrum(0.0); + return {testRayItsOpt, tr}; + } + + auto testRayIts = testRayItsOpt.value(); + + // corner case: non-null surface + if (testRayIts.material != nullptr) { + if (!testRayIts.material->getBxDF(testRayIts)->isNull()) { + if (currentMedium != nullptr) + tr *= currentMedium->evalTransmittance2(testRayIts.position, lastScatteringPoint, &transientMeidumState); + return {testRayItsOpt, tr}; + } + } + + // hit a null surface, calculate tr + if (currentMedium != nullptr) + tr *= currentMedium->evalTransmittance2(testRayIts.position, lastScatteringPoint, &transientMeidumState); + + // update medium + currentMedium = getTargetMedium(testRayIts, dir); + transientMeidumState.reset(); + + // update ray and intersection point. + marchRay.origin = testRayIts.position + dir * eps; + lastScatteringPoint = testRayIts.position; + testRayItsOpt = scene->intersect(marchRay); + } + return {testRayItsOpt, tr}; +} diff --git a/src/FunctionLayer/Integrator/VolPathIntegrator.h b/src/FunctionLayer/Integrator/VolPathIntegrator.h index 2184505b..5607fa57 100644 --- a/src/FunctionLayer/Integrator/VolPathIntegrator.h +++ b/src/FunctionLayer/Integrator/VolPathIntegrator.h @@ -4,8 +4,8 @@ * @brief Path Integrator with volume * @version 0.1 * @date 2022-09-22 - * - * @copyright NJUMeta (c) 2022 + * + * @copyright NJUMeta (c) 2022 * www.njumeta.com */ @@ -77,6 +77,23 @@ class VolPathIntegrator : public AbstractPathIntegrator { const Ray &ray, std::shared_ptr medium) const; + // some meidum need its state to eval transmittance,we provide the another version for that purpose + Spectrum evalTransmittance2(std::shared_ptr scene, + const Intersection &its, + Point3d pointOnLight, + const MediumState *mediumState) const; + + PathIntegratorLocalRecord sampleDirectLighting2(std::shared_ptr scene, + const Intersection &its, + const Ray &ray, + const MediumState *mediumState); + + std::pair, Spectrum> + intersectIgnoreSurface2(std::shared_ptr scene, + const Ray &ray, + std::shared_ptr medium, + const MediumState *mediumState) const; + protected: const int nPathLengthLimit = 64; const double pRussianRoulette = 0.95; diff --git a/src/FunctionLayer/Material/BxDF/LambertainBxDF.cpp b/src/FunctionLayer/Material/BxDF/LambertainBxDF.cpp index b1493f3d..5149ab41 100644 --- a/src/FunctionLayer/Material/BxDF/LambertainBxDF.cpp +++ b/src/FunctionLayer/Material/BxDF/LambertainBxDF.cpp @@ -33,6 +33,5 @@ BxDFSampleResult LambertainBxDF::sample(const Vec3d & wo, const Point2d &sample result.bxdfSampleType = BXDFType(BXDF_DIFFUSE | BXDF_REFLECTION); result.pdf = pdf(wo, wi); result.s = f(wo, wi); - return result; } diff --git a/src/FunctionLayer/Material/BxDF/MicrofacetDistribution.cpp b/src/FunctionLayer/Material/BxDF/MicrofacetDistribution.cpp index 1dafee8c..1df77525 100644 --- a/src/FunctionLayer/Material/BxDF/MicrofacetDistribution.cpp +++ b/src/FunctionLayer/Material/BxDF/MicrofacetDistribution.cpp @@ -184,8 +184,9 @@ double MicrograinDistribution::D(const Vec3d &wh, const Vec2d &alphaXY) const { double cos2 = cos * cos; double tan2 = (1.f - cos2) / cos2; double tmp = beta2 + tan2; - double num = beta2 * fm::log(1. - tau0) * fm::pow(1. - tau0, tan2 / tmp); + double num = beta2 * fm::log(1. - tau0) * fm::pow(1. - tau0, tan2 / tmp); float denum = tau0 * fm::pi_d * tmp * tmp * cos2 * cos2; + return -num / denum; } diff --git a/src/FunctionLayer/Material/BxDF/PorousLayerBxDF.cpp b/src/FunctionLayer/Material/BxDF/PorousLayerBxDF.cpp index 999372b5..eda54a2a 100644 --- a/src/FunctionLayer/Material/BxDF/PorousLayerBxDF.cpp +++ b/src/FunctionLayer/Material/BxDF/PorousLayerBxDF.cpp @@ -36,27 +36,32 @@ BxDFSampleResult PourousLayerBxDF::sample(const Vec3d &out, const Point2d &sampl double resample = sample[0] / sampleWeight; result = micrograinBRDF->sample(out, {resample, sample[1]}); } else { - result = bulkBxDF->sample(out, sample); + double resample = (1.-sample[0]) / (1.-sampleWeight); + result = bulkBxDF->sample(out, {resample, sample[1]}); chooseBulkLobe = true; } double w = getMicrograinWeight(tau0, beta, CosTheta(result.directionIn), CosTheta(out)); + if (chooseBulkLobe) { - result.s = (1 - w) * result.s + w*micrograinBRDF->f(out, result.directionIn); + Spectrum micrograin_f = micrograinBRDF->f(out, result.directionIn); + result.s = (1 - w) * result.s + w*micrograin_f; + double micrograin_pdf = micrograinBRDF->pdf(out, result.directionIn); + if (micrograin_pdf > 0) { + double misw = 1 - w; + result.s *= misw; + } + result.pdf *= 1.-sampleWeight; } else { - result.s = w * result.s + (1 - w) * bulkBxDF->f(out, result.directionIn); + Spectrum bulk_f = bulkBxDF->f(out, result.directionIn); + result.s = w * result.s + (1 - w) * bulk_f; + double bulk_pdf = bulkBxDF->pdf(out, result.directionIn); + if (bulk_pdf > 0) { + double misw = w; + result.s *= misw; + } + result.pdf *= sampleWeight; } - result.bxdfSampleType = BXDFType(BXDF_REFLECTION | BXDF_GLOSSY); - result.pdf = pdf(out, result.directionIn); + return result; - /*double theta = fm::acos(fm::sqrt(1 - sample[0])); - double phi = 2 * fm::pi_d * sample[1]; - Vec3d in = {fm::cos(phi) * fm::sin(theta), fm::sin(phi) * fm::sin(theta), fm::cos(theta)}; - - BxDFSampleResult result; - result.pdf = pdf(out, in); - result.s = f(out, in); - result.bxdfSampleType = BXDFType(BXDF_GLOSSY | BXDF_REFLECTION); - result.directionIn = in; - return result;*/ } diff --git a/src/FunctionLayer/Material/BxDF/PorousLayerMicrograinBxDF.cpp b/src/FunctionLayer/Material/BxDF/PorousLayerMicrograinBxDF.cpp index 60b7f9b4..e0d9d323 100644 --- a/src/FunctionLayer/Material/BxDF/PorousLayerMicrograinBxDF.cpp +++ b/src/FunctionLayer/Material/BxDF/PorousLayerMicrograinBxDF.cpp @@ -15,7 +15,7 @@ ConductorMicrograinBxDF::ConductorMicrograinBxDF(Spectrum _R0, Vec3d _k, double } Spectrum ConductorMicrograinBxDF::f(const Vec3d &out, const Vec3d &in) const { - if (CosTheta(out) < 1e-5 || CosTheta(in) < 1e-5) { + if (CosTheta(out) < 1e-3 || CosTheta(in) < 1e-3) { return 0.; } Vec3d eta = R0ToEta(R0); @@ -41,14 +41,14 @@ BxDFSampleResult ConductorMicrograinBxDF::sample(const Vec3d &out, const Point2d } Vec3d wh = distrib->Sample_wh(out, sample, {tau0, beta}); Vec3d in = Frame::reflect(out, wh); + result.pdf = pdf(out, in); + result.directionIn = in; + result.bxdfSampleType = BXDFType(BXDF_REFLECTION | BXDF_GLOSSY); if (dot(wh, out) < 0 || CosTheta(in) < 0) { - result.s = Spectrum(0); + result.s = Spectrum{0}; return result; } - result.pdf = pdf(out, in); result.s = f(out, in); - result.directionIn = in; - result.bxdfSampleType = BXDFType(BXDF_REFLECTION | BXDF_GLOSSY); return result; } @@ -63,7 +63,7 @@ Spectrum PlasticMicrograinBxDF::f(const Vec3d &out, const Vec3d &in) const { } Vec3d eta = R0ToEta(R0); Vec3d half = normalize(in + out); - if (dot(half, out) < 0) { + if (dot(half, out) <= 0) { return 0.; } Spectrum Fr = Fresnel::dielectricReflectance(eta, dot(half, in)); @@ -91,13 +91,13 @@ BxDFSampleResult PlasticMicrograinBxDF::sample(const Vec3d &out, const Point2d & } Vec3d wh = distrib->Sample_wh(out, sample, {tau0, beta}); Vec3d in = Frame::reflect(out, wh); + result.pdf = pdf(out, in); + result.directionIn = in; + result.bxdfSampleType = BXDFType(BXDF_REFLECTION | BXDF_GLOSSY); if (dot(wh, out) < 0 || CosTheta(in) < 0) { - result.s = Spectrum(0); + result.s = Spectrum{0}; return result; } - result.pdf = pdf(out, in); result.s = f(out, in); - result.directionIn = in; - result.bxdfSampleType = BXDFType(BXDF_REFLECTION | BXDF_GLOSSY); return result; } diff --git a/src/FunctionLayer/Medium/Beerslaw.cpp b/src/FunctionLayer/Medium/Beerslaw.cpp index f7f6d5f1..6cc7fe5f 100644 --- a/src/FunctionLayer/Medium/Beerslaw.cpp +++ b/src/FunctionLayer/Medium/Beerslaw.cpp @@ -1,22 +1,20 @@ #include "Beerslaw.h" bool BeerslawMedium::sampleDistance(MediumSampleRecord *mRec, - const Ray &ray, - const Intersection &its, - Point2d sample) const -{ + const Ray &ray, + const Intersection &its, + Point2d sample) const { //* No scattering, only absorbtion mRec->marchLength = its.t; mRec->pdf = 1; mRec->tr = evalTransmittance(ray.origin, its.position); - return false; + return false; } -Spectrum BeerslawMedium::evalTransmittance(Point3d from, - Point3d end) const -{ +Spectrum BeerslawMedium::evalTransmittance(Point3d from, + Point3d end) const { double distance = (end - from).length(); return exp(mDensity * -distance); -} \ No newline at end of file +} \ No newline at end of file diff --git a/src/FunctionLayer/Medium/Beerslaw.h b/src/FunctionLayer/Medium/Beerslaw.h index 470f8d0b..607ed2cb 100644 --- a/src/FunctionLayer/Medium/Beerslaw.h +++ b/src/FunctionLayer/Medium/Beerslaw.h @@ -25,6 +25,7 @@ class BeerslawMedium : public Medium { Point2d sample) const; virtual Spectrum evalTransmittance(Point3d from, Point3d dest) const; + private: Spectrum mDensity; }; \ No newline at end of file diff --git a/src/FunctionLayer/Medium/GPISMedium.cpp b/src/FunctionLayer/Medium/GPISMedium.cpp new file mode 100644 index 00000000..c18b7bbb --- /dev/null +++ b/src/FunctionLayer/Medium/GPISMedium.cpp @@ -0,0 +1,148 @@ +#include "GPISMedium.h" +#include "GPISPhase.h" +#include "FunctionLayer/Sampler/Independent.h" +#include "FunctionLayer/GaussianProcess/GaussianProcessFactory.h" + +GPISMedium::GPISMedium(const Json &json) : Medium(std::make_shared(json["phase"])) { + marchingNumSamplePoints = getOptional(json, "marching_num_sample_points", 8); + if (marchingNumSamplePoints < 2) marchingNumSamplePoints = 2; + marchingStepSize = getOptional(json, "marching_step_size", 0.); + marchingDesiredCov = getOptional(json, "marching_desired_cov", 0.); + + gaussianProcess = GaussianProcessFactory::LoadGaussianProcessFromJson(json["gaussian_process"]); +} + +bool GPISMedium::sampleDistance(MediumSampleRecord *mRec, const Ray &ray, const Intersection &its, Point2d sample) const { +#if defined(ENABLE_GPISMEDIUM) + const double eps = 1e-6; + GPRealization &gpRealization = mRec->mediumState->realization; + Sampler &sampler = mRec->mediumState->sampler; + + Ray r = ray; + r.timeMax = std::min(r.timeMax, r.timeMin + 200); + r.timeMax = std::min(its.t, r.timeMax); + + double t = r.timeMin; + bool intersected = false; + do { + intersected = intersectGP(r, gpRealization, t, sampler); + if (t < r.timeMax) { + Point3d point = r.origin + t * r.direction; + Vec3d grad = gpRealization.sampleGradient(point, r.direction, sampler); + if (intersected) { + mRec->aniso = normalize(grad); + mRec->marchLength = t; + mRec->scatterPoint = point; + } + gpRealization.applyMemoryModel(r.direction, memoryModel); + } + } while (!intersected && r.timeMax - t > eps); + mRec->sigmaS = 1.; + mRec->sigmaA = 0.; + mRec->pdf = 1.; + mRec->tr = 1.; + mRec->needAniso = true; + return intersected; +#else + return false; +#endif +} + +Spectrum GPISMedium::evalTransmittance(Point3d from, Point3d dest) const { + // TODO(Cchen77): + // limited by vol path tracer framework,we just have a naive version,which just like applying Renewal memory model before cast the shadowray + // need a more elegant way + IndependentSampler transientSampler(1, 5);// 1,5 is meaningless + + Vec3d direction = (dest - from); + if (direction.isZero()) { + return 1.; + } + direction = normalize(direction); + Ray ray{from, direction}; + MediumState mediumState{transientSampler}; + mediumState.realization.reset(); + mediumState.realization.points.push_back(dest); + mediumState.realization.values.push_back(0); + mediumState.realization.derivativeTypes.push_back(DerivativeType::None); + mediumState.realization.derivativeDirections.push_back({}); + + MediumSampleRecord sampleRecord{}; + sampleRecord.mediumState = &mediumState; + + Intersection its; + its.t = (dest - from)[0] / direction[0]; + bool shadowed = sampleDistance(&sampleRecord, ray, its, {}); + + return 1 - shadowed; +} + +Spectrum GPISMedium::evalTransmittance2(Point3d from, Point3d dest, MediumState *mediumState) const { + Vec3d direction = (dest - from); + if (direction.length() < 1e-4) { + return 1.; + } + direction = normalize(direction); + Ray ray{from, direction}; + + MediumSampleRecord sampleRecord{}; + sampleRecord.mediumState = mediumState; + + Intersection its; + its.t = (dest - from).length(); + bool shadowed = sampleDistance(&sampleRecord, ray, its, {}); + + return 1. - shadowed; +} + +bool GPISMedium::intersectGP(const Ray &ray, GPRealization &gpRealization, double &t, Sampler &sampler) const { + double maxDistance = ray.timeMax - t; + double determinedStepSize = maxDistance / (marchingNumSamplePoints - 1); + if (marchingStepSize < determinedStepSize) { + determinedStepSize = marchingStepSize; + } + + std::vector points; + std::vector derivativeTypes; + std::vector ts; + + Point3d p = ray.origin + ray.direction * (determinedStepSize * 0.1 + t); + points.push_back(p); + derivativeTypes.push_back(DerivativeType::None); + ts.push_back(determinedStepSize * 0.1 + t); + + for (int i = 1; i < marchingNumSamplePoints; ++i) { + Point3d p = ray.origin + ray.direction * (i * determinedStepSize + t); + points.push_back(p); + derivativeTypes.push_back(DerivativeType::None); + ts.push_back(i * determinedStepSize + t); + } + + // it's the first time we have a realization,so we can sampling without condition + if (gpRealization.isEmpty()) { + gpRealization = gaussianProcess->sample(points.data(), derivativeTypes.data(), nullptr, marchingNumSamplePoints, {}, sampler); + } + // use last realization as conditon + else { + gpRealization = gaussianProcess->sampleCond( + points.data(), derivativeTypes.data(), nullptr, marchingNumSamplePoints, {}, EXPAND_GPREALIZATION_WITH_VALUE(gpRealization), sampler); + } + + double lastV = gpRealization.values[0]; + double lastT = ray.timeMin; + t = ts[0]; + for (int i = 1; i < marchingNumSamplePoints; ++i) { + double curV = gpRealization.values[i]; + double curT = ts[i]; + if (lastV * curV <= 0.) { + double offset = lastV / (lastV - curV); + gpRealization.makeIntersection(i, offset); + t = lerp(lastT, curT, offset); + return true; + } + t = ts[i]; + lastV = curV; + lastT = curT; + } + return false; +} diff --git a/src/FunctionLayer/Medium/GPISMedium.h b/src/FunctionLayer/Medium/GPISMedium.h new file mode 100644 index 00000000..ac41bf38 --- /dev/null +++ b/src/FunctionLayer/Medium/GPISMedium.h @@ -0,0 +1,35 @@ +#pragma once + +#include "Medium.h" +#include "CoreLayer/Adapter/JsonUtil.h" +#include "FunctionLayer/GaussianProcess/GaussianProcess.h" + +class GPISMedium : public Medium { +public: + GPISMedium() = default; + GPISMedium(const Json &json); + + virtual bool sampleDistance(MediumSampleRecord *mRec, + const Ray &ray, + const Intersection &its, + Point2d sample) const override; + + virtual Spectrum evalTransmittance(Point3d from, + Point3d dest) const override; + + virtual Spectrum evalTransmittance2(Point3d from, + Point3d dest, + MediumState *mediumState) const override; + + bool intersectGP(const Ray &ray, GPRealization &gpRealization, double &t, Sampler &sampler) const; + +private: + std::shared_ptr gaussianProcess; + + int marchingNumSamplePoints; + double marchingStepSize; + double marchingDesiredCov; + + MemoryModel memoryModel = MemoryModel::RenewalPlus; + +}; \ No newline at end of file diff --git a/src/FunctionLayer/Medium/GPISPhase.cpp b/src/FunctionLayer/Medium/GPISPhase.cpp new file mode 100644 index 00000000..fbadb6a1 --- /dev/null +++ b/src/FunctionLayer/Medium/GPISPhase.cpp @@ -0,0 +1,21 @@ +#include "GPISPhase.h" +#include "FunctionLayer/Material/MaterialFactory.h" +#include "FunctionLayer/Intersection.h" + +GPISPhase::GPISPhase(const Json &json) { + auto transientMaterial = MaterialFactory::LoadMaterialFromJson(json); + Intersection pseudoIts{}; + innerBxDF = transientMaterial->getBxDF(pseudoIts); +} + +std::tuple GPISPhase::evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { + + auto phaseValue = innerBxDF->f(wo, wi) * std::max(0., CosTheta(wi)); + double phasePdf = innerBxDF->pdf(wo, wi); + return {phaseValue, phasePdf, false}; +} + +std::tuple GPISPhase::samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const { + auto bxdfSampleResult = innerBxDF->sample(wo, sample); + return {bxdfSampleResult.directionIn, bxdfSampleResult.s * std::max(0., CosTheta(bxdfSampleResult.directionIn)), bxdfSampleResult.pdf, BxDF::MatchFlags(bxdfSampleResult.bxdfSampleType, BXDF_SPECULAR)}; +} diff --git a/src/FunctionLayer/Medium/GPISPhase.h b/src/FunctionLayer/Medium/GPISPhase.h new file mode 100644 index 00000000..ae52b4e3 --- /dev/null +++ b/src/FunctionLayer/Medium/GPISPhase.h @@ -0,0 +1,25 @@ +#pragma once + +#include "Phase.h" +#include "CoreLayer/Adapter/JsonUtil.h" +#include "FunctionLayer/Material/Material.h" +#include "FunctionLayer/Material/BxDF/BxDF.h" + +// MARK(Cchen77): +// we treat GPIS as a special medium and its phase function can be any bxdf. +// our implmentation now is a little bit tricky now,we just utilize the material factory to produce material which then produce bxdf for evaling&sampling +// since we have not suitable mapping system for meidum,do not use materials with texture. +class GPISPhase : public PhaseFunction { +public: + GPISPhase() = default; + GPISPhase(const Json &json); + + virtual std::tuple evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const override; + + virtual std::tuple + samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const override; + +protected: + + std::shared_ptr innerBxDF; +}; \ No newline at end of file diff --git a/src/FunctionLayer/Medium/HGPhase.cpp b/src/FunctionLayer/Medium/HGPhase.cpp index 9d4690ff..cf4b1a42 100644 --- a/src/FunctionLayer/Medium/HGPhase.cpp +++ b/src/FunctionLayer/Medium/HGPhase.cpp @@ -1,7 +1,7 @@ #include "HGPhase.h" #include"CoreLayer/ColorSpace/Color.h" -std::tuple +std::tuple HGPhase::evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { Vec3d w = wo * wi; float cosTheta = (w[0] + w[1] + w[2]) / wo.length() * wi.length(); @@ -11,7 +11,7 @@ HGPhase::evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { return {phaseValue,phasePdf, false}; } -std::tuple +std::tuple HGPhase::samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const { float cosTheta; if (std::abs(g) < 1e-3) { diff --git a/src/FunctionLayer/Medium/HGPhase.h b/src/FunctionLayer/Medium/HGPhase.h index 3b1ea881..6db9e9f3 100644 --- a/src/FunctionLayer/Medium/HGPhase.h +++ b/src/FunctionLayer/Medium/HGPhase.h @@ -16,10 +16,10 @@ class HGPhase : public PhaseFunction { public: HGPhase(float g) : g(g) {} - virtual std::tuple + virtual std::tuple evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const override; - virtual std::tuple + virtual std::tuple samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const override; public: diff --git a/src/FunctionLayer/Medium/Heterogeneous.cpp b/src/FunctionLayer/Medium/Heterogeneous.cpp index 7ef158f9..76f76ff7 100644 --- a/src/FunctionLayer/Medium/Heterogeneous.cpp +++ b/src/FunctionLayer/Medium/Heterogeneous.cpp @@ -213,10 +213,10 @@ bool HeterogeneousMedium::sampleDistance(MediumSampleRecord *mRec, const Ray &ra Point3d idx = worldToIndex(mRec->scatterPoint); mRec->marchLength = t_world; - mRec->sigmaA = Spectrum(.0f);//TODO + mRec->sigmaA = Spectrum(.0f);// TODO mRec->sigmaS = Spectrum(density); mRec->tr = Spectrum(fm::exp(-thick)); - mRec->pdf = mRec->tr[0] * density;//TODO + mRec->pdf = mRec->tr[0] * density;// TODO return true; } diff --git a/src/FunctionLayer/Medium/IsotropicPhase.cpp b/src/FunctionLayer/Medium/IsotropicPhase.cpp index e9039a7d..431223a6 100644 --- a/src/FunctionLayer/Medium/IsotropicPhase.cpp +++ b/src/FunctionLayer/Medium/IsotropicPhase.cpp @@ -2,13 +2,13 @@ #include "CoreLayer/Math/Warp.h" -std::tuple +std::tuple IsotropicPhase::evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { return {0.25 * INV_PI, 0.25 * INV_PI, false}; } -std::tuple +std::tuple IsotropicPhase::samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const { return {SquareToUniformSphere(sample), 0.25 * INV_PI, 0.25 * INV_PI, false}; diff --git a/src/FunctionLayer/Medium/IsotropicPhase.h b/src/FunctionLayer/Medium/IsotropicPhase.h index 6178cd4c..19472908 100644 --- a/src/FunctionLayer/Medium/IsotropicPhase.h +++ b/src/FunctionLayer/Medium/IsotropicPhase.h @@ -16,10 +16,10 @@ class IsotropicPhase : public PhaseFunction { public: IsotropicPhase() = default; - virtual std::tuple + virtual std::tuple evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const override; - virtual std::tuple + virtual std::tuple samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const override; }; \ No newline at end of file diff --git a/src/FunctionLayer/Medium/Medium.h b/src/FunctionLayer/Medium/Medium.h index 339fee02..4fba33af 100644 --- a/src/FunctionLayer/Medium/Medium.h +++ b/src/FunctionLayer/Medium/Medium.h @@ -6,7 +6,7 @@ * @version 0.1 * @date 2022-04-30 * - * @copyright NJUMeta (c) 2022 + * @copyright NJUMeta (c) 2022 * www.njumeta.com * */ @@ -18,47 +18,88 @@ #include "CoreLayer/Ray/Ray.h" #include "FunctionLayer/Intersection.h" +#include "FunctionLayer/GaussianProcess/GaussianProcess.h" + +#include + +struct MediumState { + // mainly for gpis medium which need the sampler + Sampler &sampler; +#if defined(ENABLE_GPISMEDIUM) + GPRealization realization; +#endif + void reset() { +#if defined(ENABLE_GPISMEDIUM) + realization.reset(); +#endif + } +}; + struct MediumSampleRecord { double marchLength; - double pdf; + double pdf; Point3d scatterPoint; Vec3d wi; Spectrum tr; Spectrum sigmaA; Spectrum sigmaS; -}; + Spectrum emission; + Vec3d aniso; + bool needAniso; + + MediumState *mediumState; +}; class Medium { public: + /** + * @brief Sample a distance, the ray will transport in media without collision until reach the distance + * @return + * - true, if the ray will endure a collision in medium + * - false, if the ray will pass through the media without collision + */ - /** - * @brief Sample a distance, the ray will transport in media without collision until reach the distance - * @return - * - true, if the ray will endure a collision in medium - * - false, if the ray will pass through the media without collision - */ + Medium(std::shared_ptr phase) : mPhase(phase) {} - Medium(std::shared_ptr phase) : mPhase(phase) { } + virtual bool sampleDistance(MediumSampleRecord *mRec, + const Ray &ray, + const Intersection &its, + Point2d sample) const = 0; - virtual bool sampleDistance(MediumSampleRecord *mRec, - const Ray &ray, - const Intersection &itsOpt, - Point2d sample) const = 0; + bool sampleDistanceSafe(MediumSampleRecord *mRec, + const Ray &ray, + const std::optional &itsOpt, + Point2d sample) { - virtual Spectrum evalTransmittance (Point3d from, Point3d dest) const = 0; + // this should not happen actually since our medium should be in a bounding box. + // but with testing,we found that sometimes ray can not in and out the medium and hit the boundary correctly because of the accuracy + // we have already chosen a small 'eps' to avoid the problem,but we don't know if that will happen again. + // provided a safe version so that the program won't crash + if (!itsOpt) { + mRec->tr = 1.; + mRec->pdf = 1.; + return false; + } + return sampleDistance(mRec, ray, itsOpt.value(), sample); + } - auto evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { - return mPhase->evalPhase(wo, wi, scatterPoint); - } + virtual Spectrum evalTransmittance(Point3d from, Point3d dest) const = 0; - auto samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const { - return mPhase->samplePhase(wo, scatterPoint, sample); - } + // interface for medium those who need state. + virtual Spectrum evalTransmittance2(Point3d from, Point3d dest, MediumState *mediumState) const { + return evalTransmittance(from, dest); + } -protected: + auto evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const { + return mPhase->evalPhase(wo, wi, scatterPoint); + } - std::shared_ptr mPhase; + auto samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const { + return mPhase->samplePhase(wo, scatterPoint, sample); + } +protected: + std::shared_ptr mPhase; }; \ No newline at end of file diff --git a/src/FunctionLayer/Medium/MediumFactory.cpp b/src/FunctionLayer/Medium/MediumFactory.cpp index 0f608f35..5cce3063 100644 --- a/src/FunctionLayer/Medium/MediumFactory.cpp +++ b/src/FunctionLayer/Medium/MediumFactory.cpp @@ -3,6 +3,7 @@ #include "Heterogeneous.h" #include "IsotropicPhase.h" #include "Beerslaw.h" +#include "GPISMedium.h" namespace MediumFactory { std::shared_ptr LoadMediumFromJson(const Json json) { @@ -31,6 +32,8 @@ std::shared_ptr LoadMediumFromJson(const Json json) { EulerType::EULER_YZX); return std::make_shared(filepath, std::make_shared(), transformMatrix, sigmaScale); + } else if (medium_type == "gpis") { + return std::make_shared(json); } return nullptr; } diff --git a/src/FunctionLayer/Medium/Phase.h b/src/FunctionLayer/Medium/Phase.h index 7078e25a..07642410 100644 --- a/src/FunctionLayer/Medium/Phase.h +++ b/src/FunctionLayer/Medium/Phase.h @@ -4,8 +4,8 @@ * @brief Abstraction for phase function * @version 0.1 * @date 2022-09-25 - * - * @copyright NJUMeta (c) 2022 + * + * @copyright NJUMeta (c) 2022 * www.njumeta.com */ @@ -13,34 +13,34 @@ #include #include "CoreLayer/Geometry/Geometry.h" +#include "CoreLayer/ColorSpace/Color.h" class PhaseFunction { public: PhaseFunction() = default; /** - * @brief Given the direction of wi/wo (both in world), + * @brief Given the direction of wi/wo (both in world), * evaluate the phase function - * + * * @param wo The direction towards light/next intersection * @param wi The direction towards camera/previous intersection * @param scatterPoint The position where scatter occurs - * @return std::tuple + * @return std::tuple * phaseValue, phasePdf and whether delta distribution */ - virtual std::tuple + virtual std::tuple evalPhase(Vec3d wo, Vec3d wi, Point3d scatterPoint) const = 0; /** * @brief Given the direction of wi (in world), * sample a scatter direction - * - * @param wo - * @param scatterPoint + * + * @param wo + * @param scatterPoint * @return std::tuple - * wi, phaseValue, phasePdf and whether delta distribution + * wi, phaseValue, phasePdf and whether delta distribution */ - virtual std::tuple + virtual std::tuple samplePhase(Vec3d wo, Point3d scatterPoint, Point2d sample) const = 0; - }; \ No newline at end of file diff --git a/src/FunctionLayer/SDFFunction/SdfFunctions.cpp b/src/FunctionLayer/SDFFunction/SdfFunctions.cpp new file mode 100644 index 00000000..bca3b1e1 --- /dev/null +++ b/src/FunctionLayer/SDFFunction/SdfFunctions.cpp @@ -0,0 +1,137 @@ +#include "SdfFunctions.h" +#include "Eigen/Dense" +#include "FastMath/FastMath.h" +#include "CoreLayer/Math/Common.h" +/* + Copyright 2020 Towaki Takikawa @yongyuanxi + The MIT License + Link: N/A + */ + +/****************************************************************************** + * The MIT License (MIT) + * Copyright (c) 2021, NVIDIA CORPORATION. + * Permission is hereby granted, free of charge, to any person obtaining a copy of + * this software and associated documentation files (the "Software"), to deal in + * the Software without restriction, including without limitation the rights to + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of + * the Software, and to permit persons to whom the Software is furnished to do so, + * subject to the following conditions: + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR + * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + * IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + ******************************************************************************/ + +///////////////////////////////////////////////////////////////////////////////////////////////////////////// +// distance functions +// taken from https://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm +///////////////////////////////////////////////////////////////////////////////////////////////////////////// + +double sdSphere(Eigen::Vector3d v, double r) { + return v.norm() - r; +} + +double sdTorus(Eigen::Vector3d p, Eigen::Vector2d t) { + Eigen::Vector2d q = Eigen::Vector2d(fm::sqrt(p.x() * p.x() + p.z() * p.z()) - t.x(), p.y()); + return q.norm() - t.y(); +} + +double sdCone(Eigen::Vector3d p, Eigen::Vector2d c) { + double q = (p.head<2>().norm()); + return c.dot(Eigen::Vector2d(q, p.z())); +} + +double sdCappedCylinder(Eigen::Vector3d p, double h, double r) { + Eigen::Vector2d d = (Eigen::Vector2d(p.head<2>().norm(), p.y()).cwiseAbs() - Eigen::Vector2d(h, r)); + return std::min(std::max(d.x(), d.y()), 0.) + d.cwiseMax(Eigen::Vector2d(0.,0.)).norm(); +} + +double sdTriPrism(Eigen::Vector3d p, Eigen::Vector2d h) { + Eigen::Vector3d q = p.cwiseAbs(); + return std::max(q.z() - h.y(), std::max(q.x() * 0.866025 + p.y() * 0.5, -p.y()) - h.x() * 0.5); +} + +double opSmoothUnion(double d1, double d2, double k) { + double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0); + return lerp(d2, d1, h) - k * h * (1.0 - h); +} + +double ssub(double d1, double d2, double k) { + double h = std::clamp(0.5 - 0.5 * (d2 + d1) / k, 0.0, 1.0); + return lerp(d2, -d1, h) + k * h * (1.0 - h); +} + +///////////////////////////////////////////////////////////////////////////////////////////////////////////// +// actual distance functions +///////////////////////////////////////////////////////////////////////////////////////////////////////////// + +double sdBase(Eigen::Vector3d p) { + + double base = opSmoothUnion(sdCone((Eigen::AngleAxisd(-90. * fm::pi_d / 180.0, Eigen::Vector3d(1.0, 0.0, 0.0))).toRotationMatrix() * (p + Eigen::Vector3d(0.0, 0.9, 0.0)), + Eigen::Vector2d(fm::pi_d / 3.0, fm::pi_d / 3.0)), + sdCone((Eigen::AngleAxisd(90. * fm::pi_d / 180.0, Eigen::Vector3d(1.0, 0.0, 0.0))).toRotationMatrix() * (p - Eigen::Vector3d(0.0, 0.9, 0.0)), + Eigen::Vector2d(fm::pi_d / 3.0, fm::pi_d / 3.0)), + 0.02); + base = std::max(base, sdCappedCylinder(p, 1.1, 0.25)) * 0.7; + base = std::max(-sdCappedCylinder(p, 0.6, 0.3), base); + base = std::max(-sdTriPrism((Eigen::AngleAxisd(90. * fm::pi_d / 180.0, Eigen::Vector3d(1.0, 0.0, 0.0))).toRotationMatrix() * (p + Eigen::Vector3d(0.0, 0.0, -1.0)), Eigen::Vector2d(1.2, 0.3)), base); + return base; +} + +double sdKnob(Eigen::Vector3d p) { + double sphere = sdSphere(p, 1.0); + double cutout = sdSphere(p - Eigen::Vector3d(0.0, 0.5, 0.5), 0.7); + double cutout_etch = sdTorus((Eigen::AngleAxisd(-45 * fm::pi_d / 180.0, Eigen::Vector3d(1.0, 0.0, 0.0))).toRotationMatrix() * (p - Eigen::Vector3d(0.0, 0.2, 0.2)), Eigen::Vector2d(1.0, 0.05)); + double innersphere = sdSphere(p - Eigen::Vector3d(0.0, 0.0, 0.0), 0.75); + + double d = ssub(cutout, sphere, 0.1); + d = std::min(d, innersphere); + d = std::max(-cutout_etch, d); + d = std::min(ssub(sphere, sdBase(p - Eigen::Vector3d(0.0, -0.775, 0.0)), 0.1), d); + return d; +} + +double sdKnobInner(Eigen::Vector3d p) { + return sdSphere(p - Eigen::Vector3d(0.0, 0.0, 0.0), 0.75); +} + +double sdKnobOuter(Eigen::Vector3d p) { + double sphere = sdSphere(p, 1.0); + double cutout = sdSphere(p - Eigen::Vector3d(0.0, 0.5, 0.5), 0.7); + double cutout_etch = sdTorus((Eigen::AngleAxisd(-45. * fm::pi_d / 180.0, Eigen::Vector3d(1.0, 0.0, 0.0))).toRotationMatrix() * (p - Eigen::Vector3d(0.0, 0.2, 0.2)), Eigen::Vector2d(1.0, 0.05)); + double innersphere = sdSphere(p - Eigen::Vector3d(0.0, 0.0, 0.0), 0.75); + + double d = ssub(cutout, sphere, 0.1); + d = std::max(d, -innersphere); + d = std::max(-cutout_etch, d); + d = std::min(ssub(sphere, sdBase(p - Eigen::Vector3d(0.0, -0.775, 0.0)), 0.1), d); + return d; +} + +double SdfFunctions::knob(Point3d p) { + const double scale = 0.8; + p *= 1.0 / scale; + return sdKnob({p.x, p.y, p.z}) * scale; +} + +double SdfFunctions::knob_inner(Point3d p) { + const double scale = 0.8; + p *= 1.0 / scale; + return sdKnobInner({p.x, p.y, p.z}) * scale; +} + +double SdfFunctions::knob_outer(Point3d p) { + const float scale = 0.8; + p *= 1. / scale; + return sdKnobOuter({p.x, p.y, p.z}) * scale; +} + +double SdfFunctions::two_spheres(Point3d p) { + Eigen::Vector3d eigen_p = {p.x, p.y, p.z}; + return std::min((eigen_p - Eigen::Vector3d(0., 10., 0.)).norm() - 9.5, (eigen_p - Eigen::Vector3d(0., -10., 0.)).norm() - 9.5); +} \ No newline at end of file diff --git a/src/FunctionLayer/SDFFunction/SdfFunctions.h b/src/FunctionLayer/SDFFunction/SdfFunctions.h new file mode 100644 index 00000000..1f920911 --- /dev/null +++ b/src/FunctionLayer/SDFFunction/SdfFunctions.h @@ -0,0 +1,77 @@ +#pragma once +#include "CoreLayer/Geometry/Geometry.h" +class SdfFunctions { +public: + enum class Function { + Knob, + KnobInner, + KnobOuter, + TwoSpheres + }; + + static std::string funcEnumToString(Function func) { + switch (func) { + case Function::Knob: + return "knob"; + case Function::KnobInner: + return "knob_inner"; + case Function::KnobOuter: + return "knob_outer"; + case Function::TwoSpheres: + return "two_spheres"; + default: + return ""; + } + } + static Function funcStringToEnum(const std::string &func) { + if (func == "knob") { + return Function::Knob; + } else if (func == "knob_inner") { + return Function::KnobInner; + } else if (func == "knob_outer") { + return Function::KnobOuter; + } else if (func == "two_spheres") { + return Function::TwoSpheres; + } + return Function::Knob; + } + + static double knob(Point3d p); + static double knob_inner(Point3d p); + static double knob_outer(Point3d p); + static double two_spheres(Point3d p); + + static double eval(Function fn, Point3d p) { + switch (fn) { + case Function::Knob: return knob(p); + case Function::KnobInner: return knob_inner(p); + case Function::KnobOuter: return knob_outer(p); + case Function::TwoSpheres: return two_spheres(p); + } + } + + /*template + static Vec3d grad(sdf func, Point3d p) { + constexpr double eps = 0.001; + + std::array vals = { + func(p + Vec3d(eps, 0., 0.)), + func(p + Vec3d(0., eps, 0.)), + func(p + Vec3d(0., 0., eps)), + func(p)}; + + return Vec3d( + vals[0] - vals[3], + vals[1] - vals[3], + vals[2] - vals[3]) / + eps; + } + static Vec3d grad(Function fn, Point3d p) { + switch (fn) { + case Function::Knob: return grad(knob, p); + case Function::KnobInner: return grad(knob_inner, p); + case Function::KnobOuter: return grad(knob_outer, p); + case Function::TwoSpheres: return grad(two_spheres, p); + } + }*/ +}; \ No newline at end of file