diff --git a/algoim/factorial.hpp b/algoim/factorial.hpp index b8a574d..e69a3eb 100644 --- a/algoim/factorial.hpp +++ b/algoim/factorial.hpp @@ -9,21 +9,16 @@ 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); + for (int i = cache.size(); i <= n; ++i) cache.push_back(cache[i - 1] * i); } return cache[n]; } }; -std::vector Factorial::cache; // 。 +std::vector Factorial::cache = {1}; // 。 } // namespace algoim \ No newline at end of file diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 81e576c..2e41173 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -61,7 +61,7 @@ int numOfZero(const uvector& v) } template -int binaryToDecimal(const uvector& v, int zeroRep = 0) +int binary2Decimal(const uvector& v, int zeroRep = 0) { int res = 0; int base = 1; @@ -160,7 +160,7 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, if (zeroIdx == -1) { tensor3 halfCellPoly(nullptr, poly.ext()); algoim_spark_alloc(real, halfCellPoly); - int subIdx = detail::binaryToDecimal(mark, -1); + int subIdx = detail::binary2Decimal(mark, -1); auto& subNode = subNodes[subIdx]; bernstein::deCasteljau(poly, subNode.aabb.min, subNode.aabb.max, halfCellPoly); // 求1/8空间下的表达 if (bernstein::uniformSign(halfCellPoly) != 1) { @@ -256,7 +256,7 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector sideMark = mark; int zeroIdx = detail::replaceFirst(sideMark, 0, 1); - int sideIdx1 = detail::binaryToDecimal(sideMark, -1); + int sideIdx1 = detail::binary2Decimal(sideMark, -1); sideMark(zeroIdx) = -1; - int sideIdx2 = detail::binaryToDecimal(sideMark, -1); + int sideIdx2 = detail::binary2Decimal(sideMark, -1); subNodes[sideIdx1].polyIntersectIndices.emplace_back(polyIntersectIndex); subNodes[sideIdx2].polyIntersectIndices.emplace_back(polyIntersectIndex); for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { @@ -306,7 +306,7 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector& polyIntersectIndices = node.polyIntersectIndices; - if (polyIntersectIndices.size() <= 3) { + if (polyIntersectIndices.size() <= 4) { leaves.emplace_back(node); return; } diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 5620744..0022731 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -152,8 +152,8 @@ 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)); + // int maxDegree = sum(originalPower.ext()) - 3 + 1; + // assert(all(res.ext() == maxDegree)); // 获得旋转矩阵 tensor2 rotation(nullptr, uvector2(3, 3)); algoim_spark_alloc(real, rotation); @@ -178,11 +178,11 @@ void powerRotation(const uvector3& vec, const real theta, const tensor3& origina item *= pow(rotation.m(uvector2(1, dim)), jy(dim)); item *= pow(rotation.m(uvector2(2, dim)), jz(dim)); } - item *= Factorial::calculate(exp(0)) + item *= Factorial::calculate(exps(0)) / (Factorial::calculate(jx(0)) * Factorial::calculate(jx(1)) * Factorial::calculate(jx(2))); - item *= Factorial::calculate(exp(1)) + item *= Factorial::calculate(exps(1)) / (Factorial::calculate(jy(0)) * Factorial::calculate(jy(1)) * Factorial::calculate(jy(2))); - item *= Factorial::calculate(exp(2)) + item *= Factorial::calculate(exps(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; diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 6a607b3..6624ee9 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -10,6 +11,7 @@ #include #include #include "bernstein.hpp" +#include "factorial.hpp" #include "multiloop.hpp" #include "quadrature_multipoly.hpp" #include "binomial.hpp" @@ -203,6 +205,109 @@ void caseScene() quadratureScene(primitiveDescriptions, uvector3(-0.8, -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree); } +void caseScene30Objs() +{ + const int PRIMITIVE_CNT = 30; + std::vector> primitiveDescriptions(PRIMITIVE_CNT); + uvector3 sceneSize = 0.8, sceneCenter = 0.; + for (MultiLoop<3> i(0, 2); ~i; ++i) { + uvector3 center = sceneSize * i() - 0.5 * sceneSize + sceneCenter; + int idx = organizer::detail::binary2Decimal(i()); + primitiveDescriptions[2 * idx] = std::make_shared(CuboidDesc(center, 0.16)); + primitiveDescriptions[2 * idx + 1] = std::make_shared(SphereDesc(0.15, center, 1.)); + } // 16 primitives + for (int dim = 0; dim < 3; ++dim) { + for (MultiLoop<2> i(0, 2); ~i; ++i) { + uvector3 bottomCenter = add_component(i(), dim, 0) * sceneSize - 0.5 * sceneSize + sceneCenter; + primitiveDescriptions[organizer::detail::binary2Decimal(i()) + 4 * dim + 16] = + std::make_shared(CylinderDesc(bottomCenter, 0.1, sceneSize(dim), dim)); + } + } // 12 primitives + uvector3 coneBottom = uvector3(1, 0, 0) * sceneSize - 0.5 * sceneSize + sceneCenter; + primitiveDescriptions[28] = std::make_shared(coneBottom, 0.2, 0.5, 1); + uvector3 pyramidBottom = uvector3(0, 0, 1) * sceneSize - 0.5 * sceneSize + sceneCenter; + std::vector pyramidBottomVertices = { + 0.2 * pyramidBottom + uvector3{-1, 0, -1}, + 0.2 * pyramidBottom + uvector3{1, 0, -1}, + 0.2 * pyramidBottom + uvector3{1, 0, 1 }, + 0.2 * pyramidBottom + uvector3{-1, 0, 1 } + }; + primitiveDescriptions[29] = std::make_shared(pyramidBottomVertices, pyramidBottom + uvector3{0, 0.7, 0}); + + organizer::BlobTree blobTree; + blobTree.structure.resize(PRIMITIVE_CNT * 2 - 1); + blobTree.primitiveNodeIdx.resize(PRIMITIVE_CNT); + blobTree.structure[0] = {1, 0, 0, 0, 1, 2}; // cube1 + blobTree.structure[1] = {1, 0, 0, 0, 0, 0}; // sphere1 + blobTree.structure[2] = {0, OP_DIFFERENCE, 0, 0, 1, 6 - 2}; // Dif of = opNode1 + blobTree.structure[3] = {1, 0, 0, 0, 1, 5 - 3}; // cube2 + blobTree.structure[4] = {1, 0, 0, 0, 0, 0}; // sphere2 + blobTree.structure[5] = {0, OP_DIFFERENCE, 0, 0, 0, 0}; // Dif of = opNode2 + blobTree.structure[6] = {0, OP_UNION, 0, 0, 1, 14 - 6}; // Union of = opNode3 + blobTree.structure[7] = {1, 0, 0, 0, 1, 9 - 7}; // cube3 + blobTree.structure[8] = {1, 0, 0, 0, 0, 0}; // sphere3 + blobTree.structure[9] = {0, OP_DIFFERENCE, 0, 0, 1, 13 - 9}; // Dif of = opNode4 + blobTree.structure[10] = {1, 0, 0, 0, 1, 12 - 10}; // cube4 + blobTree.structure[11] = {1, 0, 0, 0, 0, 0}; // sphere4 + blobTree.structure[12] = {0, OP_DIFFERENCE, 0, 0, 0, 0}; // Dif of = opNode5 + blobTree.structure[13] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode6 + blobTree.structure[14] = {0, OP_UNION, 0, 0, 1, 30 - 14}; // Union of = opNode7 + blobTree.structure[15] = {1, 0, 0, 0, 1, 17 - 15}; // cube5 + blobTree.structure[16] = {1, 0, 0, 0, 0, 0}; // sphere5 + blobTree.structure[17] = {0, OP_DIFFERENCE, 0, 0, 1, 21 - 17}; // Dif of = opNode8 + blobTree.structure[18] = {1, 0, 0, 0, 1, 20 - 18}; // cube6 + blobTree.structure[19] = {1, 0, 0, 0, 0, 0}; // sphere6 + blobTree.structure[20] = {0, OP_DIFFERENCE, 0, 0, 0, 0}; // Dif of = opNode9 + blobTree.structure[21] = {0, OP_UNION, 0, 0, 1, 29 - 21}; // Union of = opNode10 + blobTree.structure[22] = {1, 0, 0, 0, 1, 24 - 22}; // cube7 + blobTree.structure[23] = {1, 0, 0, 0, 0, 0}; // sphere7 + blobTree.structure[24] = {0, OP_DIFFERENCE, 0, 0, 1, 28 - 24}; // Dif of = opNode11 + blobTree.structure[25] = {1, 0, 0, 0, 1, 27 - 25}; // cube8 + blobTree.structure[26] = {1, 0, 0, 0, 0, 0}; // sphere8 + blobTree.structure[27] = {0, OP_DIFFERENCE, 0, 0, 0, 0}; // Dif of = opNode12 + blobTree.structure[28] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode13 + blobTree.structure[29] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode14 + blobTree.structure[30] = {0, OP_UNION, 0, 0, 1, 32 - 30}; // Dif of = opNode15 + + blobTree.structure[31] = {1, 0, 0, 0, 1, 33 - 31}; // cylinder1 + blobTree.structure[32] = {1, 0, 0, 0, 0, 0}; // cylinder2 + blobTree.structure[33] = {0, OP_UNION, 0, 0, 1, 37 - 33}; // Union of = opNode16 + blobTree.structure[34] = {1, 0, 0, 0, 1, 36 - 34}; // cylinder3 + blobTree.structure[35] = {1, 0, 0, 0, 0, 0}; // cylinder4 + blobTree.structure[36] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode17 + blobTree.structure[37] = {0, OP_UNION, 0, 0, 1, 45 - 37}; // Union of = opNode18 + blobTree.structure[38] = {1, 0, 0, 0, 1, 40 - 38}; // cylinder5 + blobTree.structure[39] = {1, 0, 0, 0, 0, 0}; // cylinder6 + blobTree.structure[40] = {0, OP_UNION, 0, 0, 1, 44 - 40}; // Union of = opNode19 + blobTree.structure[41] = {1, 0, 0, 0, 1, 43 - 41}; // cylinder7 + blobTree.structure[42] = {1, 0, 0, 0, 0, 0}; // cylinder8 + blobTree.structure[43] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode20 + blobTree.structure[44] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode21 + blobTree.structure[45] = {0, OP_UNION, 0, 0, 1, 53 - 45}; // Union of = opNode22 + blobTree.structure[46] = {1, 0, 0, 0, 1, 48 - 46}; // cylinder9 + blobTree.structure[47] = {1, 0, 0, 0, 0, 0}; // cylinder10 + blobTree.structure[48] = {0, OP_UNION, 0, 0, 1, 52 - 48}; // Union of = opNode23 + blobTree.structure[49] = {1, 0, 0, 0, 1, 51 - 49}; // cylinder11 + blobTree.structure[50] = {1, 0, 0, 0, 0, 0}; // cylinder12 + blobTree.structure[51] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode24 + blobTree.structure[52] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode25 + blobTree.structure[53] = {0, OP_UNION, 0, 0, 0, 0}; // Union of = opNode26 + + blobTree.structure[54] = {0, OP_UNION, 0, 0, 1, 58 - 54}; // UNION of = opNode27 + + blobTree.structure[55] = {1, 0, 0, 0, 1, 57 - 55}; // cone + blobTree.structure[56] = {1, 0, 0, 0, 0, 0}; // pyramid + blobTree.structure[57] = {0, OP_UNION, 0, 0, 1, 59 - 57}; // UNION of = opNode28 + blobTree.structure[58] = {0, OP_UNION, 0, 0, 0, 0}; // UNION of = opNode29 + + int leafCnt = 0; + for (int i = 0; i < blobTree.structure.size(); ++i) { + if (blobTree.structure[i].isPrimitive) { blobTree.primitiveNodeIdx[leafCnt++] = i; } + } + assert(leafCnt == PRIMITIVE_CNT); + quadratureScene(primitiveDescriptions, sceneCenter - 0.5 * sceneSize - 0.2, sceneCenter + 0.5 * sceneSize + 0.2, blobTree); +} + void caseScene1() { const int PRIMITIVE_CNT = 2; @@ -413,7 +518,7 @@ void testTensorInverse() void testRotation() { - auto desc = std::make_shared(SphereDesc(0.1, uvector3(0.1, 0.2, 0.1))); + auto desc = std::make_shared(SphereDesc(0.1, uvector3(0, 0., 0.))); tensor3 power(nullptr, 3); algoim_spark_alloc(real, power); VisiblePrimitiveRep visiblePrimitiveRep{{power}, AABB{}, BlobTree()}; @@ -422,18 +527,18 @@ void testRotation() 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; + tensor3 rotatedPower(nullptr, sum(power.ext())); + 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, 1e5 * 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() @@ -443,8 +548,10 @@ void testPrimitive() // casePolyhedronSphere(); // testSubDivideWithDeCasteljau(); // testBlob(); + // testRotation(); // caseScene(); - testRotation(); + + caseScene30Objs(); // caseScene1(); // caseScene2(); // caseCone();