From 1aec62db6e96234a18208d8c175c8a31e5ebbb57 Mon Sep 17 00:00:00 2001 From: Cchen77 <1204918763@qq.com> Date: Mon, 7 Oct 2024 14:31:49 +0800 Subject: [PATCH] add global conditioning for gpis --- scenes/gpis-1/scene.json | 4 +- scenes/gpis-2/scene.json | 21 +++- .../GaussianProcess/GaussianProcess.cpp | 107 +++++++++++++----- .../GaussianProcess/GaussianProcess.h | 36 +++--- .../GaussianProcessFactory.cpp | 21 +++- src/FunctionLayer/Medium/GPISMedium.cpp | 2 +- 6 files changed, 142 insertions(+), 49 deletions(-) diff --git a/scenes/gpis-1/scene.json b/scenes/gpis-1/scene.json index a93f59f8..0b3e0eab 100644 --- a/scenes/gpis-1/scene.json +++ b/scenes/gpis-1/scene.json @@ -1,6 +1,6 @@ { "renderer" : { - "spp" : 32, + "spp" : 1, "output_file" : "gpis" }, @@ -8,7 +8,7 @@ { "name" : "gpis", "type" : "gpis", - "num_sample_points": 32, + "num_sample_points": 8, "marching_step_size": 10000, "gaussian_process":{ "type":"std", diff --git a/scenes/gpis-2/scene.json b/scenes/gpis-2/scene.json index 76a02d45..e217bc4c 100644 --- a/scenes/gpis-2/scene.json +++ b/scenes/gpis-2/scene.json @@ -1,6 +1,6 @@ { "renderer" : { - "spp" : 64, + "spp" : 32, "output_file" : "gpis" }, @@ -25,6 +25,25 @@ "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":{ diff --git a/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp b/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp index 0dc7a24e..c0815e40 100644 --- a/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp +++ b/src/FunctionLayer/GaussianProcess/GaussianProcess.cpp @@ -27,7 +27,7 @@ void GPRealization::makeIntersection(size_t p, double offset) { points.push_back(zeroCrossing); derivativeTypes.push_back(DerivativeType::None); - values.push_back(lerp(preV,curV,offset)); + values.push_back(lerp(preV, curV, offset)); // place holder derivativeDirections.push_back({}); @@ -124,7 +124,34 @@ void GPRealization::applyMemoryModel(Vec3d rayDir, MemoryModel memoryModel) { } } -GaussianProcess::GaussianProcess(std::shared_ptr _mean, std::shared_ptr _cov) : meanFunction(_mean), covFunction(_cov) { +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 { @@ -134,7 +161,6 @@ GPRealization GaussianProcess::sample(const Point3d *points, const DerivativeTyp auto v = mvn.sample(sampler); for (int i = 0; i < numPoints; ++i) { values.push_back(v(i)); - //values.push_back((*meanFunction)(DerivativeType::None, points[i], {})); } return GPRealization(this, points, derivativeTypes, derivativeDirs, values.data(), numPoints, derivativeDir); @@ -147,11 +173,11 @@ GPRealization GaussianProcess::sampleCond(const Point3d *points, const Derivativ return sample(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir, sampler); } - Eigen::MatrixXd kCC = covSym(pointsCond, derivativeTypesCond, derivativeDirsCond, derivativeDirCond, numPointsCond); - Eigen::MatrixXd kCx = cov(pointsCond, derivativeTypesCond, derivativeDirsCond, derivativeDirCond, numPointsCond, - points, derivativeTypes, derivativeDirs, derivativeDir, numPoints); + 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 method from https://github.com/daseyb/gpis-light-transport/blob/main/src/core/math/GaussianProcess.cpp + // 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) { @@ -169,7 +195,7 @@ GPRealization GaussianProcess::sampleCond(const Point3d *points, const Derivativ Eigen::BDCSVD solver; solver.compute(kCC.triangularView()); - if (solver.info() != Eigen::ComputationInfo::Success) {; + if (solver.info() != Eigen::ComputationInfo::Success) { std::cerr << "Conditioning decomposition failed (BDCSVD)!\n"; } @@ -179,16 +205,17 @@ GPRealization GaussianProcess::sampleCond(const Point3d *points, const Derivativ } } - Eigen::VectorXd _mean = mean(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir) + + 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 = covSym(points, derivativeTypes, derivativeDirs, derivativeDir, numPoints) - + 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)); - //values.push_back((*meanFunction)(DerivativeType::None, points[i], {})); } return GPRealization(this, points, derivativeTypes, derivativeDirs, values.data(), numPoints, derivativeDir); } @@ -202,8 +229,8 @@ Eigen::VectorXd GaussianProcess::meanPrior(const Point3d *points, const Derivati return _mean; } -Eigen::MatrixXd GaussianProcess::covPrior(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, const Vec3d &derivativeDirX, size_t numPointsX, - const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, const Vec3d &derivativeDirY, size_t numPointsY) const { +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; @@ -215,7 +242,7 @@ Eigen::MatrixXd GaussianProcess::covPrior(const Point3d *pointsX, const Derivati return _cov; } -Eigen::MatrixXd GaussianProcess::covPriorSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, const Vec3d &derivativeDir, size_t numPoints) const { +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; @@ -229,23 +256,41 @@ Eigen::MatrixXd GaussianProcess::covPriorSym(const Point3d *points, const Deriva return _cov; } -// TODO(Cchen77): applying global condition,for now we just use prior version Eigen::VectorXd GaussianProcess::mean(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir) const { - return meanPrior(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); -} + if (globalCondition.size() == 0) { + return meanPrior(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); + } + auto mean_prior = meanPrior(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir); -Eigen::MatrixXd GaussianProcess::cov(const Point3d *pointsX, const DerivativeType *derivativeTypesX, const Vec3d *derivativeDirsX, const Vec3d &derivativeDirX, size_t numPointsX, - const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, const Vec3d &derivativeDirY, size_t numPointsY) const { - return covPrior(pointsX, derivativeTypesX, derivativeDirsX, derivativeDirX, numPointsX, - pointsY, derivativeTypesY, derivativeDirsY, derivativeDirY, numPointsY); + 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)))); } -double GaussianProcess::goodStepSize(Point3d p, Vec3d rd, double desiredCov) const{ - return 0.; +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, const Vec3d &derivativeDir, size_t numPoint) const { - return covPriorSym(points, derivativeTypes, derivativeDirs, derivativeDir, numPoint); +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 { @@ -261,6 +306,16 @@ std::tuple GaussianProcess::meanAndCov(const P _cov(i, j) = _cov(j, i) = cov_i_j; } } - // TODO(Cchen77): consider gloabl conditioning when construct mean and cov + 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 index f349cc45..3befab43 100644 --- a/src/FunctionLayer/GaussianProcess/GaussianProcess.h +++ b/src/FunctionLayer/GaussianProcess/GaussianProcess.h @@ -5,6 +5,8 @@ #include "Eigen/Dense" #include "GPFunctions.h" + +#include enum class MemoryModel { None, GlobalN, @@ -23,7 +25,7 @@ struct GPRealization { std::vector values; // there is a zero-crossing between values[p-1] ~ values[p] - virtual void makeIntersection(size_t p,double offset); + 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); @@ -46,38 +48,36 @@ struct GPRealization { } }; struct GaussianProcess { - GaussianProcess(std::shared_ptr _mean, std::shared_ptr _cov); + 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; - GPRealization sampleCond(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, size_t numPoints, const Vec3d &derivativeDir, - const GPRealization &condRealization, - Sampler &sampler) const { - return sampleCond(points, derivativeTypes, derivativeDirs, numPoints, derivativeDir, - condRealization.points.data(), condRealization.derivativeTypes.data(), condRealization.derivativeDirections.data(), condRealization.values.data(), 0, {}, - sampler); - } - std::shared_ptr meanFunction; std::shared_ptr covFunction; virtual double goodStepSize(Point3d p, Vec3d rd, double desiredCov) const; - // TODO(Cchen77): Global Conditioning +protected: GPRealization globalCondition; + std::variant, Eigen::BDCSVD> globalCondtionSolver; + void initGlobalCondition(const GPRealization &_globalCondition); -protected: // 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, const Vec3d &derivativeDirX, size_t numPointsX, - const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, const Vec3d &derivativeDirY, size_t numPointsY) const; - Eigen::MatrixXd covPriorSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, const Vec3d &derivativeDir, size_t numPoints) 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, const Vec3d &derivativeDirX, size_t numPointsX, - const Point3d *pointsY, const DerivativeType *derivativeTypesY, const Vec3d *derivativeDirsY, const Vec3d &derivativeDirY, size_t numPointsY) const; - Eigen::MatrixXd covSym(const Point3d *points, const DerivativeType *derivativeTypes, const Vec3d *derivativeDirs, const Vec3d &derivativeDir, size_t numPoint) 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 index 3d4ea317..b42b2a58 100644 --- a/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.cpp +++ b/src/FunctionLayer/GaussianProcess/GaussianProcessFactory.cpp @@ -5,7 +5,26 @@ std::shared_ptr GaussianProcessFactory::LoadGaussianProcessFrom if (gpType == "std") { auto meanFunction = MeanFunctionFactory::LoadMeanFunctionFromJson(json["mean"]); auto covarianceFunction = CovarianceFunctionFactory::LoadCovarianceFunctionFromJson(json["covariance"]); - return std::make_shared(meanFunction, covarianceFunction); + 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; diff --git a/src/FunctionLayer/Medium/GPISMedium.cpp b/src/FunctionLayer/Medium/GPISMedium.cpp index e203f106..c18b7bbb 100644 --- a/src/FunctionLayer/Medium/GPISMedium.cpp +++ b/src/FunctionLayer/Medium/GPISMedium.cpp @@ -125,7 +125,7 @@ bool GPISMedium::intersectGP(const Ray &ray, GPRealization &gpRealization, doubl // use last realization as conditon else { gpRealization = gaussianProcess->sampleCond( - points.data(), derivativeTypes.data(), nullptr, marchingNumSamplePoints, {}, gpRealization, sampler); + points.data(), derivativeTypes.data(), nullptr, marchingNumSamplePoints, {}, EXPAND_GPREALIZATION_WITH_VALUE(gpRealization), sampler); } double lastV = gpRealization.values[0];