diff --git a/include/common.hpp b/include/common.hpp
new file mode 100644
index 0000000..6f70f09
--- /dev/null
+++ b/include/common.hpp
@@ -0,0 +1 @@
+#pragma once
diff --git a/include/line.hpp b/include/line.hpp
index c65dcae..6819ead 100644
--- a/include/line.hpp
+++ b/include/line.hpp
@@ -32,50 +32,37 @@ class ILine {
 public:
     virtual ~ILine() = default;
 
-    virtual Vec3 eval(double param) const = 0;
+    virtual Vec3 eval(real param) = 0;
 
-    virtual Vec3 der1(double param) const = 0;
+    virtual Vec3 der1(real param) = 0;
 
-    virtual Vec3 der2(double param) const = 0;
+    virtual Vec3 der2(real param) = 0;
 
-    virtual ClosestDescOnSeg getClosestParam(const Vec3 &p) const = 0;
+    virtual ClosestDescOnSeg getClosestParam(const Vec3 &p) = 0;
 };
 
-template<size_t N>
-using PtArray = std::vector<Vec<N> >;
-using Pt3Array = PtArray<3>;
-using Pt2Array = PtArray<2>;
-
-
 class Polyline : public ILine {
 public:
-    using Point = Vec<3>;
+    using Point = Vec3;
 
-    Polyline(const Pt3Array &points, const std::vector<double> &bugles, const Vec3 &normal, bool closed = false)
-            : _points(points), _bugles(bugles), _closed(closed), _normal(normal.normalize()) {
-        assert(points.size() >= 2);
+    Polyline(Pt3Array points, std::vector<real> bugles, const Vec3 &normal, bool closed = false)
+            : _points(std::move(points)), _bugles(std::move(bugles)), _closed(closed), _normal(normal.normalize()) {
+        assert(_points.size() >= 2);
         if (closed) {
-            assert(points.size() == bugles.size());
+            assert(_points.size() == _points.size());
         } else {
-            assert(points.size() - 1 == bugles.size());
-        }
-        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()];
-            Vec3 AB = B - A;
-            Vec3 ABNorm = AB.normalize();
-            Vec3 QONorm = normal.cross(ABNorm) * (abs(bugles[i]) > 1 ? -1 : 1);
-            real theta = std::atan(bugles[i]) * 4;
-            circularArcs[i].h = AB.length() * 0.5 * std::tan(theta * 0.5);
-            circularArcs[i].center = A + AB * 0.5 + QONorm * circularArcs[i].h;
-            circularArcs[i].theta = theta;
-            circularArcs[i].radius = (circularArcs[i].center - A).length();
-            circularArcs[i].u = (A - circularArcs[i].center).normalize();
-            circularArcs[i].v = ABNorm.cross(circularArcs[i].u);
+            assert(_points.size() - 1 == _points.size());
         }
     }
 
+    [[nodiscard]] const Pt3Array &getPoints() const { return _points; }
+
+    [[nodiscard]] const std::vector<real> &getBugles() const { return _bugles; }
+
+    [[nodiscard]] const Vec3 &getNormal() const { return _normal; }
+
+    [[nodiscard]] bool isClosed() const { return _closed; }
+
 private:
     Pt3Array _points;
     std::vector<real> _bugles;
@@ -94,8 +81,27 @@ private:
     std::vector<CircularArc> circularArcs;
 
 public:
-    Vec3 eval(real param) const override {
+    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()];
+            Vec3 ABHalf = (B - A) * 0.5;
+            Vec3 ABNorm = ABHalf.normalize();
+            Vec3 QONorm = _normal.cross(ABNorm) * (abs(_bugles[i]) > 1 ? -1 : 1);
+            real theta = std::atan(_bugles[i]) * 4;
+            circularArcs[i].h = ABHalf.norm() * std::tan(theta * 0.5);
+            circularArcs[i].center = A + ABHalf + QONorm * circularArcs[i].h;
+            circularArcs[i].theta = theta;
+            circularArcs[i].radius = (circularArcs[i].center - A).norm();
+            circularArcs[i].u = (A - circularArcs[i].center).normalize();
+            circularArcs[i].v = ABNorm.cross(circularArcs[i].u);
+        }
+    }
+
+    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;
         const auto &arc = circularArcs[seg];
@@ -103,8 +109,9 @@ public:
         return arc.center + arc.radius * (arc.u * std::cos(phi) + arc.v * std::sin(phi));
     }
 
-    Vec3 der1(real param) const override {
+    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;
         const auto &arc = circularArcs[seg];
@@ -112,8 +119,9 @@ public:
         return arc.radius * (arc.u * -std::sin(phi) + arc.v * std::cos(phi));
     }
 
-    Vec3 der2(real param) const override {
+    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;
         const auto &arc = circularArcs[seg];
@@ -121,7 +129,7 @@ public:
         return -arc.radius * (arc.u * std::cos(phi) + arc.v * std::cos(phi));
     }
 
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) const override {
+    ClosestDescOnSeg getClosestParam(const Vec3 &p) override {
         real closestDis = std::numeric_limits<real>::max();
         real closestParam;
         for (int i = 0; i < _bugles.size(); ++i) {
@@ -130,19 +138,19 @@ public:
             const auto &arc = circularArcs[i];
             real dis2Seg = segPtDist(p, A, B).dis;
             if (dis2Seg - arc.h > closestDis) continue;
-            if ((A - p).length() < closestDis) {
-                closestDis = (A - p).length();
+            if ((A - p).norm() < closestDis) {
+                closestDis = (A - p).norm();
                 closestParam = i;
             }
-            if ((B - p).length() < closestDis) {
-                closestDis = (B - p).length();
+            if ((B - p).norm() < closestDis) {
+                closestDis = (B - p).norm();
                 closestParam = i + 1;
             }
             int segInsertedCnt = arc.theta / DISC_ARC_ANGLE;
             for (int j = 0; j < segInsertedCnt; ++j) {
                 real insertParam = i + j * DISC_ARC_ANGLE / arc.theta;
                 const Vec3 insertPt = eval(insertParam);
-                real dis2InsertPt = (p - insertPt).length();
+                real dis2InsertPt = (p - insertPt).norm();
                 if (dis2InsertPt < closestDis) {
                     closestDis = dis2InsertPt;
                     closestParam = insertParam;
@@ -167,17 +175,18 @@ public:
             printf("After iter %d, dL is %lf\n", iter, lDer1);
             if (closestParam < seg - std::numeric_limits<real>::epsilon()) {
                 closestParam = seg;
-                closestDis = (_points[seg] - p).length();
+                closestDis = (_points[seg] - p).norm();
                 break;
             } else if (closestParam > seg + 1 + std::numeric_limits<real>::epsilon()) {
                 closestParam = seg + 1;
-                closestDis = (_points[(seg + 1) % _points.size()] - p).length();
+                closestDis = (_points[(seg + 1) % _points.size()] - p).norm();
                 break;
             }
         }
         return {closestParam, closestDis};
     }
 
+
     void print() const {
         if (_closed) printf("Closed Polyline: \n");
         else printf("Open Polyline: \n");
@@ -202,18 +211,18 @@ private:
         Vec3 AB = B - A;
         Vec3 AP = p - A;
         real h = std::clamp(AP.dot(AB) / AB.dot(AB), 0., 1.);
-        return {h, (AP - AB * h).length()};
+        return {h, (AP - AB * h).norm()};
     }
 
 };
 
 class PolynomialLine : public ILine {
 public:
-    Vec3 eval(double param) const override { return {}; };
+    Vec3 eval(real param) override { return {}; };
 
-    Vec3 der1(double param) const override { return {}; };
+    Vec3 der1(real param) override { return {}; };
 
-    Vec3 der2(double param) const override { return {}; };
+    Vec3 der2(real param) override { return {}; };
 
-    ClosestDescOnSeg getClosestParam(const Vec3 &p) const override { return {}; };
+    ClosestDescOnSeg getClosestParam(const Vec3 &p) override { return {}; };
 };
diff --git a/include/solid.hpp b/include/solid.hpp
index 4b85627..e586614 100644
--- a/include/solid.hpp
+++ b/include/solid.hpp
@@ -1,27 +1,70 @@
 #pragma once
+
 #include <real.hpp>
+#include <utility>
 #include <vec.hpp>
 #include <line.hpp>
 
 class ISolid {
 public:
     virtual ~ISolid() = default;
-
-    virtual real sdf(Vec3 p);
+    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 {
-protected:
-    Pt3Array _profile;
-    real rScale;
+public:
+    Polyline _profile;
+    real _rScale;
+public:
+    IExtrudedSolid(Polyline profile, real rScale) : _profile(std::move(profile)), _rScale(rScale) {
+    }
 };
 
-class extrudedSolidPolyline : public IExtrudedSolid {
-protected:
+class ExtrudedSolidPolyline : public IExtrudedSolid {
+private:
     Polyline _axis;
+    Pt2Array _localProfile2D;
+public:
+    ExtrudedSolidPolyline(Polyline profile, Polyline axis, real rScale) :
+            IExtrudedSolid(std::move(profile), rScale), _axis(std::move(axis)) {
+        // 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 (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);
+        }
+    }
+
+    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
+        {}
+        return 0;
+    }
 };
 
-class extrudedSolidPolynomialLine : public IExtrudedSolid {
+class ExtrudedSolidPolynomialLine : public IExtrudedSolid {
 protected:
     PolynomialLine _axis;
 
diff --git a/include/vec.hpp b/include/vec.hpp
index 4b6917a..a50ae7c 100644
--- a/include/vec.hpp
+++ b/include/vec.hpp
@@ -1,67 +1,75 @@
 #pragma once
+
 #include "real.hpp"
 #include <cmath>
 #include <array>
 #include <assert.h>
+#include <vector>
 
-template<size_t N>
-class Vec {
+template<typename Derived, size_t N>
+class VecBase {
 public:
     std::array<real, N> data;
 
-    Vec() {
+    VecBase() {
         data.fill(0);
     }
 
-    Vec(std::initializer_list<real> values) {
+    VecBase(std::initializer_list<real> values) {
         std::copy(values.begin(), values.end(), data.begin());
     }
 
-    Vec(real... args) : data{static_cast<real>(args)...} {}
+    template<typename... Args>
+    explicit VecBase(Args... args) : data{static_cast<real>(args)...} {
+        static_assert(sizeof...(args) == N, "Argument count must match vector size.");
+    }
 
-    Vec(const Vec<N> &v) {
+    explicit VecBase(const Derived &v) {
         data = v.data;
     }
 
-    Vec<N> &operator=(const Vec<N> &v) {
+    VecBase& operator=(const VecBase &v) {
+        if (this == &v) return *this;
         data = v.data;
         return *this;
     }
 
     real &operator[](size_t index) {
+        if (index >= N) throw std::out_of_range("Index out of range");
         return data[index];
     }
 
     const real &operator[](size_t index) const {
+        if (index >= N) throw std::out_of_range("Index out of range");
         return data[index];
     }
 
-    Vec<N> operator+(const Vec<N> &v) const {
-        Vec<N> result;
+    Derived operator+(const VecBase &v) const {
+        Derived result;
         for (size_t i = 0; i < N; ++i) {
             result[i] = data[i] + v[i];
         }
         return result;
     }
 
-    Vec<N> operator-(const Vec<N> &v) const {
-        Vec<N> result;
+    Derived operator-(const VecBase &v) const {
+        Derived result;
         for (size_t i = 0; i < N; ++i) {
             result[i] = data[i] - v[i];
         }
         return result;
     }
 
-    Vec<N> operator*(real s) const {
-        Vec<N> result;
+    Derived operator*(real s) const {
+        Derived result;
         for (size_t i = 0; i < N; ++i) {
             result[i] = data[i] * s;
         }
         return result;
     }
 
-    friend Vec<N> operator*(real s, const Vec<N> &v) {
-        Vec<N> result;
+    friend Derived operator*(real s, const VecBase &v) {
+        Derived result;
         for (size_t i = 0; i < N; ++i) {
             result[i] = s * v[i];
         }
@@ -69,15 +77,15 @@ public:
     }
 
 
-    Vec<N> operator/(real s) const {
-        Vec<N> result;
+    Derived operator/(real s) const {
+        Derived result;
         for (size_t i = 0; i < N; ++i) {
             result[i] = data[i] / s;
         }
         return result;
     }
 
-    real dot(const Vec<N> &v) const {
+    real dot(const VecBase &v) const {
         real sum = 0;
         for (size_t i = 0; i < N; ++i) {
             sum += data[i] * v[i];
@@ -89,96 +97,67 @@ public:
         return std::sqrt(dot(*this));
     }
 
-    Vec<N> normalize() const {
-        return *this / norm();
+    Derived normalize() const {
+        return static_cast<Derived>(*this / norm());
     }
 
-    Vec<N> reflect(const Vec<N> &n) const {
-        return *this - n * 2 * dot(n);
+    Derived reflect(const Derived &n) const {
+        return static_cast<Derived>(*this - n * 2 * dot(n));
     }
 };
 
 
-// specialize template class Vec<3>;
-template<>
-class Vec<3> {
+// Vec<2> 特化
+class Vec2 : public VecBase<Vec2, 2> {
 public:
     union {
-        std::array<real, 3> data;
-
         struct {
-            real x, y, z;
+            real &x, &y;
         };
-
         struct {
-            real u, v, w;
+            real &u, &v;
+        };
+        struct {
+            real &a, &b;
         };
     };
 
-    Vec() : data{0, 0, 0} {}
-
-    Vec(real x, real y, real z) : data{x, y, z} {}
-
-    Vec(const Vec &v) : data{v.data[0], v.data[1], v.data[2]} {}
-
-    Vec &operator=(const Vec &v) {
-        data[0] = v.data[0];
-        data[1] = v.data[1];
-        data[2] = v.data[2];
-        return *this;
-    }
-
-    real operator[](size_t index) const {
-        return data[index];
-    }
-
-    real &operator[](size_t index) {
-        return data[index];
-    }
-
-    Vec operator+(const Vec &v) const {
-        return {x + v.x, y + v.y, z + v.z};
-    }
-
-    Vec operator-(const Vec &v) const {
-        return {x - v.x, y - v.y, z - v.z};
-    }
+    Vec2(real xVal, real yVal) : VecBase<Vec2, 2>(xVal, yVal), x(data[0]), y(data[1]) {}
 
-    Vec operator*(real s) const {
-        return {x * s, y * s, z * s};
-    }
+    Vec2() : VecBase<Vec2, 2>(), x(data[0]), y(data[1]) {}
 
-    friend Vec operator*(real s, const Vec &v) {
-        return {s * v.x, s * v.y, s * v.z};
-    }
+    Vec2(const Vec2 &v) : VecBase<Vec2, 2>(v), x(data[0]), y(data[1]) {}
+};
 
 
-    Vec operator/(real s) const {
-        assert(s != 0);
-        real inv = 1 / s;
-        return *this * inv;
-    }
+// specialize template class Vec<3>;
+class Vec3 : public VecBase<Vec3, 3> {
+public:
+    union {
+        struct {
+            real &x, &y, &z;
+        };
+        struct {
+            real &u, &v, &w;
+        };
+        struct {
+            real &a, &b, &c;
+        };
+    };
 
-    real dot(const Vec &v) const {
-        return x * v.x + y * v.y + z * v.z;
-    }
+    Vec3() : VecBase<Vec3, 3>(), x(data[0]), y(data[1]), z(data[2]) {}
 
-    real length() const {
-        return std::sqrt(dot(*this));
-    }
+    Vec3(real xVal, real yVal, real zVal) : VecBase<Vec3, 3>(xVal, yVal, zVal), x(data[0]), y(data[1]), z(data[2]) {}
 
-    Vec normalize() const {
-            return *this / length();
-    }
+    Vec3(const Vec3 &vec) : VecBase<Vec3, 3>(vec), x(data[0]), y(data[1]), z(data[2]) {}
 
-    Vec cross(const Vec &v) const {
+    Vec3 cross(const Vec3 &v) const {
         return {y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x};
     }
 };
 
-using Vec2 = Vec<2>;
-using Vec3 = Vec<3>;
-
+using Pt3Array = std::vector<Vec3>;
+using Pt2Array = std::vector<Vec2>;
 
 //class Vec3 {
 //public: