PMClassifier/include/solid.hpp

#pragma once

#include <real.hpp>
#include <utility>
#include <vec.hpp>
#include <line.hpp>

class ISolid {
public:
    virtual ~ISolid() = default;

    virtual real sdf(const Vec3 &p) = 0;
};

Vec2 get2DRepOf3DPt(const Vec3 &pt3D, const Vec3 &u, const Vec3 &v, const Vec3 &localO) {
    Vec3 OP = pt3D - localO;
    return {OP.dot(u), OP.dot(v)};
}

class IExtrudedSolid : public ISolid {
public:
    Polyline _profile;
    real _rScale;

public:
    IExtrudedSolid(Polyline profile, real rScale) : _profile(std::move(profile)), _rScale(rScale) {
    }
};

/**
*  calculate winding number of a point w.r.t. a segment ab
*/
real unsignedWindingNumberSegment(const Vec3 &p, const Vec3 &a, const Vec3 &b, const Vec3 &refNormal) {
    Vec3 pa = a - p;
    Vec3 pb = b - p;
    return std::acos(std::clamp(pa.dot(pb) / (pa.norm() * pb.norm()), static_cast<real>(-1.),
                                static_cast<real>(1.))) / (std::numbers::pi * 2);
}

class ExtrudedSolidPolyline : public IExtrudedSolid {
private:
    Polyline _axis;
    Pt2Array _localProfile2D;
    Pt2Array _localCircleCenter2D;

public:
    ExtrudedSolidPolyline(Polyline profile, Polyline axis, real rScale) : IExtrudedSolid(std::move(profile), rScale),
                                                                          _axis(std::move(axis)) {
        assert(_profile.isClosed());
        // TODO: project profile at st point to 2D
        Vec3 T = _axis.der1(0).normalize();
        Vec3 N = _axis.der2(0).normalize();
        Vec3 B = T.cross(N);
        Vec3 Q = _axis.eval(0);
        for (int i = 0; i < _profile.getPoints().size(); ++i) {
            _localProfile2D.emplace_back(get2DRepOf3DPt(_profile.getPoints()[i] - Q, N, B, Q));
            _localCircleCenter2D.emplace_back(get2DRepOf3DPt(_profile.getCircularArcs()[i].center - Q, N, B, Q));
        }
    }

    real sdf(const Vec3 &p) override {
        ClosestDescOnSeg closestDesc = _axis.getClosestParam(p);
        // TNB coordinate system
        auto t = closestDesc.t;
        Vec3 T = _axis.der1(t).normalize();
        Vec3 N = _axis.der2(t).normalize();
        Vec3 B = T.cross(N);
        Vec3 Q = _axis.eval(t);
        Vec3 QP = p - Q;
        auto uv = get2DRepOf3DPt(QP, N, B, Q);
        //test it
        {
        }
        // TODO: to test if uv is in _localProfile2D
        for (auto i = 0; i < _localProfile2D.size(); ++i) {
        }
        return 0;
    }

private:
    /**
     * in + in = out
     * in + out = in
     * out + in = in
     * out + out = out
     */
    bool isInside2DPolyline(const Polyline& outline, const Vec3& p3D, const Vec2& p2D) {
        assert(outline.isClosed());
        int intersectionCount = 0, segCount = outline.getBugles().size();
        int onSegIdx = -1;
        constexpr int numRays = 3; // 射线数量
        int majorityIn = 0; // 在多边形内的射线计数
        int majorityOut = 0; // 在多边形外的射线计数

        for (int rayIdx = 0; rayIdx < numRays && onSegIdx == -1; ++rayIdx) {
            double angle = (2.0 * std::numbers::pi * rayIdx) / numRays;
            Vec2 rayDir(cos(angle), sin(angle));
            int crossings = 0;

            for (int i = 0; i < segCount; ++i) {
                const Vec2 &a = _localProfile2D[i];
                const Vec2 &b = _localProfile2D[(i + 1) % segCount];
                if (isPointOnSegment(p2D, a, b)) {
                    onSegIdx = i;
                    break;
                }
                // 使用向量方法计算射线和边的交点
                double dx1 = b[0] - a[0];
                double dy1 = b[1] - a[1];
                double dx2 = rayDir[0];
                double dy2 = rayDir[1];
                double determinant = dx1 * dy2 - dy1 * dx2;

                // 如果determinant为0，则射线和边平行，不计算交点
                if (isEqual(determinant, 0)) continue;

                double t1 = ((p2D[0] - a[0]) * dy2 - (p2D[1] - a[1]) * dx2) / determinant;
                double t2 = ((p2D[0] - a[0]) * dy1 - (p2D[1] - a[1]) * dx1) / determinant;

                // 检查交点是否在边上（0 <= t1 <= 1）且射线上（t2 >= 0）
                if (t1 >= 0 && t1 <= 1 && t2 >= 0) {
                    crossings++;
                }
            }
            if (crossings % 2 == 0) {
                majorityOut++;
            } else {
                majorityIn++;
            }
        }
        // 判断是否在扇内
        bool inFan = false;
        for (int i = 0; i < segCount; ++i) {
            const Vec2& a = _localProfile2D[i];
            const Vec2& b = _localProfile2D[(i + 1) % segCount];
            real po = (p2D - _localCircleCenter2D[i]).norm();
            if (po == _profile.getCircularArcs()[i].radius) {
                // TODO
            }
            if ((po < _profile.getCircularArcs()[i].radius) {
                if ((p2D - a).dot(b - a) > 0) {
                    inFan = true;
                    break;
                }
            }
        }
    }

    bool isPointOnSegment(const Vec2 &p, const Vec2 &a, const Vec2 &b) {
        // check collinearity
        double crossProduct = (p[1] - a[1]) * (b[0] - a[0]) - (p[0] - a[0]) * (b[1] - a[1]);
        if (!isEqual(crossProduct, 0)) return false; // Not collinear

        // Check if point is within segment bounds
        return (p[0] >= std::min(a[0], b[0]) && p[0] <= std::max(a[0], b[0]) &&
                p[1] >= std::min(a[1], b[1]) && p[1] <= std::max(a[1], b[1]));
    }

    real wnCircularArc(const Vec3 &p, const Vec3 &a, const Vec3 &b, const Vec3 &plgNormal,
                       const Polyline::CircularArc &arc, int dir) {
        Vec3 pa = a - p;
        Vec3 pb = b - p;
        real wn = std::acos(std::clamp(pa.dot(pb) / (pa.norm() * pb.norm()), static_cast<real>(-1.),
                                       static_cast<real>(1.))) / (std::numbers::pi * 2);
        auto inOutCircle = arc.inCircleCheck(p);
        if (inOutCircle == PtBoundaryRelation::Outside || pa.cross(pb).dot(plgNormal) < 0) {
            // outside
            // pa.cross(pb).dot(plgNormal) 不会 == 0
            return -wn * dir;
        }
        if (inOutCircle == PtBoundaryRelation::Inside) {
            return wn * dir;
        }

        return 0;
    }
};

class ExtrudedSolidPolynomialLine : public IExtrudedSolid {
protected:
    PolynomialLine _axis;
};