Browse Source

closest pt; 2d projection; vec class refoctor

master
Dtouch 4 months ago
parent
commit
690cc7ce13
  1. 1
      include/common.hpp
  2. 103
      include/line.hpp
  3. 59
      include/solid.hpp
  4. 143
      include/vec.hpp

1
include/common.hpp

@ -0,0 +1 @@
#pragma once

103
include/line.hpp

@ -32,50 +32,37 @@ class ILine {
public: public:
virtual ~ILine() = default; 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 { class Polyline : public ILine {
public: public:
using Point = Vec<3>; using Point = Vec3;
Polyline(const Pt3Array &points, const std::vector<double> &bugles, const Vec3 &normal, bool closed = false) Polyline(Pt3Array points, std::vector<real> bugles, const Vec3 &normal, bool closed = false)
: _points(points), _bugles(bugles), _closed(closed), _normal(normal.normalize()) { : _points(std::move(points)), _bugles(std::move(bugles)), _closed(closed), _normal(normal.normalize()) {
assert(points.size() >= 2); assert(_points.size() >= 2);
if (closed) { if (closed) {
assert(points.size() == bugles.size()); assert(_points.size() == _points.size());
} else { } else {
assert(points.size() - 1 == bugles.size()); assert(_points.size() - 1 == _points.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);
} }
} }
[[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: private:
Pt3Array _points; Pt3Array _points;
std::vector<real> _bugles; std::vector<real> _bugles;
@ -94,8 +81,27 @@ private:
std::vector<CircularArc> circularArcs; std::vector<CircularArc> circularArcs;
public: 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()); assert(param >= 0 && param <= _bugles.size());
if (circularArcs.empty()) initSegInfo();
int seg = static_cast<int>(param); int seg = static_cast<int>(param);
real tOnSeg = param - seg; real tOnSeg = param - seg;
const auto &arc = circularArcs[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)); 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()); assert(param >= 0 && param <= _bugles.size());
if (circularArcs.empty()) initSegInfo();
int seg = static_cast<int>(param); int seg = static_cast<int>(param);
real tOnSeg = param - seg; real tOnSeg = param - seg;
const auto &arc = circularArcs[seg]; const auto &arc = circularArcs[seg];
@ -112,8 +119,9 @@ public:
return arc.radius * (arc.u * -std::sin(phi) + arc.v * std::cos(phi)); 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()); assert(param >= 0 && param <= _bugles.size());
if (circularArcs.empty()) initSegInfo();
int seg = static_cast<int>(param); int seg = static_cast<int>(param);
real tOnSeg = param - seg; real tOnSeg = param - seg;
const auto &arc = circularArcs[seg]; const auto &arc = circularArcs[seg];
@ -121,7 +129,7 @@ public:
return -arc.radius * (arc.u * std::cos(phi) + arc.v * std::cos(phi)); 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 closestDis = std::numeric_limits<real>::max();
real closestParam; real closestParam;
for (int i = 0; i < _bugles.size(); ++i) { for (int i = 0; i < _bugles.size(); ++i) {
@ -130,19 +138,19 @@ public:
const auto &arc = circularArcs[i]; const auto &arc = circularArcs[i];
real dis2Seg = segPtDist(p, A, B).dis; real dis2Seg = segPtDist(p, A, B).dis;
if (dis2Seg - arc.h > closestDis) continue; if (dis2Seg - arc.h > closestDis) continue;
if ((A - p).length() < closestDis) { if ((A - p).norm() < closestDis) {
closestDis = (A - p).length(); closestDis = (A - p).norm();
closestParam = i; closestParam = i;
} }
if ((B - p).length() < closestDis) { if ((B - p).norm() < closestDis) {
closestDis = (B - p).length(); closestDis = (B - p).norm();
closestParam = i + 1; closestParam = i + 1;
} }
int segInsertedCnt = arc.theta / DISC_ARC_ANGLE; int segInsertedCnt = arc.theta / DISC_ARC_ANGLE;
for (int j = 0; j < segInsertedCnt; ++j) { for (int j = 0; j < segInsertedCnt; ++j) {
real insertParam = i + j * DISC_ARC_ANGLE / arc.theta; real insertParam = i + j * DISC_ARC_ANGLE / arc.theta;
const Vec3 insertPt = eval(insertParam); const Vec3 insertPt = eval(insertParam);
real dis2InsertPt = (p - insertPt).length(); real dis2InsertPt = (p - insertPt).norm();
if (dis2InsertPt < closestDis) { if (dis2InsertPt < closestDis) {
closestDis = dis2InsertPt; closestDis = dis2InsertPt;
closestParam = insertParam; closestParam = insertParam;
@ -167,17 +175,18 @@ public:
printf("After iter %d, dL is %lf\n", iter, lDer1); printf("After iter %d, dL is %lf\n", iter, lDer1);
if (closestParam < seg - std::numeric_limits<real>::epsilon()) { if (closestParam < seg - std::numeric_limits<real>::epsilon()) {
closestParam = seg; closestParam = seg;
closestDis = (_points[seg] - p).length(); closestDis = (_points[seg] - p).norm();
break; break;
} else if (closestParam > seg + 1 + std::numeric_limits<real>::epsilon()) { } else if (closestParam > seg + 1 + std::numeric_limits<real>::epsilon()) {
closestParam = seg + 1; closestParam = seg + 1;
closestDis = (_points[(seg + 1) % _points.size()] - p).length(); closestDis = (_points[(seg + 1) % _points.size()] - p).norm();
break; break;
} }
} }
return {closestParam, closestDis}; return {closestParam, closestDis};
} }
void print() const { void print() const {
if (_closed) printf("Closed Polyline: \n"); if (_closed) printf("Closed Polyline: \n");
else printf("Open Polyline: \n"); else printf("Open Polyline: \n");
@ -202,18 +211,18 @@ private:
Vec3 AB = B - A; Vec3 AB = B - A;
Vec3 AP = p - A; Vec3 AP = p - A;
real h = std::clamp(AP.dot(AB) / AB.dot(AB), 0., 1.); 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 { class PolynomialLine : public ILine {
public: 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 {}; };
}; };

59
include/solid.hpp

@ -1,27 +1,70 @@
#pragma once #pragma once
#include <real.hpp> #include <real.hpp>
#include <utility>
#include <vec.hpp> #include <vec.hpp>
#include <line.hpp> #include <line.hpp>
class ISolid { class ISolid {
public: public:
virtual ~ISolid() = default; virtual ~ISolid() = default;
virtual real sdf(const Vec3 &p) = 0;
virtual real sdf(Vec3 p);
}; };
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 { class IExtrudedSolid : public ISolid {
protected: public:
Pt3Array _profile; Polyline _profile;
real rScale; real _rScale;
public:
IExtrudedSolid(Polyline profile, real rScale) : _profile(std::move(profile)), _rScale(rScale) {
}
}; };
class extrudedSolidPolyline : public IExtrudedSolid { class ExtrudedSolidPolyline : public IExtrudedSolid {
protected: private:
Polyline _axis; 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: protected:
PolynomialLine _axis; PolynomialLine _axis;

143
include/vec.hpp

@ -1,67 +1,75 @@
#pragma once #pragma once
#include "real.hpp" #include "real.hpp"
#include <cmath> #include <cmath>
#include <array> #include <array>
#include <assert.h> #include <assert.h>
#include <vector>
template<size_t N> template<typename Derived, size_t N>
class Vec { class VecBase {
public: public:
std::array<real, N> data; std::array<real, N> data;
Vec() { VecBase() {
data.fill(0); data.fill(0);
} }
Vec(std::initializer_list<real> values) { VecBase(std::initializer_list<real> values) {
std::copy(values.begin(), values.end(), data.begin()); 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; data = v.data;
} }
Vec<N> &operator=(const Vec<N> &v) { VecBase& operator=(const VecBase &v) {
if (this == &v) return *this;
data = v.data; data = v.data;
return *this; return *this;
} }
real &operator[](size_t index) { real &operator[](size_t index) {
if (index >= N) throw std::out_of_range("Index out of range");
return data[index]; return data[index];
} }
const real &operator[](size_t index) const { const real &operator[](size_t index) const {
if (index >= N) throw std::out_of_range("Index out of range");
return data[index]; return data[index];
} }
Vec<N> operator+(const Vec<N> &v) const { Derived operator+(const VecBase &v) const {
Vec<N> result; Derived result;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
result[i] = data[i] + v[i]; result[i] = data[i] + v[i];
} }
return result; return result;
} }
Vec<N> operator-(const Vec<N> &v) const { Derived operator-(const VecBase &v) const {
Vec<N> result; Derived result;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
result[i] = data[i] - v[i]; result[i] = data[i] - v[i];
} }
return result; return result;
} }
Vec<N> operator*(real s) const { Derived operator*(real s) const {
Vec<N> result; Derived result;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
result[i] = data[i] * s; result[i] = data[i] * s;
} }
return result; return result;
} }
friend Vec<N> operator*(real s, const Vec<N> &v) { friend Derived operator*(real s, const VecBase &v) {
Vec<N> result; Derived result;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
result[i] = s * v[i]; result[i] = s * v[i];
} }
@ -69,15 +77,15 @@ public:
} }
Vec<N> operator/(real s) const { Derived operator/(real s) const {
Vec<N> result; Derived result;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
result[i] = data[i] / s; result[i] = data[i] / s;
} }
return result; return result;
} }
real dot(const Vec<N> &v) const { real dot(const VecBase &v) const {
real sum = 0; real sum = 0;
for (size_t i = 0; i < N; ++i) { for (size_t i = 0; i < N; ++i) {
sum += data[i] * v[i]; sum += data[i] * v[i];
@ -89,96 +97,67 @@ public:
return std::sqrt(dot(*this)); return std::sqrt(dot(*this));
} }
Vec<N> normalize() const { Derived normalize() const {
return *this / norm(); return static_cast<Derived>(*this / norm());
} }
Vec<N> reflect(const Vec<N> &n) const { Derived reflect(const Derived &n) const {
return *this - n * 2 * dot(n); return static_cast<Derived>(*this - n * 2 * dot(n));
} }
}; };
// specialize template class Vec<3>; // Vec<2> 特化
template<> class Vec2 : public VecBase<Vec2, 2> {
class Vec<3> {
public: public:
union { union {
std::array<real, 3> data;
struct { struct {
real x, y, z; real &x, &y;
}; };
struct { struct {
real u, v, w; real &u, &v;
};
struct {
real &a, &b;
}; };
}; };
Vec() : data{0, 0, 0} {} Vec2(real xVal, real yVal) : VecBase<Vec2, 2>(xVal, yVal), x(data[0]), y(data[1]) {}
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};
}
Vec operator*(real s) const { Vec2() : VecBase<Vec2, 2>(), x(data[0]), y(data[1]) {}
return {x * s, y * s, z * s};
}
friend Vec operator*(real s, const Vec &v) { Vec2(const Vec2 &v) : VecBase<Vec2, 2>(v), x(data[0]), y(data[1]) {}
return {s * v.x, s * v.y, s * v.z}; };
}
Vec operator/(real s) const { // specialize template class Vec<3>;
assert(s != 0); class Vec3 : public VecBase<Vec3, 3> {
real inv = 1 / s; public:
return *this * inv; union {
} struct {
real &x, &y, &z;
};
struct {
real &u, &v, &w;
};
struct {
real &a, &b, &c;
};
};
real dot(const Vec &v) const { Vec3() : VecBase<Vec3, 3>(), x(data[0]), y(data[1]), z(data[2]) {}
return x * v.x + y * v.y + z * v.z;
}
real length() const { Vec3(real xVal, real yVal, real zVal) : VecBase<Vec3, 3>(xVal, yVal, zVal), x(data[0]), y(data[1]), z(data[2]) {}
return std::sqrt(dot(*this));
}
Vec normalize() const { Vec3(const Vec3 &vec) : VecBase<Vec3, 3>(vec), x(data[0]), y(data[1]), z(data[2]) {}
return *this / length();
}
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}; return {y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x};
} }
}; };
using Vec2 = Vec<2>; using Pt3Array = std::vector<Vec3>;
using Vec3 = Vec<3>; using Pt2Array = std::vector<Vec2>;
//class Vec3 { //class Vec3 {
//public: //public:

Loading…
Cancel
Save