From 60f1bdb1f8d2bd570fa97390e2c638e607ce6193 Mon Sep 17 00:00:00 2001 From: gjj Date: Sat, 28 Sep 2024 20:37:28 +0800 Subject: [PATCH] rotation --- algoim/factorial.hpp | 29 +++++++++++++++++++ algoim/organizer/primitive.hpp | 51 +++++++++++++++++++++++++++++++++- algoim/quaternion.hpp | 45 ++++++++++++++++++++++++++++++ gjj/primitiveDebug.hpp | 31 ++++++++++++++++++++- 4 files changed, 154 insertions(+), 2 deletions(-) create mode 100644 algoim/factorial.hpp create mode 100644 algoim/quaternion.hpp diff --git a/algoim/factorial.hpp b/algoim/factorial.hpp new file mode 100644 index 0000000..b8a574d --- /dev/null +++ b/algoim/factorial.hpp @@ -0,0 +1,29 @@ +#pragma once +#include "real.hpp" +#include + +namespace algoim +{ +class Factorial +{ +public: + static std::vector cache; + + Factorial() + { + if (cache.empty()) cache.push_back(1); + } + + static real calculate(int n) + { + if (n == 0) return 1; + if (cache.size() <= n) { + for (int i = cache.size(); i <= n; ++i) cache.push_back(cache.back() * i); + } + return cache[n]; + } +}; + +std::vector Factorial::cache; // 。 + +} // namespace algoim \ No newline at end of file diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 14088dd..5620744 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -14,6 +14,8 @@ #include "binomial.hpp" #include "bernstein.hpp" #include "blobtree.hpp" +#include "quaternion.hpp" +#include "factorial.hpp" #include namespace algoim::organizer @@ -85,6 +87,7 @@ void power2BernsteinTensor(xarray& tensor) } } +// 先缩放后平移，非原地修改 void powerTransformation(const uvector& scale, const uvector& bias, const xarray& origin, @@ -117,6 +120,7 @@ void powerTransformation(const uvector& scale, } } +// 先缩放后平移，原地修改 void powerTransformation(const uvector& scale, const uvector& bias, xarray& phiPower) { std::vector>> dimOrderExpansion; @@ -145,6 +149,51 @@ void powerTransformation(const uvector& scale, const uvector& } } +// Note: 旋转张量似乎会提高结果张量的阶次，需要旋转后再降阶 +void powerRotation(const uvector3& vec, const real theta, const tensor3& originalPower, tensor3& res) +{ + int maxDegree = sum(originalPower.ext()) - 3 + 1; + assert(all(res.ext() == maxDegree)); + // 获得旋转矩阵 + tensor2 rotation(nullptr, uvector2(3, 3)); + algoim_spark_alloc(real, rotation); + Quaternion q(vec, theta); + q.getRotation(rotation); + // 旋转张量 + xarrayInit(res); + for (auto i = originalPower.loop(); ~i; ++i) { + real base = originalPower.l(i); + if (base == 0) continue; + auto exps = i(); + for (MultiLoop<3> jx(0, exps(0) + 1); ~jx; ++jx) { + if (sum(jx()) != exps(0)) continue; + for (MultiLoop<3> jy(0, exps(1) + 1); ~jy; ++jy) { + if (sum(jy()) != exps(1)) continue; + for (MultiLoop<3> jz(0, exps(2) + 1); ~jz; ++jz) { + if (sum(jz()) != exps(2)) continue; + uvector3 sumExps = jx() + jy() + jz(); + real item = base; + for (int dim = 0; dim < 3; ++dim) { + item *= pow(rotation.m(uvector2(0, dim)), jx(dim)); + item *= pow(rotation.m(uvector2(1, dim)), jy(dim)); + item *= pow(rotation.m(uvector2(2, dim)), jz(dim)); + } + item *= Factorial::calculate(exp(0)) + / (Factorial::calculate(jx(0)) * Factorial::calculate(jx(1)) * Factorial::calculate(jx(2))); + item *= Factorial::calculate(exp(1)) + / (Factorial::calculate(jy(0)) * Factorial::calculate(jy(1)) * Factorial::calculate(jy(2))); + item *= Factorial::calculate(exp(2)) + / (Factorial::calculate(jz(0)) * Factorial::calculate(jz(1)) * Factorial::calculate(jz(2))); + assert(sumExps(0) < res.ext(0) && sumExps(1) < res.ext(1) && sumExps(2) < res.ext(2)); + res.m(sumExps) += item; + } + } + } + } +} + +void affineTransformation() {} + uvector3i getCompositeExt(const std::vector& tensors) { uvector3i resExt = 1; @@ -536,7 +585,7 @@ public: uvector3 center; uvector3 amplitude; - SphereDesc(real r_, const uvector3& c_, const uvector3& a_) : PrimitiveDesc(), radius(r_), center(c_), amplitude(a_) {} + SphereDesc(real r_, const uvector3& c_, const uvector3& a_ = 1) : PrimitiveDesc(), radius(r_), center(c_), amplitude(a_) {} void print() override { std::cout << "Sphere Description" << std::endl; } }; diff --git a/algoim/quaternion.hpp b/algoim/quaternion.hpp new file mode 100644 index 0000000..cb3edac --- /dev/null +++ b/algoim/quaternion.hpp @@ -0,0 +1,45 @@ +#pragma once +#include "real.hpp" +#include "uvector.hpp" +#include "xarray.hpp" + +namespace algoim +{ + +class Quaternion +{ +public: + real aw, ai, aj, ak; + + Quaternion(real aw_, real ai_, real aj_, real ak_) : aw(aw_), ai(ai_), aj(aj_), ak(ak_) {} + + Quaternion(real ai_, real aj_, real ak_) : aw(0), ai(ai_), aj(aj_), ak(ak_) {} + + Quaternion(const uvector3& axis, real angle) + { + uvector3 axisNorm = axis / norm(axis); + real halfAngle = angle / 2; + aw = cos(halfAngle); + real sinHalfAngle = sin(halfAngle); + ai = axisNorm(0) * sinHalfAngle; + aj = axisNorm(1) * sinHalfAngle; + ak = axisNorm(2) * sinHalfAngle; + } + + void getRotation(tensor2& rotation) const + { + assert(all(rotation.ext() == uvector2(3))); + rotation.m(uvector2(0, 0)) = aw * aw + ai * ai - aj * aj - ak * ak; + rotation.m(uvector2(0, 1)) = 2 * (ai * aj - aw * ak); + rotation.m(uvector2(0, 2)) = 2 * (ai * ak + aw * aj); + + rotation.m(uvector2(1, 0)) = 2 * (ai * aj + aw * ak); + rotation.m(uvector2(1, 1)) = aw * aw - ai * ai + aj * aj - ak * ak; + rotation.m(uvector2(1, 2)) = 2 * (aj * ak - aw * ai); + + rotation.m(uvector2(2, 0)) = 2 * (ai * ak - aw * aj); + rotation.m(uvector2(2, 1)) = 2 * (aj * ak + aw * ai); + rotation.m(uvector2(2, 2)) = aw * aw - ai * ai - aj * aj + ak * ak; + } +}; +} // namespace algoim diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index d13d522..6a607b3 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -16,6 +16,7 @@ #include "real.hpp" #include "sparkstack.hpp" +#include "utility.hpp" #include "uvector.hpp" #include "vector" #include "xarray.hpp" @@ -406,6 +407,33 @@ void testTensorInverse() inverseTensor(tensor); evalX = evalPower(tensor, testX); std::cout << "after inverse, evalX = " << evalX << std::endl; + // free + for (auto &ptr : sparkStackPtrs) delete ptr; +} + +void testRotation() +{ + auto desc = std::make_shared(SphereDesc(0.1, uvector3(0.1, 0.2, 0.1))); + tensor3 power(nullptr, 3); + algoim_spark_alloc(real, power); + VisiblePrimitiveRep visiblePrimitiveRep{{power}, AABB{}, BlobTree()}; + makeSphere(*desc, visiblePrimitiveRep); + uvector3 testX(0.5, 0.7, 0.2); + real evalX = evalPower(power, testX); + std::cout << "before rotation, evalX = " << evalX << std::endl; + + // tensor3 rotatedPower(nullptr, sum(power.ext()) - 3 + 1); + // algoim_spark_alloc(real, rotatedPower); + // uvector3 axis = uvector3(4, 2, 3); + // axis = axis / norm(axis); + + // organizer::detail::powerRotation(axis, util::pi / 2, power, rotatedPower); + // evalX = evalPower(rotatedPower, testX); + // std::cout << "after rotation, evalX = " << evalX << std::endl; + // bool b = bernstein::autoReduction(rotatedPower, 1e4 * std::numeric_limits::epsilon()); + // std::cout << "auto reduction = " << b << std::endl; + // evalX = evalPower(rotatedPower, testX); + // std::cout << "after reduction, evalX = " << evalX << std::endl; } void testPrimitive() @@ -415,7 +443,8 @@ void testPrimitive() // casePolyhedronSphere(); // testSubDivideWithDeCasteljau(); // testBlob(); - caseScene(); + // caseScene(); + testRotation(); // caseScene1(); // caseScene2(); // caseCone();