winding number of arc

include/common.hpp

@@ -1 +1,7 @@
 #pragma once
+
+enum PtBoundaryRelation {
+    Inside = 1,
+    OnBoundary,
+    Outside = -1
+};

include/line.hpp

@@ -8,6 +8,7 @@
 #include <cmath>
 #include <algorithm>
 #include <numbers>
+#include "common.hpp"
 
 // class ILineParam {
 // public:
@@ -53,6 +54,7 @@ public:
         } else {
             assert(_points.size() - 1 == _points.size());
         }
+        circularArcs.resize(_bugles.size());
     }
 
     [[nodiscard]] const Pt3Array &getPoints() const { return _points; }
@@ -63,12 +65,6 @@ public:
 
     [[nodiscard]] bool isClosed() const { return _closed; }
 
-private:
-    Pt3Array _points;
-    std::vector<real> _bugles;
-    Vec3 _normal;
-    bool _closed;
-
     struct CircularArc {
         Vec3 center;
         real radius;
@@ -76,13 +72,23 @@ private:
         real h;
         Vec3 u;  // dir of OA
         Vec3 v;  // u X v = normal
+
+        PtBoundaryRelation inCircleCheck(const Vec3 &pt) const {
+            real d = (pt - center).norm();
+            return d < radius ? Inside : d > radius ? Outside : OnBoundary;
+        }
     };
 
+private:
+    Pt3Array _points;
+    std::vector<real> _bugles;
+    Vec3 _normal;
+    bool _closed;
+
     std::vector<CircularArc> circularArcs;
 
 public:
     void initSegInfo() {
-        circularArcs.resize(_bugles.size());
         for (size_t i = 0; i < _bugles.size(); ++i) {
             const Point &A = _points[i];
             const Point &B = _points[(i + 1) % _points.size()];
@@ -102,7 +108,6 @@ public:
     }
 
     Vec3 eval(real param) override {
-        assert(param >= 0 && param <= _bugles.size());
         if (circularArcs.empty()) initSegInfo();
         int seg = static_cast<int>(param);
         real tOnSeg = param - seg;
@@ -112,7 +117,6 @@ public:
     }
 
     Vec3 der1(real param) override {
-        assert(param >= 0 && param <= _bugles.size());
         if (circularArcs.empty()) initSegInfo();
         int seg = static_cast<int>(param);
         real tOnSeg = param - seg;
@@ -122,7 +126,6 @@ public:
     }
 
     Vec3 der2(real param) override {
-        assert(param >= 0 && param <= _bugles.size());
         if (circularArcs.empty()) initSegInfo();
         int seg = static_cast<int>(param);
         real tOnSeg = param - seg;
@@ -188,6 +191,7 @@ public:
         return {closestParam, closestDis};
     }
 
+    const std::vector<CircularArc> &getCircularArcs() const { return circularArcs; }
 
     void print() const {
         if (_closed) printf("Closed Polyline: \n");

include/solid.hpp

@@ -8,6 +8,7 @@
 class ISolid {
 public:
     virtual ~ISolid() = default;
+
     virtual real sdf(const Vec3 &p) = 0;
 };
 
@@ -20,31 +21,44 @@ 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 wnSegment(const Vec3 &p, const Vec3 &a, const Vec3 &b, const Vec3 &refNormal) {
+    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.)));
+    return wn * (refNormal.dot(pa.cross(pb)) > 0 ? 1 : -1);
+}
+
 class ExtrudedSolidPolyline : public IExtrudedSolid {
 private:
     Polyline _axis;
     Pt2Array _localProfile2D;
+
 public:
-    ExtrudedSolidPolyline(Polyline profile, Polyline axis, real rScale) :
+    ExtrudedSolidPolyline(Polyline profile, Polyline axis, real rScale) : IExtrudedSolid(std::move(profile), rScale),
-            IExtrudedSolid(std::move(profile), rScale), _axis(std::move(axis)) {
+                                                                          _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);
-        const Pt3Array& profilePts = _profile.getPoints();
+        for (const auto &P: _profile.getPoints()) {
-        for (int i = 0; i < profilePts.size(); ++i) {
-            const Vec3& P = profilePts[i];
             Vec3 QP = P - Q;
             auto uv = get2DRepOf3DPt(QP, N, B, Q);
             // test it
-            {}
+            {
-            _localProfile2D.emplace_back(uv.x(), uv.y());
+            }
+            _localProfile2D.emplace_back(uv.u(), uv.v());
         }
     }
 
@@ -58,14 +72,30 @@ public:
         Vec3 Q = _axis.eval(t);
         Vec3 QP = p - Q;
         auto uv = get2DRepOf3DPt(QP, N, B, Q);
-        // test it
+        //test it
+        {
+        }
         // TODO: to test if uv is in _localProfile2D
+        for (auto i = 0; i < _localProfile2D.size(); ++i) {
+        }
         return 0;
     }
+private:
+    real wnCircularArc(const Vec3& p, const Vec3& a, const Vec3& b, const Vec3& refNormal, const Polyline::CircularArc& arc) {
+        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.)));
+        real dir = refNormal.dot(pb.cross(pa)) > 0 ? 1 : -1;  // 注意这里是pb x pa， 不是pa x pb
+        auto inOutCircle = arc.inCircleCheck(p);
+        if (inOutCircle == OnBoundary) {
+            // TODO
+        }
+        return (inOutCircle - wn) * dir;
+    }
 };
 
 class ExtrudedSolidPolynomialLine : public IExtrudedSolid {
 protected:
     PolynomialLine _axis;
-
 };