simple test case & debug & luanch.json

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.25)
 project(PMClassifier)
 
-set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 
 include_directories(./include)

include/common.hpp

@@ -1,12 +1,11 @@
 #pragma once
 #include "real.hpp"
 #include <limits>
-#include <numbers>
 #include <cmath>
 
 enum PtBoundaryRelation { Inside = -1, OnBoundary, Outside = 1 };
 
-const real PI = std::numbers::pi;
+const real PI = 3.141592653589793238462643383279502884Q;
 const real PI2 = 2 * PI;
 const real EPS = std::numeric_limits<real>::epsilon() * 1e2;
 const real EPS_END_PARAM = std::numeric_limits<real>::epsilon() * 1e3;

include/line.hpp

@@ -10,8 +10,6 @@
 #include <iostream>
 #include <limits.h>
 #include <limits>
-#include <memory>
-#include <numbers>
 #include <vector>
 
 // class ILineParam {
@@ -45,17 +43,17 @@ public:
     ILine(ILine &&) = default;
     ILine &operator=(ILine &&) = default;
 
-    virtual Vec3 eval(real t) = 0;
+    virtual Vec3 eval(real t) const = 0;
 
-    virtual Vec3 der1(real t) = 0;
+    virtual Vec3 der1(real t) const = 0;
 
-    virtual Vec3 der2(real t) = 0;
+    virtual Vec3 der2(real t) const = 0;
 
-    virtual Vec3 tangent(real t) = 0;
+    virtual Vec3 tangent(real t) const = 0;
 
-    virtual Vec3 normal(real t, const Vec3 &tan = -1.) = 0;
+    virtual Vec3 normal(real t, const Vec3 &tan = -1.) const = 0;
 
-    virtual ClosestDescOnSeg getClosestParam(const Vec3 &p) = 0;
+    virtual ClosestDescOnSeg getClosestParam(const Vec3 &p) const = 0;
 
     [[nodiscard]] virtual bool isEndParam(real t) const = 0;
 };
@@ -96,7 +94,7 @@ struct AA {
     void print() const { std::cout << a << " " << b << std::endl; }
 };
 
-const real DISC_ARC_ANGLE = std::numbers::pi * 0.125;
+const real DISC_ARC_ANGLE = PI * 0.125;
 
 class Polyline : public ILine {
 public:
@@ -105,11 +103,12 @@ public:
           _refNormal(refNormal.normalize()) {
         assert(_points.size() >= 2);
         if (closed) {
-            assert(_points.size() == _points.size());
+            assert(_points.size() == _bugles.size());
         } else {
-            assert(_points.size() - 1 == _points.size());
+            assert(_points.size() - 1 == _bugles.size());
         }
         circularArcs.resize(_bugles.size());
+        initSegInfo();
     }
 
     [[nodiscard]] const Pt3Array &getPoints() const { return _points; }
@@ -139,12 +138,12 @@ public:
         }
     }
 
-    Vec3 eval(real t) override {
-        if (circularArcs.empty())
-            initSegInfo();
+    Vec3 eval(real t) const override {
+        // if (circularArcs.empty())
+        //     initSegInfo();
         int seg = static_cast<int>(t);
         if (isEqual(_bugles[seg], 0)) {
-            return _points[seg] + (_points[(seg + 1) % _bugles.size()] - _points[seg]) * (t - seg);
+            return _points[seg] + (_points[(seg + 1) % _points.size()] - _points[seg]) * (t - seg);
         }
         real tOnSeg = t - seg;
         const auto &arc = circularArcs[seg];
@@ -152,12 +151,10 @@ public:
         return arc.center + arc.radius * (arc.u * std::cos(phi) + arc.v * std::sin(phi));
     }
 
-    Vec3 der1(real t) override {
-        if (circularArcs.empty())
-            initSegInfo();
+    Vec3 der1(real t) const override {
         int seg = static_cast<int>(t);
         if (isEqual(_bugles[seg], 0)) {
-            return _points[(seg + 1) % _bugles.size()] - _points[seg];
+            return _points[(seg + 1) % _points.size()] - _points[seg];
         }
         real tOnSeg = t - seg;
         const auto &arc = circularArcs[seg];
@@ -165,21 +162,22 @@ public:
         return arc.radius * (arc.u * -std::sin(phi) + arc.v * std::cos(phi));
     }
 
-    Vec3 der2(real t) override {
-        if (circularArcs.empty())
-            initSegInfo();
+    Vec3 der2(real t) const override {
         int seg = static_cast<int>(t);
-        assert(!isEqual(_bugles[seg], 0));
+        if (isEqual(_bugles[seg], 0)) {
+            int aaa = 1;
+        }
+        // assert(!isEqual(_bugles[seg], 0));
         real tOnSeg = t - seg;
         const auto &arc = circularArcs[seg];
         real phi = tOnSeg * arc.theta;
         return -arc.radius * (arc.u * std::cos(phi) + arc.v * std::cos(phi));
     }
 
-    Vec3 tangent(real t) override { return der1(t).normalize(); }
+    Vec3 tangent(real t) const override { return der1(t).normalize(); }
     // TODO: 试试https://www.jianshu.com/p/9e4877e3965e算出来的结果
-    Vec3 normal(real t, [[maybe_unused]] const Vec3 &tan = -1.) override {
-        if (!isEqual(_bugles[static_cast<int>(t)], 0)) {
+    Vec3 normal(real t, [[maybe_unused]] const Vec3 &tan = -1.) const override {
+        if (isEqual(_bugles[static_cast<int>(t)], 0)) {
             return -circularArcs[static_cast<int>(t)].inCircleDir;
         }
         // 只有对于圆弧这样的特殊曲线是这样
@@ -247,9 +245,7 @@ public:
     //     return {closestParam, closestDis};
     // }
 
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) override {
-        if (circularArcs.empty())
-            initSegInfo();
+    ClosestDescOnSeg getClosestParam(const Vec3 &p) const override {
         ClosestDescOnSeg closestDes{};
         for (int i = 0; i < _bugles.size(); ++i) {
             const Vec3 &a = _points[i];
@@ -280,7 +276,7 @@ public:
                     real cosTheta = (oa).dot(oClsPt) / R2;
                     real theta = std::acos(cosTheta); // [0, pi]
                     if ((oa.cross(oClsPt)).dot(_refNormal) < 0) {
-                        theta = 2 * std::numbers::pi - theta;
+                        theta = PI2 - theta;
                     }
                     closestDes.t = i + theta / arc.theta;
                 }
@@ -370,24 +366,24 @@ public:
         // _k = _4pi2r / (advancePerRound * advancePerRound + _4pi2r * _r);
         _arcDeltaMaxFactor = _4pi2r / (advancePerRound * advancePerRound + _4pi2r * _r) * ONE_EIGHT;
     }
-    Vec3 eval(real t) override {
+    Vec3 eval(real t) const override {
         real theta = _frequency * t;
         return _axisStart + _axisDir * t + (_u * std::cos(theta) + _v * std::sin(theta)) * _r;
     };
 
-    Vec3 der1(real param) override {
+    Vec3 der1(real param) const override {
         real theta = _frequency * param;
         return _axisDir + _2pir_p * (_v * std::cos(theta) - _u * std::sin(theta));
     };
 
-    Vec3 der2(real param) override {
+    Vec3 der2(real param) const override {
         real theta = _frequency * param;
         return -_4pi2r_p2 * (_u * std::cos(theta) + _v * std::sin(theta));
     };
 
-    Vec3 tangent(real t) override { return der1(t).normalize(); }
+    Vec3 tangent(real t) const override { return der1(t).normalize(); }
 
-    Vec3 normal(real t, const Vec3 &tan = -1.) override {
+    Vec3 normal(real t, const Vec3 &tan = -1.) const override {
         Vec3 der2Vec = this->der2(t);
         if (tan == -1.) {
             Vec3 realTan = tangent(t);
@@ -396,7 +392,7 @@ public:
         return (der2Vec - der2Vec.dot(tan) * tan).normalize();
     }
 
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) override {
+    ClosestDescOnSeg getClosestParam(const Vec3 &p) const override {
         // discretization and traversal
         real startT = 0;
         real endT = SEG_T;
@@ -490,7 +486,6 @@ public:
     ArcLine(const Vec3 &a, const Vec3 &b, real bugle, const Vec3 &refNormal)
         // : _a(a), _b(b), _bugle(bugle), _refNormal(refNormal.normalize()) {}
         : Polyline(Pt3Array{{a, b}}, std::vector<real>{bugle}, refNormal, false) {}
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) override { return {}; };
     [[nodiscard]] const std::vector<real> &getBugles() const = delete;
     [[nodiscard]] const real &getBugle() const { return _bugles[0]; }
 
@@ -502,15 +497,15 @@ public:
 
 class PolynomialLine : public ILine {
 public:
-    Vec3 eval(real t) override { return {}; };
+    Vec3 eval(real t) const override { return {}; };
 
-    Vec3 der1(real t) override { return {}; };
+    Vec3 der1(real t) const override { return {}; };
 
-    Vec3 der2(real t) override { return {}; };
+    Vec3 der2(real t) const override { return {}; };
 
-    Vec3 tangent(real t) override { return {}; }
+    Vec3 tangent(real t) const override { return {}; }
 
-    Vec3 normal(real t, const Vec3 &tan = -1.) override { return {}; }
+    Vec3 normal(real t, const Vec3 &tan = -1.) const override { return {}; }
 
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) override { return {}; };
+    ClosestDescOnSeg getClosestParam(const Vec3 &p) const override { return {}; };
 };

include/solid.hpp

@@ -19,7 +19,7 @@ public:
     ISolid(ISolid &&) = default;
     ISolid &operator=(ISolid &&) = default;
 
-    virtual real sdf(const Vec3 &p) = 0;
+    virtual real sdf(const Vec3 &p) const = 0;
 };
 inline Vec2 get2DRepOf3DPt(const Vec3 &pt3D, const Vec3 &u, const Vec3 &v, const Vec3 &localO) {
     Vec3 OP = pt3D - localO;
@@ -74,10 +74,10 @@ public:
             for (int j = 0; j < segCount; ++j) {
                 // TODO:
                 _localProfiles2D[i][j] =
-                    get2DRepOf3DPt(profile.getPoints()[j] - q, normal, _biNormalStartPt, q);
+                    get2DRepOf3DPt(profile.getPoints()[j], normal, _biNormalStartPt, q);
                 auto &arc2d = _localArcs2d[i][j];
                 const auto &arc3d = profile.getCircularArcs()[j];
-                arc2d.center = get2DRepOf3DPt(arc3d.center - q, normal, _biNormalStartPt, q);
+                arc2d.center = get2DRepOf3DPt(arc3d.center, normal, _biNormalStartPt, q);
                 arc2d.inCircleDir = get2DRepOf3DDir(arc3d.inCircleDir, normal, _biNormalStartPt);
                 arc2d.radius = arc3d.radius;
                 arc2d.theta = arc3d.theta;
@@ -86,7 +86,7 @@ public:
         }
     }
 
-    real sdf(const Vec3 &p) override {
+    real sdf(const Vec3 &p) const override {
         ClosestDescOnSeg closestDescToAxis = _axis.getClosestParam(p);
         // TNB coordinate system
         auto t = closestDescToAxis.t;
@@ -118,9 +118,9 @@ public:
     }
 
 private:
-    virtual std::array<Vec3, 3> getTBN(const Vec3 &p, const Vec3 &q, real t) = 0;
+    virtual std::array<Vec3, 3> getTBN(const Vec3 &p, const Vec3 &q, real t) const = 0;
 
-    inline ClosestDescOnSeg disProfile2D(const Vec2 &p2D) {
+    inline ClosestDescOnSeg disProfile2D(const Vec2 &p2D) const {
         ClosestDescOnSeg closestDescToProfile{};
         for (int i = 0; i < _localArcs2d.size(); ++i) {
             ClosestDescOnSeg closestDescTemp = disLoop2D(p2D, _localProfiles2D[i], _localArcs2d[i]);
@@ -131,7 +131,7 @@ private:
         return closestDescToProfile;
     }
 
-    inline PtBoundaryRelation pmcProfile2d(const Vec2 &p2D) {
+    inline PtBoundaryRelation pmcProfile2d(const Vec2 &p2D) const {
         // PMC
         for (int i = 0; i < _localArcs2d.size(); ++i) {
             PtBoundaryRelation relationTmp = pmcLoop2d(p2D, _localProfiles2D[i], _localArcs2d[i]);
@@ -205,7 +205,7 @@ private:
                 for (int i = 0; i < segCount; ++i) {
                     const Vec2 &a = loop2D[i];
                     const Vec2 &b = loop2D[(i + 1) % segCount];
-                    assert(isPointOnSegment(p, a, b));
+                    assert(!isPointOnSegment(p, a, b));
                     // if (isPointOnSegment(p2D, a, b))
                     // {
                     //     onSegIdx = i;
@@ -327,7 +327,7 @@ public:
         assert(_biNormalStartPt.isParallel(_axis.getRefNormal()));
     }
 
-    std::array<Vec3, 3> getTBN(const Vec3 &p, const Vec3 &q, real t) override {
+    std::array<Vec3, 3> getTBN(const Vec3 &p, const Vec3 &q, real t) const override {
         if (!_axis.isEndParam(t) && std::abs(t - std::round(t)) < EPS) {
             // 衔接点处（注意排除断点处）
             // p到圆弧平面的投影
@@ -387,7 +387,7 @@ public:
         : IExtrudedSolidBase(std::move(profiles), std::move(axis), rScale) {}
 
     std::array<Vec3, 3> getTBN([[maybe_unused]] const Vec3 &a, [[maybe_unused]] const Vec3 &b,
-                               real t) override {
+                               real t) const override {
         Vec3 tangent = _axis.tangent(t);
         Vec3 normal = _axis.normal(t);
         return {tangent, normal, tangent.cross(normal)};

include/vec.hpp

@@ -5,9 +5,10 @@
 #include <cmath>
 #include <array>
 #include <assert.h>
+#include <cstddef>
+#include <ostream>
 #include <vector>
 #include <stdexcept>
-#include <concepts>
 
 template <typename Derived, size_t N>
 class VecBase {
@@ -110,6 +111,16 @@ public:
         return true;
     }
 
+    friend std::ostream& operator<<(std::ostream& os, const VecBase& vec) {
+        os << "Vec" << N << "(";
+        for (size_t i = 0; i < N; ++i) {
+            os<<vec[i];
+            if (i<N-1) os<<", ";
+        }
+        os<<")";
+        return os;
+    }
+
     Derived operator-() const { return static_cast<Derived>(*this * -1); }
 };
 
@@ -181,19 +192,6 @@ public:
 using Pt3Array = std::vector<Vec3>;
 using Pt2Array = std::vector<Vec2>;
 
-template <typename T>
-concept IsVec2 = requires(T t) {
-    t.x();
-    t.y();
-};
-
-template <typename T>
-concept IsVec3 = requires(T t) {
-    t.x();
-    t.y();
-    t.z();
-};
-
 //class Vec3 {
 //public:
 //    real x, y, z;

main.cpp

@@ -1,18 +1,20 @@
 #include <iostream>
+#include "line.hpp"
 #include "solid.hpp"
 
-// TIP To <b>Run</b> code, press <shortcut actionId="Run"/> or
-// click the <icon src="AllIcons.Actions.Execute"/> icon in the gutter.
-int main() {
-    Polyline axis{{}, {}, Vec3()};
-    ExtrudedSolidPolyline a({axis}, axis);
-    const Vec3 queryP = Vec3();
-    auto f = a.sdf(queryP);
+void testPrint(const ISolid &solid, const Vec3 &queryPt) {
+    std::cout << "sdf on " << queryPt << " = " << solid.sdf(queryPt) << std::endl;
+}
 
-    return 0;
+void caseCyliner() {
+    Polyline profile{{{-1, 0., 0}, {1., 0., 0.}}, {1., 1.}, {0, 0, -1}, true};
+    Polyline axis{{{0, 0, 0}, {0, 0, 2}}, {0}, {0, 1, 0}};
+    ExtrudedSolidPolyline extrusion({profile}, axis);
+    testPrint(extrusion, {1, 2, 1});
+    testPrint(extrusion, {1, 2, -1});
 }
 
-// TIP See CLion help at <a
-// href="https://www.jetbrains.com/help/clion/">jetbrains.com/help/clion/</a>.
-//  Also, you can try interactive lessons for CLion by selecting
-//  'Help | Learn IDE Features' from the main menu.
+int main() {
+    caseCyliner();
+    return 0;
+}