#pragma once
#include <array>
#include <bitset>
#include <cassert>
#include <iostream>
#include <booluarray.hpp>
#include <algorithm>
#include <cstddef>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include "bernstein.hpp"
#include "multiloop.hpp"
#include "organizer/blobtree.hpp"
#include "quadrature_multipoly.hpp"
#include "binomial.hpp"

#include "organizer/minimalPrimitive.hpp"

#include "real.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp"
#include "vector"
#include "xarray.hpp"
#include <chrono>
#include <cmath>
#include <memory>

#include "organizer/primitive.hpp"

namespace algoim::organizer
{


namespace detail
{
void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {}

void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res)
{
    assert(0 <= k && k < 3);
    assert(all(res.ext() == remove_component(phi.ext(), k)));
    // assert(0 < facePos && facePos < 1);
    real* basisAlongK;
    int   P = phi.ext(k);
    algoim_spark_alloc(real, &basisAlongK, P);
    bernstein::evalBernsteinBasis(facePos, P, basisAlongK);

    for (auto i = res.loop(); ~i; ++i) {
        real evalAlongKAtFacePos = 0;
        for (int j = 0; j < P; ++j) { evalAlongKAtFacePos += basisAlongK[j] * phi.m(add_component(i(), k, j)); }
        res.l(i) = evalAlongKAtFacePos;
    }
}

template <int N>
int numOfZero(const uvector<int, N>& v)
{
    int res = 0;
    for (int i = 0; i < N; ++i) {
        if (v(i) == 0) { res++; }
    }
    return res;
}

template <int N>
int binary2Decimal(const uvector<int, N>& v, int zeroRep = 0)
{
    int res  = 0;
    int base = 1;
    for (int i = N - 1; i >= 0; --i) {
        if (v(i) != zeroRep) { res += base; }
        base *= 2;
    }
    return res;
}

template <int N>
int replaceFirst(uvector<int, N>& v, int oldVal, int newVal, int startIdx = 0)
{
    for (int i = startIdx; i < N; ++i) {
        if (v(i) == oldVal) {
            v(i) = newVal;
            return i;
        }
    }
    return -1;
}

template <int N>
int findFirst(const uvector<int, N>& v, int val, int startIdx = 0)
{
    for (int i = startIdx; i < N; ++i) {
        if (v(i) == val) return i;
    }
    return -1;
}

AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector<int, 3> side, int negativeRep = -1)
{
    AABB     res = fatherAABB;
    uvector3 mid = res.center();
    for (int i = 0; i < 3; ++i) {
        if (side(i) == negativeRep) {
            res.max(i) = mid(i);
        } else {
            res.min(i) = mid(i);
        }
    }
    return res;
}
} // namespace detail

struct OcTreeNode {
    std::vector<int> polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中
    // std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
    // std::array<int, CHILD_NUM> children;        // octree
    AABB             aabb;
    organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递

    // Node() { children.fill(-1); }
    OcTreeNode() = default;

    OcTreeNode(const uvector3& min_, const uvector3& max_, const organizer::BlobTree& blobTree_)
        : aabb(min_, max_), blobTree(blobTree_) {};
    OcTreeNode(const AABB& aabb_, const organizer::BlobTree& blobTree_) : aabb(aabb_), blobTree(blobTree_) {};

    // 深拷贝
    OcTreeNode(const OcTreeNode& node)
        : polyIntersectIndices(node.polyIntersectIndices),
          // polyFullyContainedIndices(node.polyFullyContainedIndices),
          aabb(node.aabb),
          blobTree(node.blobTree)
    {
        int a = 1;
        int b = 2;
    }
};

bool keepQuadraturePoint(const std::vector<tensor3>& tensors, const OcTreeNode& ocTreeNode, const uvector3& point)
{
    // 只需要考虑intersect polys，不用考虑fully contained polys
    const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices;
    assert(tensors.size() == polyIntersectIndices.size());
    // if (polyIntersectIndices.size() == 0) {
    //     assert(ocTreeNode.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
    //     return ocTreeNode.blobTree.structure.back().inOut == organizer::NODE_IN;
    // }
    std::vector<char> primitiveInOuts(polyIntersectIndices.size());
    for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) {
        // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt);
        primitiveInOuts[i] = isInsideBernstein(tensors[i], point);
    }
    if (primitiveInOuts.size() == 3 && primitiveInOuts[0] == NODE_IN && primitiveInOuts[1] == NODE_IN
        && primitiveInOuts[2] == NODE_IN) {
        int aaa = 1;
        int bbb = 1;
    }
    // 这里blobTree是拷贝传参
    auto blobTree = ocTreeNode.blobTree;
    int  res      = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts);
    assert(res == organizer::NODE_IN || res == organizer::NODE_OUT);
    return res == organizer::NODE_IN;
}

// std::vector<std::shared_ptr<PrimitiveDesc>> primitives;

// BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};

const int CHILD_NUM = 8;

const std::array<int, 2> sides = {-1, 1};

struct BasicTaskRes {
    real volume;
    real surf;
};

// 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合
// TODO: 参数太多了，考虑换用std::function + lambda
void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid,
                                  const int       polyIntersectIndex,
                                  const tensor3&  poly,

                                  // std::array<OcTreeNode, CHILD_NUM>& subNodes,
                                  std::vector<OcTreeNode>& subNodes,
                                  int                      startIdxToCheck,
                                  uvector<int, 3>&         mark)
{
    int zeroIdx = detail::findFirst(mark, 0, startIdxToCheck);
    if (zeroIdx == -1) {
        tensor3 halfCellPoly(nullptr, poly.ext());
        algoim_spark_alloc(real, halfCellPoly);
        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) {
            subNode.polyIntersectIndices.emplace_back(polyIntersectIndex);
        } else {
            organizer::traverse(subNode.blobTree, polyIntersectIndex, false);
        }
    } else {
        for (auto side : sides) {
            mark(zeroIdx) = side;
            visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, zeroIdx + 1, mark);
        }
        mark(zeroIdx) = 0;
    }
}

// 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历
// 所以树结构是不需要的记录的
// void buildOcTree(const Scene& scene, std::vector<Node>& intermediateNodes, std::vector<Node> leaves, int nowNodeIdx)
void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves)
{
    const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
    if (polyIntersectIndices.size() <= 4) {
        leaves.emplace_back(node);
        return;
    }
    std::vector<OcTreeNode> subNodes(CHILD_NUM);
    // std::array<OcTreeNode, CHILD_NUM> subNodes;
    // intermediateNodes.resize(lastIdx + 8);
    int                     subIdx = 0;
    for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
        subNodes[subIdx].aabb     = detail::getOneEightCellAABB(node.aabb, j(), 0);
        subNodes[subIdx].blobTree = node.blobTree;
    }
    uvector3 nodeMid = node.aabb.center();
    for (int i = 0; i < polyIntersectIndices.size(); ++i) {
        const int       polyIntersectIndex = polyIntersectIndices[i];
        const auto&     poly               = scene.minimalReps[polyIntersectIndex].tensor;
        uvector<int, 3> mark(0, 0, 0);
        for (int faceAxis = 0; faceAxis < 3; ++faceAxis) {
            real    centerPlane = nodeMid(faceAxis);
            tensor2 restrictToCenterPlane(nullptr, remove_component(poly.ext(), faceAxis));
            algoim_spark_alloc(real, restrictToCenterPlane);
            // NOTE: 这里restrict To Face用的是合成的tensor，这可能导致产生很多实际不需要考虑的交线
            detail::restrictToInnerFace(poly, faceAxis, centerPlane, restrictToCenterPlane);
            int signOnHalfPlane = bernstein::uniformSign(restrictToCenterPlane);
            if (signOnHalfPlane == -1) {
                // primitive contain the centerface
                mark(faceAxis) = 2;
            } else if (signOnHalfPlane == 1) {
                // primitive intersects either side or both sides of the centerface
                // deCasteljau to transformation
                AABB halfCell          = node.aabb;
                halfCell.max(faceAxis) = nodeMid(faceAxis);
                tensor3 halfCellPoly(nullptr, poly.ext());
                algoim_spark_alloc(real, halfCellPoly);
                bernstein::deCasteljau(poly, halfCell.min, halfCell.max, halfCellPoly);
                if (bernstein::uniformSign(halfCellPoly) != 1) {
                    // 负空间有
                    mark(faceAxis) = -1;
                }
                halfCell.max(faceAxis) = node.aabb.max(faceAxis);
                halfCell.min(faceAxis) = nodeMid(faceAxis);
                bernstein::deCasteljau(poly, halfCell.min, halfCell.max, halfCellPoly);
                if (bernstein::uniformSign(halfCellPoly) != 1) {
                    // 正空间有
                    mark(faceAxis) += 1; // 当正负空间都有，记0
                }
            }
        }

        if (any(mark == uvector3(2, 2, 2))) {
            if (all(mark == uvector3(2, 2, 2))) {
                // fully containing cases
                for (int i = 0; i < CHILD_NUM; ++i) {
                    tensor3 subPoly(nullptr, poly.ext());
                    bernstein::deCasteljau(poly, subNodes[i].aabb.min, subNodes[i].aabb.max, subPoly);
                    if (bernstein::uniformSign(subPoly) == -1) {
                        // subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex);
                        organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, true);
                    } else {
                        subNodes[i].polyIntersectIndices.emplace_back(polyIntersectIndex);
                    }
                }
                continue;
            }
            for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
                subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
            }
            continue;
        }
        int zeroCount = detail::numOfZero<3>(mark);
        if (zeroCount == 0) {
            // mark has -1 or 1 only
            real nodeMidVal = bernstein::evalBernsteinPoly(poly, nodeMid);
            int  markIdx    = detail::binary2Decimal(mark, -1);
            for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
                if (subIdx == markIdx)
                    subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
                else
                    organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
            }
        } else if (zeroCount == 1) {
            // poly related to 2 subcells
            uvector<int, 3> sideMark = mark;
            int             zeroIdx  = detail::replaceFirst(sideMark, 0, 1);
            int             sideIdx1 = detail::binary2Decimal(sideMark, -1);
            sideMark(zeroIdx)        = -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) {
                if (subIdx != sideIdx1 && subIdx != sideIdx2) {
                    organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
                }
            }
        } else {
            // zeroCount == 2 or 3
            visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, 0, mark);
        }
    }
    // launch subdivision in subcells
    for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
        // subNodes[subIdx].polyFullyContainedIndices.resize(subNodes[subIdx].polyFullyContainedIndices.size()
        //                                                   + node.polyFullyContainedIndices.size());
        // subNodes[subIdx].polyFullyContainedIndices.insert(subNodes[subIdx].polyFullyContainedIndices.end(),
        //                                                   node.polyFullyContainedIndices.begin(),
        //                                                   node.polyFullyContainedIndices.end());
        buildOcTreeV1(scene, subNodes[subIdx], leaves);
    }
}

static int deepest = 0;

void buildOcTreeV0(const Scene&             scene,
                   const OcTreeNode&        node,
                   std::vector<OcTreeNode>& leaves,
                   int                      depth,
                   int&                     cnt,
                   int                      depthWithSamePolyCnt,
                   int                      fatherPolyCnt)
{
    if (deepest < depth) {
        deepest = depth;
        std::cout << "depth: " << deepest << std::endl;
    }
    if (cnt == 2) {
        int aaa = 1;
        int bbb = 1;
    }
    const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
    if (polyIntersectIndices.size() <= 6) {
        leaves.emplace_back(node);
        return;
    }
    int d = deepest;
    if (d == 3 && depth == 1) {
        int aaa = 1;
        int bbb = 1;
    }

    // 限制无限递归：
    if (fatherPolyCnt == polyIntersectIndices.size()) {
        depthWithSamePolyCnt++;
    } else {
        depthWithSamePolyCnt = 1;
    }
    if (depthWithSamePolyCnt == 2) {
        leaves.emplace_back(node);
        return;
    }

    if (depth >= 8) {
        leaves.emplace_back(node);
        return;
    }

    // std::array<OcTreeNode, CHILD_NUM> subNodes;
    std::vector<OcTreeNode> subNodes(CHILD_NUM);
    // intermediateNodes.resize(lastIdx + 8);
    int                     subIdx = 0;
    for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
        subNodes[subIdx].aabb                 = detail::getOneEightCellAABB(node.aabb, j(), 0);
        subNodes[subIdx].blobTree             = node.blobTree;
        subNodes[subIdx].polyIntersectIndices = std::vector<int>();
    }
    if (node.blobTree.structure[4].inOut != NODE_IN_OUT_UNKNOWN
        && std::find(polyIntersectIndices.begin(), polyIntersectIndices.end(), 3) != polyIntersectIndices.end()) {
        int aaa = 1;
        int bbb = 1;
    }
    if (cnt == 514) {
        int aaa = 1;
        int bbb = 1;
    }
    uvector3 nodeMid = node.aabb.center();
    for (int i = 0; i < polyIntersectIndices.size(); ++i) {
        const int   polyIntersectIndex = polyIntersectIndices[i];
        const auto& minimalRep         = scene.minimalReps[polyIntersectIndex];
        subIdx                         = 0;

        for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
            if (subIdx == 0) {
                int aaa = 1;
                int bbb = 1;
            }
            if (!minimalRep.aabb.isIntersect(subNodes[subIdx].aabb.transform01ToGlobal(scene.boundary))) {
                // out of the subcell
                organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
                continue;
            }
            tensor3 subcellPoly(nullptr, minimalRep.tensor.ext());
            algoim_spark_alloc(real, subcellPoly);
            bernstein::deCasteljau(minimalRep.tensor, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly);
            int sign = bernstein::uniformSign(subcellPoly);
            if (sign == 1) {
                organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
            } else if (sign == -1) {
                // // 采样一个点，带入所有原始primitive看是否是冗余部分
                // uvector3 sampleX = subNodes[subIdx].aabb.center() * scene.boundary.size() + scene.boundary.min;
                // if (isInsidePowers(poly.rawPower, sampleX))
                organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, true);
                // else
                //     organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
            } else {
                subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
            }
        }
    }
    if (subNodes[1].blobTree.structure[4].inOut != NODE_IN_OUT_UNKNOWN
        && std::find(subNodes[1].polyIntersectIndices.begin(), subNodes[1].polyIntersectIndices.end(), 3)
               != subNodes[1].polyIntersectIndices.end()) {
        int aaa = 1;
        int bbb = 1;
    }
    for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
        buildOcTreeV0(scene, subNodes[subIdx], leaves, depth + 1, ++cnt, depthWithSamePolyCnt, polyIntersectIndices.size());
    }
}

/** single object quadrature */
void basicTask(const std::shared_ptr<PrimitiveDesc>& p, int q = 20, real xmin = -1, real xmax = 1)
{
    real     volume    = 0;
    auto     integrand = [](const uvector<real, 3>& x) { return 1.0; };
    uvector3 range     = xmax - xmin;

    if (auto pt = std::dynamic_pointer_cast<SphereDesc>(p)) {
        tensor3 tensor(nullptr, 3);
        algoim_spark_alloc(real, tensor);
        std::vector<AABB>   aabbs;
        VisiblePrimitiveRep visiblePrimitiveRep{{tensor}, aabbs, AABB(), BlobTree()};
        makeSphere(*pt, visiblePrimitiveRep);
        detail::powerTransformation(range, xmin, tensor);
        detail::power2BernsteinTensor(tensor);

        ImplicitPolyQuadrature<3> ipquad(tensor);
        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
            if (isInsideBernstein(tensor, x)) volume += w * integrand(xmin + x * range);
        });
    } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
        const int faceCount = pt->faces.size();
        assert(faceCount > 1);
        std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 2));
        algoim_spark_alloc(real, planeTensors);
        std::vector<AABB>   aabbs;
        VisiblePrimitiveRep visiblePrimitiveRep{planeTensors, aabbs, AABB(), BlobTree()};
        AABB                aabb;
        makeMesh(*pt, visiblePrimitiveRep);
        for (auto& tensor : planeTensors) {
            detail::powerTransformation(range, xmin, tensor);
            detail::power2BernsteinTensor(tensor);
        }

        ImplicitPolyQuadrature<3> ipquad(planeTensors);
        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
            if (isInsideBernsteins(planeTensors, x)) volume += w * integrand(xmin + x * range);
        });
    }
    volume *= prod(range);
    std::cout << "Volume xxx: " << volume << std::endl;
};

void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, int q = 20, real xmin = -1, real xmax = 1)
{
    std::vector<SparkStack<real>*> phiStacks;
    std::vector<tensor3>           phis;
    real                           volume;
    auto                           integrand = [](const uvector<real, 3>& x) { return 1.0; };
    uvector3                       range     = xmax - xmin;
    uvector<real, 3>               testX(0., 0.75, 0.2);
    for (int i = 0; i < primitives.size(); i++) {
        if (auto pt = std::dynamic_pointer_cast<SphereDesc>(primitives[i])) {
            tensor3 tensor(nullptr, 3);
            phiStacks.emplace_back(algoim_spark_alloc_heap(real, tensor));
            std::vector<AABB>   aabbs;
            VisiblePrimitiveRep visiblePrimitiveRep{{tensor}, aabbs, AABB(), BlobTree()};
            makeSphere(*pt, visiblePrimitiveRep);

            detail::powerTransformation(range, xmin, tensor);
            detail::power2BernsteinTensor(tensor);
            phis.emplace_back(tensor);
        } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
            const int faceCount = pt->faces.size();
            assert(faceCount > 1);
            std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 3));
            algoimSparkAllocHeapVector(phiStacks, planeTensors);
            std::vector<AABB>   aabbs;
            VisiblePrimitiveRep visiblePrimitiveRep{planeTensors, aabbs, AABB(), BlobTree()};
            makeMesh(*pt, visiblePrimitiveRep);

            for (auto& tensor : planeTensors) {
                detail::powerTransformation(range, xmin, tensor);
                detail::power2BernsteinTensor(tensor);
            }
            phis.insert(phis.end(), planeTensors.begin(), planeTensors.end());
        }
    }
    real                      testEvalBernstein = bernstein::evalBernsteinPoly(phis[0], testX);
    ImplicitPolyQuadrature<3> ipquad(phis);

    ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
        auto realX = x * range + xmin;
        if (isInsideBernsteins(phis, x)) volume += w * integrand(realX);
    });
    volume *= prod(range);
    std::cout << "Volume xxx: " << volume << std::endl;

    // free memory, further deallocating memory of xarray
    for (auto& p : phiStacks) delete p;
    // for (auto& p : originTensorStacks) delete p;
};

BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
{
    if (node.polyIntersectIndices.size() == 0) {
        assert(node.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
        if (node.blobTree.structure.back().inOut == NODE_IN) {
            return {node.aabb.volume(), 0};
        } else {
            return {0, 0};
        }
    }

    auto                           integrand = [](const uvector<real, 3>& x) { return 1.0; };
    real                           volume = 0., surf = 0.;
    auto                           range = node.aabb.size();
    std::vector<tensor3>           phis; // phis using subcell as [0,1]^3
    std::vector<SparkStack<real>*> phiStacks;

    uvector3i numSamples(20, 20, 20);
    // if (false) {
    if (node.polyIntersectIndices.size() >= 5) {
        // TODO: 蒙特卡洛
        std::cout << "Monte Carlo 1" << std::endl;
        ;
        uvector3f interval = range / (numSamples);
        int       inCnt    = 0;
        for (int i = 0; i < node.polyIntersectIndices.size(); i++) {
            const auto& polyWithinScene = scene.minimalReps[node.polyIntersectIndices[i]].tensor;
            phis.emplace_back(tensor3(nullptr, polyWithinScene.ext()));
            phiStacks.emplace_back(algoim_spark_alloc_heap(real, phis.back()));
            phis[i] = polyWithinScene; // copy, not deCasteljau
        }
        for (auto i = MultiLoop<3>(0, numSamples); ~i; ++i) {
            uvector3f x = i() * interval + node.aabb.min + interval * 0.5;
            if (keepQuadraturePoint(phis, node, x)) ++inCnt;
        }

        for (auto& p : phiStacks) delete p;
        return {node.aabb.volume() * inCnt / prod(numSamples), 0};
    }
    for (int i = 0; i < node.polyIntersectIndices.size(); i++) {
        const auto& polyWithinScene = scene.minimalReps[node.polyIntersectIndices[i]].tensor;
        phis.emplace_back(tensor3(nullptr, polyWithinScene.ext()));
        phiStacks.emplace_back(algoim_spark_alloc_heap(real, phis.back()));
        bernstein::deCasteljau(polyWithinScene, node.aabb.min, node.aabb.max, phis.back());
    }
    // ImplicitPolyQuadrature<3> ipquad(scene.minimalReps, node.polyIntersectIndices);
    ImplicitPolyQuadrature<3> ipquad(phis);
#if STOP_WHEN_BLOCKED
    // if (stopCurrentQuadrature
    //     || std::chrono::duration_cast<std::chrono::seconds>(std::chrono::high_resolution_clock::now() -
    //     timerStart).count()
    //            > MAX_DURATION) {
    if (stopCurrentQuadrature
        || std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - timerStart).count() > MAX_DURATION) {
        std::cout << "Monte Carlo 2" << std::endl;
        stopCurrentQuadrature = false;
        // std::cout << "blocked" << std::endl;
        // 注意这里tensor已经经过deCasteljau变换，所以采样点属于[0,1]^3
        uvector3f interval    = 1.f / (numSamples);
        int       inCnt       = 0;
        for (auto i = MultiLoop<3>(0, numSamples); ~i; ++i) {
            uvector3f x = i() * interval + interval * 0.5;
            if (keepQuadraturePoint(phis, node, x)) ++inCnt;
        }

        for (auto& p : phiStacks) delete p;
        return {node.aabb.volume() * inCnt / prod(numSamples), 0};
    }
#endif
    ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
        auto realX = x * range + node.aabb.min; // 这里realX应该是最原始空间下的点？不过因为算体积，所以不影响
        if (keepQuadraturePoint(phis, node, x)) { volume += w * integrand(realX); }
    });
#if STOP_WHEN_BLOCKED
    // if (stopCurrentQuadrature
    //     || std::chrono::duration_cast<std::chrono::seconds>(std::chrono::high_resolution_clock::now() -
    //     timerStart).count()
    //            > MAX_DURATION) {
    if (stopCurrentQuadrature
        || std::chrono::duration<double>(std::chrono::high_resolution_clock::now() - timerStart).count() > MAX_DURATION) {
        std::cout << "Monte Carlo 2" << std::endl;
        stopCurrentQuadrature = false;
        // std::cout << "blocked" << std::endl;
        // 注意这里tensor已经经过deCasteljau变换，所以采样点属于[0,1]^3
        uvector3f interval    = 1.f / (numSamples);
        int       inCnt       = 0;
        for (auto i = MultiLoop<3>(0, numSamples); ~i; ++i) {
            uvector3f x = i() * interval + interval * 0.5;
            if (keepQuadraturePoint(phis, node, x)) ++inCnt;
        }

        for (auto& p : phiStacks) delete p;
        return {node.aabb.volume() * inCnt / prod(numSamples), 0};
    }
#endif

    for (auto& p : phiStacks) delete p;
    double aaa = volume * node.aabb.volume();
    return {volume * node.aabb.volume(), surf};
}

void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives,
                     const uvector3&                                    xmin,
                     const uvector3&                                    xmax,
                     organizer::BlobTree&                               blobTree,
                     const real                                         q = 10)
{
    std::vector<OcTreeNode> leaves;

    std::vector<SparkStack<real>*>   tensorStacks;
    std::vector<VisiblePrimitiveRep> visiblePrimitiveReps(primitives.size());

    real             volume    = 0;
    auto             integrand = [](const uvector<real, 3>& x) { return 1.0; };
    uvector3         range     = xmax - xmin;
    uvector<real, 3> testX(0.8, 0.8, 0.8);
    for (int i = 0; i < primitives.size(); i++) {
        if (auto pt = std::dynamic_pointer_cast<SphereDesc>(primitives[i])) {
            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3)};
            auto& phi                       = visiblePrimitiveReps[i].tensors[0];
            tensorStacks.emplace_back(algoim_spark_alloc_heap(real, phi));

            makeSphere(*pt, visiblePrimitiveReps[i]);

            detail::powerTransformation(range, xmin, phi);
            detail::power2BernsteinTensor(phi);
            visiblePrimitiveReps[i].aabb.normalize(range, xmin);
        } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
            const int faceCount = pt->faces.size();
            assert(faceCount > 1);
            visiblePrimitiveReps[i].tensors = std::vector<tensor3>(faceCount, tensor3(nullptr, 3));
            auto& planeTensors              = visiblePrimitiveReps[i].tensors;
            algoimSparkAllocHeapVector(tensorStacks, planeTensors);

            makeMesh(*pt, visiblePrimitiveReps[i]);

            for (auto& tensor : planeTensors) {
                detail::powerTransformation(range, xmin, tensor);
                detail::power2BernsteinTensor(tensor);
            }
            visiblePrimitiveReps[i].aabb.normalize(range, xmin);
        } else if (auto pt = std::dynamic_pointer_cast<CylinderDesc>(primitives[i])) {
            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)};
            auto& tensors                   = visiblePrimitiveReps[i].tensors;
            tensors[0].ext_(pt->alignAxis)  = 1;
            algoimSparkAllocHeapVector(tensorStacks, tensors);

            makeCylinder(*pt, visiblePrimitiveReps[i]);

            for (auto& tensor : tensors) {
                detail::powerTransformation(range, xmin, tensor);
                detail::power2BernsteinTensor(tensor);
            }
            visiblePrimitiveReps[i].aabb.normalize(range, xmin);
        } else if (auto pt = std::dynamic_pointer_cast<ConeDesc>(primitives[i])) {
            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)};
            auto& tensors                   = visiblePrimitiveReps[i].tensors;
            algoimSparkAllocHeapVector(tensorStacks, tensors);

            makeCone(*pt, visiblePrimitiveReps[i]);

            for (auto& tensor : tensors) {
                detail::powerTransformation(range, xmin, tensor);
                detail::power2BernsteinTensor(tensor);
            }
            visiblePrimitiveReps[i].aabb.normalize(range, xmin);
        } else if (auto pt = std::dynamic_pointer_cast<HalfPlaneDesc>(primitives[i])) {
            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3)};
            auto& tensor                    = visiblePrimitiveReps[i].tensors[0];
            tensorStacks.emplace_back(algoim_spark_alloc_heap(real, tensor));

            makeHalfPlane(*pt, visiblePrimitiveReps[i]);

            detail::powerTransformation(range, xmin, tensor);
            detail::power2BernsteinTensor(tensor);

            visiblePrimitiveReps[i].aabb.normalize(range, xmin);
        } else {
            std::cerr << "unrecognized primitive type." << std::endl;
        }
    }

    std::vector<MinimalPrimitiveRep> minimalReps;
    /*** merge subtrees to main tree ***/
    // std::vector<int>                 realLeafIndices;
    // for (int i = 0; i < visiblePrimitiveReps.size() - 1; ++i) {
    //     // 必须以自左向右的顺序访问所有叶节点
    //     assert(blobTree.primitiveNodeIdx[i] < blobTree.primitiveNodeIdx[i + 1]);
    // }

    // /***注意这个写法的bug***/
    // /***blobTree.structure[j].ancestor一边被修改一边作为条件判断依据***/
    // /***需要对原始的ancestor存一份副本***/
    // // std::vector<int> oldAncestors(blobTree.structure.);
    // // for (int i = 0; i < visiblePrimitiveReps.size(); ++i) {
    // //     int originLeafIdx = blobTree.primitiveNodeIdx[i];
    // //     for (int j = 0; j < originLeafIdx; ++j) {
    // //         if (blobTree.structure[j].isLeft && blobTree.structure[j].ancestor + j > originLeafIdx) {
    // //             blobTree.structure[j].ancestor += std::max(int(visiblePrimitiveReps[i].subBlobTree.structure.size()) -
    // 1, 0);
    // //         }
    // //     }
    // // }
    // /***注意这个写法的bug***/

    // for (int i = 0; i < blobTree.structure.size(); ++i) {
    //     int oldAncestor = blobTree.structure[i].ancestor;
    //     for (int j = visiblePrimitiveReps.size() - 1; blobTree.primitiveNodeIdx[j] > i; --j) {
    //         if (blobTree.structure[i].isLeft && oldAncestor + i > blobTree.primitiveNodeIdx[j]) {
    //             blobTree.structure[i].ancestor += std::max(int(visiblePrimitiveReps[j].subBlobTree.structure.size()) - 1,
    //             0);
    //         }
    //     }
    // }
    // for (int i = 0; i < visiblePrimitiveReps.size(); ++i) {
    //     int originLeafIdx   = blobTree.primitiveNodeIdx[i];
    //     int subBlobTreeSize = visiblePrimitiveReps[i].subBlobTree.structure.size();
    //     if (visiblePrimitiveReps[i].tensors.size() != 1) {
    //         for (int j = i + 1; j < visiblePrimitiveReps.size(); ++j) {
    //             blobTree.primitiveNodeIdx[j] += std::max(int(subBlobTreeSize) - 1, 0);
    //         }
    //         blobTree.structure[originLeafIdx].isPrimitive = false;
    //         blobTree.structure[originLeafIdx].nodeOp      = visiblePrimitiveReps[i].subBlobTree.structure.back().nodeOp;
    //         blobTree.structure.insert(blobTree.structure.begin() + originLeafIdx,
    //                                   visiblePrimitiveReps[i].subBlobTree.structure.begin(),
    //                                   visiblePrimitiveReps[i].subBlobTree.structure.end() - 1);

    //         realLeafIndices.reserve(realLeafIndices.size() +
    //         visiblePrimitiveReps[i].subBlobTree.primitiveNodeIdx.size()); for (auto primitiveIdx :
    //         visiblePrimitiveReps[i].subBlobTree.primitiveNodeIdx) {
    //             realLeafIndices.push_back(primitiveIdx + originLeafIdx);
    //         }
    //         minimalReps.reserve(minimalReps.size() + visiblePrimitiveReps[i].tensors.size());
    //         const auto& aabb = visiblePrimitiveReps[i].aabb;
    //         for (const auto& tensor : visiblePrimitiveReps[i].tensors) {
    //             minimalReps.emplace_back(MinimalPrimitiveRep{tensor, aabb});
    //         }
    //     } else {
    //         blobTree.structure[originLeafIdx].isPrimitive = true;
    //         realLeafIndices.push_back(originLeafIdx);
    //         minimalReps.emplace_back(MinimalPrimitiveRep{visiblePrimitiveReps[i].tensors[0],
    //         visiblePrimitiveReps[i].aabb});
    //     }
    // }
    // blobTree.primitiveNodeIdx = realLeafIndices;
    /*** merge subtrees to main tree ***/

    mergeSubtree2Leaf(blobTree, minimalReps, visiblePrimitiveReps);
    Scene      scene{minimalReps, AABB(xmin, xmax)};
    OcTreeNode rootNode(0, 1, blobTree);
    for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
    int cnt = 1;
    buildOcTreeV0(scene, rootNode, leaves, 1, cnt, 0, -1);
    std::cout << "octree built over" << std::endl;
    // basicTask(scene, leaves[14], q);

    int i = 0;

    for (int i = 0; i < leaves.size(); ++i) {
        const auto& leaf     = leaves[i];
        auto        basicRes = basicTask(scene, leaf, q);
        if (std::isinf(basicRes.volume)) { std::cout << "inf volume when solving leaf: " << i << std::endl; }
        volume += basicRes.volume;
        std::cout << "Solved leaves: " << ++i << "/" << leaves.size() << ", primitive cnt: " << leaf.polyIntersectIndices.size()
                  << ", volume: " << basicRes.volume << std::endl;
    }

    volume *= prod(xmax - xmin);
    std::cout << "Volume xxx: " << volume << std::endl;

    // free memory, further deallocating memory of xarray
    for (auto& p : tensorStacks) delete p;
}
}; // namespace algoim::organizer

// clang-format off
/** 
	最一开始的cell应当包含所有的primitive
	func(given cell) {
		得到8个subcell
		for each primitive {
			mark = [0,0,0]
			for k in 0,1,2 {
				测试primitive和k轴向的centerface有无交
				if 无交 {0
					if primitive包裹centerface {
						mark[k] = 2
					} else {
						判断primitive位于哪一侧
						if k正方向 mark[k] = 1
						else mark[k] = -1
					}
				}
			}
			if mark 含2 {
				// TODO: 可能需要优化：在极端情况下，
				// 例如primitive的边刚好与centerface重合，
				// 此时不是所有subcell都关联primitive
				将primitive放入所有subcell
			} elif mark 没有0 {
				将primitive放入对应的subcell
			} elif mark 有1个0 {
				将primitive放入两个对应的subcell
			} else { // >= 2个0
				测试每个可能有交subcell是否与primitive相交
				相交的放入对于subcell
			}
		}
		for each subcell {
			func(subcell)
		}
	}
	
 */

// clang-format on