#ifndef MEDUSA_BITS_UTILS_GENQUAPOINTS_HPP_
#define MEDUSA_BITS_UTILS_GENQUAPOINTS_HPP_
#include <medusa/bits/utils/stepTransformer.hpp>

namespace mm
{
    /**
     * Input the path of step file, return the quadrature points calculating from Nurbs-DIVG algorithm
     *
     * @param path: step model file path
     * @param h: spacing function, currently only support constant distance
     * @param maxPointsSurf: the maxPointSize to fill the surface
     * @param maxPointsSFill: the maxPointSize to fill the whole cad model
     * @param quaPoints: the Quadrature Points to be returned
     */
    void getQuaPoints(const char *path, double h, int maxPointsSurf, int maxPointsFill, Eigen::MatrixXd &quaPoints)
    {
        Range<NURBSPatch<Vec3d, Vec2d>> patches;
        transferStepToNurbsPatches(path, patches);
        NURBSShape<Vec3d, Vec2d> shape(patches);
        shape.maxPoints(maxPointsSurf);
        shape.seed(1337).boundaryProximity(1e-11);

        // Fill
        DomainDiscretization<Vec3d> oversampled_domain = shape.discretizeBoundaryWithStep(h / 2);

        for (int i = 0; i < oversampled_domain.size(); i++)
        {
            oversampled_domain.normal(i) *= -1;
        }

        // Construct the contains function.
        KDTree<Vec3d> contains_tree;
        oversampled_domain.makeDiscreteContainsStructure(contains_tree); // contains_tree包含oversampled_domain的所有域边界点
        auto contains_function = [&](const Vec3d p)
        {
            return oversampled_domain.discreteContains(p, contains_tree);
        };

        // Fill the boundary normally.
        DomainDiscretization<Vec3d> domain = shape.discretizeBoundaryWithStep(h);

        // Fill the interior with the contains function.
        KDTreeMutable<Vec3d> tree;
        GeneralFill<Vec3d> gf;
        gf.maxPoints(maxPointsFill);
        gf(domain, h, tree, contains_function);
        quaPoints.resize(domain.positions().size(), 3);
        for (int i = 0; i < domain.positions().size(); i++)
        {
            quaPoints(i, 0) = domain.positions()[i].x();
            quaPoints(i, 1) = domain.positions()[i].y();
            quaPoints(i, 2) = domain.positions()[i].z();
        }
    }
};
#endif