From d420ac4d143e7d7a6f927a06f11617c5cc89aafe Mon Sep 17 00:00:00 2001 From: Zhicheng Wang <1627343141@qq.com> Date: Tue, 16 Jun 2026 18:41:25 +0800 Subject: [PATCH] refactor: signed pattern SDF and extrude cap improvements - Change pattern eval_sdf to return double (signed distance) instead of pair - Add reversed parameter to get_local_cap_mat for correct cap orientation - Add has_end_caps flag to extrude, handle closed axis caps properly - Add is_closed field to polyline_axis_descriptor - Use desc.is_closed instead of computing from point/bulge counts - Fix total_theta calculation in helixline (multiply by two_pi) - Fix scale vector initialization (zero-fill w component) --- frontend/src/io_axis.cpp | 5 +- frontend/src/io_primitive.cpp | 4 +- .../interface/primitive/axis/extrude_axis.hpp | 2 +- .../primitive/axis/helixline_extrude_axis.hpp | 2 +- .../primitive/axis/polyline_extrude_axis.hpp | 2 +- .../interface/primitive/extrude/extrude.hpp | 16 +- .../interface/primitive/pattern/pattern.hpp | 2 +- .../primitive/pattern/polyline_pattern.hpp | 2 +- .../interface/primitive/simple/box.hpp | 9 + .../interface/primitive_descriptor.h | 1 + primitive_process/src/base/primitive.cpp | 3 +- .../primitive/axis/helixline_extrude_axis.cpp | 117 +++++++---- .../primitive/axis/polyline_extrude_axis.cpp | 72 ++++--- .../primitive/pattern/polyline_pattern.cpp | 193 +++++++++++++----- .../src/subface/extrude/extrude_side_face.cpp | 10 +- 15 files changed, 290 insertions(+), 150 deletions(-) diff --git a/frontend/src/io_axis.cpp b/frontend/src/io_axis.cpp index b5293a6..f7f57fe 100644 --- a/frontend/src/io_axis.cpp +++ b/frontend/src/io_axis.cpp @@ -4,7 +4,6 @@ auto construct_polyline_extrude_axis(const polyline_axis_descriptor_t& desc) { - bool is_closed = desc.point_number == desc.bulge_number + 1; bool is_ccw{}; const Eigen::Vector3d vert0 = convert_vec3d(desc.points[0]); @@ -58,7 +57,7 @@ auto construct_polyline_extrude_axis(const polyline_axis_descriptor_t& desc) axis_to_world.matrix().col(2) = vec_z; axis_to_world.matrix().col(3) = vert0; - return std::make_tuple(obj, is_closed, axis_to_world); + return std::make_tuple(obj, desc.is_closed, axis_to_world); } auto construct_helixline_extrude_axis(const helixline_axis_descriptor_t& desc) @@ -67,7 +66,7 @@ auto construct_helixline_extrude_axis(const helixline_axis_descriptor_t& desc) Eigen::Vector3d vec_z = convert_vec3d(desc.axis_end) - origin; auto height = vec_z.norm(); vec_z /= height; - auto total_theta = height / desc.advance_per_round; + auto total_theta = height / desc.advance_per_round * two_pi; internal::helixline_extrude_axis obj{}; obj.scale_coeff_x = desc.radius; diff --git a/frontend/src/io_primitive.cpp b/frontend/src/io_primitive.cpp index 2b00527..cda438c 100644 --- a/frontend/src/io_primitive.cpp +++ b/frontend/src/io_primitive.cpp @@ -28,7 +28,7 @@ API primitive* create_param_primitive(const extend_axis_object_t* axis, return internal::new_param_primitive(axis->object, axis->is_closed, convert_mat3x4d(axis->axis_to_world), - pattern, + *pattern, type, static_cast(data_center)); } @@ -44,7 +44,7 @@ API void destroy_primitive(primitive* primitive_ptr) API void primitive_apply_scale(primitive* primitive_ptr, vector3d scale) { if (primitive_ptr == nullptr) return; - Eigen::Vector4d scale_eigen; + Eigen::Vector4d scale_eigen{0, 0, 0, 1}; std::move(&scale.x, &scale.x + 3, scale_eigen.data()); primitive_ptr->apply_transform(internal::transform_type::scale, scale_eigen); } diff --git a/primitive_process/interface/primitive/axis/extrude_axis.hpp b/primitive_process/interface/primitive/axis/extrude_axis.hpp index 3561037..b03482f 100644 --- a/primitive_process/interface/primitive/axis/extrude_axis.hpp +++ b/primitive_process/interface/primitive/axis/extrude_axis.hpp @@ -12,7 +12,7 @@ struct PE_API extrude_axis { virtual double get_closest_param(Eigen::Vector3d p) const = 0; // 输入的p需要在axis局部坐标系下，且输出的局部坐标需要由轴坐标系到子面坐标系 virtual std::pair to_pattern_local(double t, Eigen::Vector3d p) const = 0; - virtual Eigen::AffineCompact3d get_local_cap_mat(double t) const = 0; + virtual Eigen::AffineCompact3d get_local_cap_mat(double t, bool reversed) const = 0; virtual aabb_t get_local_aabb() const = 0; }; diff --git a/primitive_process/interface/primitive/axis/helixline_extrude_axis.hpp b/primitive_process/interface/primitive/axis/helixline_extrude_axis.hpp index 47700a5..3206da1 100644 --- a/primitive_process/interface/primitive/axis/helixline_extrude_axis.hpp +++ b/primitive_process/interface/primitive/axis/helixline_extrude_axis.hpp @@ -10,7 +10,7 @@ struct PE_API helixline_extrude_axis final : public extrude_axis { double get_closest_param(Eigen::Vector3d p) const override; std::pair to_pattern_local(double t, Eigen::Vector3d p) const override; - Eigen::AffineCompact3d get_local_cap_mat(double t) const override; + Eigen::AffineCompact3d get_local_cap_mat(double t, bool reversed) const override; aabb_t get_local_aabb() const override { return k_aabb_up_unit; } diff --git a/primitive_process/interface/primitive/axis/polyline_extrude_axis.hpp b/primitive_process/interface/primitive/axis/polyline_extrude_axis.hpp index f296463..54b1089 100644 --- a/primitive_process/interface/primitive/axis/polyline_extrude_axis.hpp +++ b/primitive_process/interface/primitive/axis/polyline_extrude_axis.hpp @@ -11,7 +11,7 @@ struct PE_API polyline_extrude_axis final : public extrude_axis { double get_closest_param(Eigen::Vector3d p) const override; std::pair to_pattern_local(double t, Eigen::Vector3d p) const override; - Eigen::AffineCompact3d get_local_cap_mat(double t) const override; + Eigen::AffineCompact3d get_local_cap_mat(double t, bool reversed) const override; aabb_t get_local_aabb() const override { diff --git a/primitive_process/interface/primitive/extrude/extrude.hpp b/primitive_process/interface/primitive/extrude/extrude.hpp index bb55712..64212f9 100644 --- a/primitive_process/interface/primitive/extrude/extrude.hpp +++ b/primitive_process/interface/primitive/extrude/extrude.hpp @@ -25,21 +25,24 @@ struct extrude_t final : public primitive { aabb.min().array() -= pattern_max_length; aabb.max().array() += pattern_max_length; - if (closed_axis) + if (closed_axis) { + has_end_caps = false; initialize( { make_pointer_wrapper(axis_transform) }, {{(void*)axis, (void*)pattern}}, aabb); - else { + } else { paired_model_matrix bottom_cap_transform{}; - bottom_cap_transform.local_to_world = axis_to_world * axis->get_local_cap_mat(0); + bottom_cap_transform.local_to_world = axis_to_world * axis->get_local_cap_mat(0, true); bottom_cap_transform.world_to_local = bottom_cap_transform.local_to_world.inverse(); paired_model_matrix top_cap_transform{}; - top_cap_transform.local_to_world = axis_to_world * axis->get_local_cap_mat(1); + top_cap_transform.local_to_world = axis_to_world * axis->get_local_cap_mat(1, false); top_cap_transform.world_to_local = top_cap_transform.local_to_world.inverse(); + has_end_caps = true; + initialize( { make_pointer_wrapper(axis_transform), @@ -55,17 +58,18 @@ struct extrude_t final : public primitive { span> get_subfaces() const override { - return {const_cast*>(subfaces.data()), subfaces[1] == nullptr ? 1u : 3u}; + return {const_cast*>(subfaces.data()), has_end_caps ? 3u : 1u}; } std::vector get_subface_types() const override { - if (subfaces[1] == nullptr) + if (!has_end_caps) return {surface_type::extrude_side}; else return {surface_type::extrude_side, surface_type::plane, surface_type::plane}; } std::array, 3> subfaces{nullptr, nullptr, nullptr}; + bool has_end_caps{}; }; } // namespace internal \ No newline at end of file diff --git a/primitive_process/interface/primitive/pattern/pattern.hpp b/primitive_process/interface/primitive/pattern/pattern.hpp index 336fc0c..5d6e31a 100644 --- a/primitive_process/interface/primitive/pattern/pattern.hpp +++ b/primitive_process/interface/primitive/pattern/pattern.hpp @@ -8,7 +8,7 @@ namespace internal struct PE_API pattern { virtual ~pattern() = default; - virtual std::pair eval_sdf(Eigen::Vector2d p) const = 0; + virtual double eval_sdf(Eigen::Vector2d p) const = 0; // virtual double eval_sdf_grad(Eigen::Vector2d p) const = 0; virtual double get_max_bounding_length() const = 0; diff --git a/primitive_process/interface/primitive/pattern/polyline_pattern.hpp b/primitive_process/interface/primitive/pattern/polyline_pattern.hpp index bed570c..cf844f6 100644 --- a/primitive_process/interface/primitive/pattern/polyline_pattern.hpp +++ b/primitive_process/interface/primitive/pattern/polyline_pattern.hpp @@ -11,7 +11,7 @@ namespace internal struct PE_API polyline_pattern final : public pattern { ~polyline_pattern() override = default; - std::pair eval_sdf(Eigen::Vector2d p) const override; + double eval_sdf(Eigen::Vector2d p) const override; double get_max_bounding_length() const override { diff --git a/primitive_process/interface/primitive/simple/box.hpp b/primitive_process/interface/primitive/simple/box.hpp index d4c8540..82d61ba 100644 --- a/primitive_process/interface/primitive/simple/box.hpp +++ b/primitive_process/interface/primitive/simple/box.hpp @@ -16,6 +16,15 @@ struct box_t final : public primitive { trans_by_target_axis_matrices_ptr}, {{}, {}, {}, {}, {}, {}}, k_aabb_unit); + + for(auto s : get_subfaces()) + { + std::cout << (s->world_to_local * Eigen::Vector4d{0, 0, 0, 1}).transpose() << std::endl; + // std::cout << s->local_to_world.matrix() << std::endl; + // std::cout << (s->local_to_world * Eigen::Vector4d{0, 0, 0, 1}).transpose() << std::endl; + // std::cout << s.get_mark() << std::endl; + s.set_mark(~s.get_mark()); + } } primitive_type get_type() const override { return PRIMITIVE_TYPE_BOX; }; diff --git a/primitive_process/interface/primitive_descriptor.h b/primitive_process/interface/primitive_descriptor.h index a180b11..c61a09a 100644 --- a/primitive_process/interface/primitive_descriptor.h +++ b/primitive_process/interface/primitive_descriptor.h @@ -30,6 +30,7 @@ typedef struct { uint32_t bulge_number; // The bulge number of the polyline double* bulges; // The bulge of each edge vector3d reference_normal; // The reference normal of the polyline + bool is_closed; } polyline_axis_descriptor_t; typedef struct { diff --git a/primitive_process/src/base/primitive.cpp b/primitive_process/src/base/primitive.cpp index 51e455f..3a51423 100644 --- a/primitive_process/src/base/primitive.cpp +++ b/primitive_process/src/base/primitive.cpp @@ -42,6 +42,7 @@ void primitive::initialize(const std::vector(reversed_flags[i])); + // std::cout << matrices[i]->local_to_world.matrix() << "\n\n"; data_center->require_surface(subface_types[i], matrices[i], additional_structures[i], subfaces[i]); } @@ -70,7 +71,7 @@ void primitive::apply_transform(internal::transform_type type, Eigen::Vector4d p switch (type) { case internal::transform_type::scale: { - auto scale = Eigen::Isometry3d(Eigen::Scaling(param)); + auto scale = Eigen::Affine3d(Eigen::Scaling(param)); for (auto &model_matrix : new_model_matrices) { model_matrix.local_to_world = scale * model_matrix.local_to_world; model_matrix.world_to_local = model_matrix.world_to_local * scale.inverse(); diff --git a/primitive_process/src/primitive/axis/helixline_extrude_axis.cpp b/primitive_process/src/primitive/axis/helixline_extrude_axis.cpp index fb447e3..111c8b8 100644 --- a/primitive_process/src/primitive/axis/helixline_extrude_axis.cpp +++ b/primitive_process/src/primitive/axis/helixline_extrude_axis.cpp @@ -1,6 +1,7 @@ #include #include +#include namespace internal { @@ -17,44 +18,68 @@ double helixline_extrude_axis::get_closest_param(Eigen::Vector3d p) const const auto theta = std::abs(this->total_theta); // 等价转化为椭圆型Kepler方程，计算相关常量 - const auto e = r_p * theta * theta; - const auto M = pi + theta * h_p - theta_p; - const auto k_range = theta * inv_two_pi; - - // 每一周期内的计算函数 - auto halley_solve = [&](auto k) { - const auto M_k = M - k * two_pi; - - auto E = (M_k + e * pi) / (1 + e); - // 经实验验证，我们只需要一步Halley迭代就能产生足够的精度 - const auto K = E - e * std::sin(E); - const auto Kp = 1 - e * std::cos(E); - const auto Kpp = e * std::sin(E); - const auto delta_E = 2 * (K - M_k) * Kp / (2 * Kp * Kp - (K - M) * Kpp); - E -= delta_E; - - const auto t = (E + theta_p + k * two_pi - pi) / theta; - - if (t >= 0 && t <= 1) - return std::array{t, distance(t, theta, r_p, theta_p, h_p)}; - else - return std::array{.0, std::numeric_limits::max()}; + const auto e = r_p * theta * theta; + auto k = std::ceil((theta * h_p - theta_p) * inv_two_pi); + auto M = pi + theta * h_p - theta_p - 2 * pi * k; + + auto sign = 1.; + if (M < 0) { + M = -M; + sign = -1.; + } + + auto markley_solve = [&] { + M = std::fmod(M, two_pi); + if (M > pi) M = two_pi - pi; + if (M < 1e-15) return .0; + + auto alpha = (3 * pi * pi + 1.6 * pi * (pi - M) / (e + 1)) / (pi * pi - 6); + auto d = 3 * (1 - e) + alpha * e; + auto q = 2 * alpha * d * (1 - e) - M * M; + auto r = 3 * alpha * d * (d - 1 + e) * M + M * M * M; + + auto w = std::pow(std::abs(r) + std::sqrt(q * q * q + r * r), 2 / 3); + auto E1 = (2 * r * w / (w * w + w * q + q * q) + M) / d; + + auto sinE = std::sin(E1), cosE = std::cos(E1); + auto f = E1 - e * sinE - M; + auto f1 = 1 - e * cosE; + auto f2 = e * sinE; + auto f3 = e * cosE; + auto f4 = -e * sinE; + + auto denom3 = f1 - .5 * f * f2 / f1; + auto d3 = -f / denom3; + auto denom4 = f1 + .5 * d3 * f2 + (d3 * d3) * f3 / 6; + auto d4 = -f / denom4; + auto denom5 = f1 + .5 * d4 * f2 + (d4 * d4) * f3 / 6 + (d4 * d4 * d4) * f4 / 24; + auto d5 = -f / denom5; + + return E1 + d5; }; - std::vector k_indices(k_range); - std::iota(k_indices.begin(), k_indices.end(), 0); - auto res = std::transform_reduce( - k_indices.begin(), - k_indices.end(), - std::array{.0, std::numeric_limits::max()}, - [](auto&& lhs, auto&& rhs) { return lhs[1] < rhs[1] ? lhs : rhs; }, - [&](auto k) { return halley_solve(k); }); - - // 额外比较与区间端点处的值 - if (auto dist0 = distance(.0, theta, r_p, theta_p, h_p); dist0 < res[1]) res = {.0, dist0}; - if (auto dist1 = distance(1., theta, r_p, theta_p, h_p); dist1 < res[1]) res = {1., dist1}; + auto halley_solve = [&] { + auto E = (M + e * pi) / (e + 1); + auto delta = std::numeric_limits::max(); + size_t step{}; + while (delta >= epsilon && step < 8) { + auto sinE = std::sin(E), cosE = std::cos(E); + auto K = E - e * sinE - M; + auto Kp = 1 - e * cosE; + auto Kpp = e * sinE; + auto denom = 2 * Kp * Kp - K * Kpp; + auto dE = denom * 1e6 > epsilon ? (2 * K * Kp) / denom : K / Kp; + + E -= dE; + delta = std::abs(dE); + } + + return E; + }; - return res[0]; + auto E = e <= 1 ? markley_solve() : halley_solve(); + E *= sign; + return (theta_p + (2 * k - 1) * pi + E) / theta; } std::pair helixline_extrude_axis::to_pattern_local(double t, Eigen::Vector3d p) const @@ -81,16 +106,20 @@ std::pair helixline_extrude_axis::to_patte }; } -Eigen::AffineCompact3d helixline_extrude_axis::get_local_cap_mat(double t) const +Eigen::AffineCompact3d helixline_extrude_axis::get_local_cap_mat(double t, bool reversed) const { - const Eigen::Vector3d o = {std::cos(this->total_theta * t), std::sin(this->total_theta * t), t}; - Eigen::Vector3d tangent = {-std::sin(this->total_theta * t), std::cos(this->total_theta * t), 1}; - tangent /= sqrt_2; - - Eigen::AffineCompact3d res{}; - auto rot = Eigen::Quaterniond::FromTwoVectors(Eigen::Vector3d{1, 0, 0}, tangent); - res.linear() = rot.toRotationMatrix(); - res.translation() = o; + auto sine = std::sin(this->total_theta * t); + auto cosine = std::cos(this->total_theta * t); + const Eigen::Vector3d o = {cosine, sine, t}; + Eigen::Vector3d tangent = {this->total_theta * -sine, this->total_theta * cosine, 1}; + tangent.normalize(); + + Eigen::AffineCompact3d res = + reversed ? internal::trans_by_target_axis_matrices_ptr->local_to_world + : internal::trans_by_target_axis_matrices_ptr->local_to_world; + // auto rot = Eigen::Quaterniond::FromTwoVectors(Eigen::Vector3d{1, 0, 0}, reversed ? -tangent : tangent); + // res.linear() = rot.toRotationMatrix(); + res.translation() = o; return res; } diff --git a/primitive_process/src/primitive/axis/polyline_extrude_axis.cpp b/primitive_process/src/primitive/axis/polyline_extrude_axis.cpp index 202ff5e..d41792a 100644 --- a/primitive_process/src/primitive/axis/polyline_extrude_axis.cpp +++ b/primitive_process/src/primitive/axis/polyline_extrude_axis.cpp @@ -41,9 +41,9 @@ double polyline_extrude_axis::get_closest_param(Eigen::Vector3d p) const const auto r = ca.norm(); // 圆弧半径 if (cp_norm < epsilon) return std::array{r, .0, static_cast(index)}; - pa.normalize(); + ca /= r; cp /= cp_norm; - const auto delta_theta = pi - std::acos(pa.dot(cp)); + const auto delta_theta = std::acos(ca.dot(cp)); const auto theta = std::atan(bulge) * 2; const auto t = std::clamp(delta_theta / theta, .0, 1.); @@ -69,16 +69,17 @@ double polyline_extrude_axis::get_closest_param(Eigen::Vector3d p) const }, [&](size_t i) { return (std::abs(this->bulges[i]) <= epsilon) ? line_cpm(i) : circle_cpm(i); }); - return n + t; + return static_cast(n + t) / bulges.size(); } std::pair polyline_extrude_axis::to_pattern_local(double t, Eigen::Vector3d p) const { const Eigen::Vector2d local_p = p.head<2>(); - size_t n = static_cast(t); - double frac = t - n; - if (n == vertices.size()) { // 重定位到合法区间的终点 + auto t_ = t * bulges.size(); + size_t n = static_cast(t_); + double frac = t_ - n; + if (n == bulges.size()) { // 重定位到合法区间的终点 n -= 1; frac = 1.; } @@ -90,14 +91,17 @@ std::pair polyline_extrude_axis::to_patter const auto line = verts[1] - verts[0]; const auto tangent = line.normalized(); normal = {tangent.y(), -tangent.x()}; - o = verts[0] + t * line; + o = verts[0] + frac * line; } else { // 圆弧段 - const Eigen::Vector2d a = verts[0]; - const Eigen::Vector2d center = verts[1]; - const auto ca = a - center; - const auto theta = std::atan(bulge) * 2; - - Eigen::Rotation2Dd rot{t * theta}; + const auto a = verts[0]; + const auto b = verts[1]; + const auto half_ab = 0.5 * (a + b); + const auto bisec_vec = 0.5 * Eigen::Vector2d{a.y() - b.y(), b.x() - a.x()}; + const auto center = half_ab + bisec_vec / bulge; + const auto ca = a - center; + const auto theta = std::atan(bulge) * 2; + + Eigen::Rotation2Dd rot{frac * theta}; const auto cp = rot * ca; normal = theta < .0 ? cp.normalized() : -cp.normalized(); @@ -109,8 +113,8 @@ std::pair polyline_extrude_axis::to_patter const auto proj_x = -p.z(); // 因为bitangent永远为(0, 0, -1) Eigen::Projective2d pattern_to_axis{}; - pattern_to_axis.linear()(2, 0) = -1; - pattern_to_axis.linear().col(1).head<2>() = normal; + pattern_to_axis.matrix()(2, 0) = -1; + pattern_to_axis.matrix().col(1).head<2>() = normal; pattern_to_axis.translation().head<2>() = o; return { @@ -119,11 +123,12 @@ std::pair polyline_extrude_axis::to_patter }; } -Eigen::AffineCompact3d polyline_extrude_axis::get_local_cap_mat(double t) const +Eigen::AffineCompact3d polyline_extrude_axis::get_local_cap_mat(double t, bool reversed) const { - size_t n = static_cast(t); - double frac = t - n; - if (n == vertices.size()) { // 重定位到合法区间的终点 + auto t_ = t * bulges.size(); + size_t n = static_cast(t_); + double frac = t_ - n; + if (n == bulges.size()) { // 重定位到合法区间的终点 n -= 1; frac = 1.; } @@ -132,26 +137,31 @@ Eigen::AffineCompact3d polyline_extrude_axis::get_local_cap_mat(double t) const const auto bulge = bulges[n]; Eigen::Vector2d tangent{}, o{}; // 法线永远指向外侧 if (std::abs(bulge) <= epsilon) { // 直线 - const auto line = verts[1] - verts[0]; - const auto tangent = line.normalized(); - o = verts[0] + t * line; + const auto line = verts[1] - verts[0]; + tangent = line.normalized(); + o = verts[0] + t * line; } else { // 圆弧段 - const Eigen::Vector2d a = verts[0]; - const Eigen::Vector2d center = verts[1]; - const auto ca = a - center; - const auto theta = std::atan(bulge) * 2; - - Eigen::Rotation2Dd rot{t * theta}; + const auto a = verts[0]; + const auto b = verts[1]; + const auto half_ab = 0.5 * (a + b); + const auto bisec_vec = 0.5 * Eigen::Vector2d{a.y() - b.y(), b.x() - a.x()}; + const auto center = half_ab + bisec_vec / bulge; + const auto ca = a - center; + const auto theta = std::atan(bulge) * 2; + + Eigen::Rotation2Dd rot{frac * theta}; auto cp = rot * ca; o = cp + center; cp.normalize(); - tangent = theta < 0 ? Eigen::Vector2d{ca.y(), -ca.x()} : Eigen::Vector2d{-ca.y(), ca.x()}; + tangent = theta < 0 ? Eigen::Vector2d{cp.y(), -cp.x()} : Eigen::Vector2d{-cp.y(), cp.x()}; } Eigen::AffineCompact3d res{}; - auto rot = Eigen::Quaterniond::FromTwoVectors(Eigen::Vector3d{1, 0, 0}, Eigen::Vector3d{tangent[0], tangent[1], 0}); - res.linear() = rot.toRotationMatrix(); + auto rot = Eigen::Quaterniond::FromTwoVectors( + Eigen::Vector3d{1, 0, 0}, + reversed ? Eigen::Vector3d{-tangent[0], -tangent[1], 0} : Eigen::Vector3d{tangent[0], tangent[1], 0}); + res.linear() = rot.toRotationMatrix(); res.translation() = Eigen::Vector3d{o[0], o[1], 0}; return res; } diff --git a/primitive_process/src/primitive/pattern/polyline_pattern.cpp b/primitive_process/src/primitive/pattern/polyline_pattern.cpp index 65efb08..5971d3b 100644 --- a/primitive_process/src/primitive/pattern/polyline_pattern.cpp +++ b/primitive_process/src/primitive/pattern/polyline_pattern.cpp @@ -2,7 +2,7 @@ namespace internal { -std::pair polyline_pattern::eval_sdf(Eigen::Vector2d p) const +double polyline_pattern::eval_sdf(Eigen::Vector2d p) const { auto line_cpm = [&](auto index) { const auto a = this->vertices[index]; @@ -22,72 +22,163 @@ std::pair polyline_pattern::eval_sdf(Eigen::Vector2d p) c const auto cross = ap.x() * bp.y() - ap.y() * bp.x(); const auto d_theta = std::atan2(cross, dot); - return std::array{dist, d_theta, proj.x(), proj.y()}; + return std::array{dist, d_theta}; }; auto circle_cpm = [&](auto index) { - const auto bulge = this->bulges[index]; - const auto a = this->vertices[index]; - const auto b = this->vertices[index + 1]; - const auto half_ab = 0.5 * (a + b); - const auto bisec_vec = 0.5 * Eigen::Vector2d{a.y() - b.y(), b.x() - a.x()}; - const auto center = half_ab + bisec_vec / bulge; - Eigen::Vector2d ca = a - center; - Eigen::Vector2d cb = b - center; - Eigen::Vector2d cp = p - center; - Eigen::Vector2d pa = a - p; - Eigen::Vector2d pb = b - p; - - const auto D = cp.norm(); // 查询点到圆心的距离 - const auto r = ca.norm(); // 圆弧半径 - if (D < epsilon) return std::array{r, std::atan(bulge) * 2, a.x(), a.y()}; + const auto bulge = this->bulges[index]; + const auto a = this->vertices[index]; + const auto b = this->vertices[index + 1]; + const Eigen::Vector2d half_ab = 0.5 * (a + b); + const Eigen::Vector2d bisec_vec = 0.5 * Eigen::Vector2d{a.y() - b.y(), b.x() - a.x()}; + const Eigen::Vector2d center = half_ab + bisec_vec / bulge; + Eigen::Vector2d ca = a - center; + Eigen::Vector2d cb = b - center; + Eigen::Vector2d cp = p - center; + Eigen::Vector2d pa = a - p; + Eigen::Vector2d pb = b - p; + + const auto D = cp.norm(); // 查询点到圆心的距离 + const auto r = ca.norm(); // 圆弧半径 + const auto half_theta = std::atan(bulge); + if (D <= epsilon) return std::array{r, half_theta * 2}; // 计算GWN - const auto term1 = std::atan(bulge); - auto dot = ca.dot(cp); - const auto cross1 = cp.x() * ca.y() - ca.x() * cp.y(); - const auto term2 = ((r - D) * cross1) / ((r + D) * (r * D + dot)); - dot = cb.dot(cp); - const auto cross2 = cp.x() * cb.y() - cb.x() * cp.y(); - const auto term3 = ((r - D) * cross2) / ((r + D) * (r * D + dot)); - const auto gwn = term1 + term3 - term2; - - cp /= D; - ca /= r; - cb /= r; - - // 原理：根据叉乘符号判断方向，不论何种情况cp都应该在ca和cb之间 - if (cross1 * cross2 > 0) { - // 在圆弧范围内 - return std::array{D - r, gwn, center.x() + cp.x(), center.y() + cp.y()}; - } else { - // 最近点为圆弧端点 - if (auto pa_norm = pa.norm(), pb_norm = pb.norm(); pa_norm < pb_norm) - return std::array{pa_norm, gwn, a.x(), a.y()}; - else - return std::array{pa_norm, gwn, b.x(), b.y()}; + const auto gwn_main = half_theta; + const auto r_mul_D = r * D, r_plus_D = r + D, abs_r_sub_D = std::abs(r - D); + auto gwn = gwn_main; + if (abs_r_sub_D > epsilon) { + auto dot = ca.dot(cp); + auto cross = cp.x() * ca.y() - ca.x() * cp.y(); + auto Y = r_plus_D * cross, X = abs_r_sub_D * (r_mul_D + dot); + auto Ga = std::atan2(Y, X); + if (std::abs(X) <= epsilon && std::abs(Y) <= epsilon) Ga = std::copysign(pi_div_2, -half_theta); + dot = cb.dot(cp); + cross = cp.x() * cb.y() - cb.x() * cp.y(); + Y = r_plus_D * cross; + X = abs_r_sub_D * (r_mul_D + dot); + auto Gb = std::atan2(Y, X); + if (std::abs(X) <= epsilon && std::abs(Y) <= epsilon) Gb = std::copysign(pi_div_2, -half_theta); + auto gwn_fix = Gb - Ga; + if (gwn_fix * half_theta < 0) gwn_fix += std::copysign(pi, half_theta); + gwn += (r > D) ? gwn_fix : -gwn_fix; } + + return std::array{std::min({abs_r_sub_D, pa.norm(), pb.norm()}), gwn}; }; // 单次遍历 std::vector indices(this->bulges.size()); std::iota(indices.begin(), indices.end(), 0); - auto [min_dist, gwn, x, y] = std::transform_reduce( + auto [min_dist, gwn] = std::transform_reduce( indices.begin(), indices.end(), - std::array{std::numeric_limits::max(), .0, .0, .0}, + std::array{std::numeric_limits::max(), .0}, [](auto&& lhs, auto&& rhs) { - const auto gwn = lhs[1] + rhs[1]; - if (lhs[0] < rhs[0]) - return std::array{lhs[0], gwn, lhs[2], lhs[3]}; - else - return std::array{rhs[0], gwn, rhs[2], rhs[3]}; + const auto gwn = lhs[1] + rhs[1]; + const auto dist = std::min(lhs[0], rhs[0]); + return std::array{dist, gwn}; }, [&](size_t i) { return (std::abs(this->bulges[i]) <= epsilon) ? line_cpm(i) : circle_cpm(i); }); - return { - std::abs(gwn) < two_pi * epsilon, - Eigen::Vector2d{x, y} - }; + // if (p.norm() > 1 && std::abs(gwn) / two_pi != .0) abort(); + return std::abs(gwn) / two_pi < 1. ? min_dist : -min_dist; } + +// double polyline_pattern::eval_sdf(Eigen::Vector2d p) const +// { +// auto line_cpm = [&](auto index) { +// const auto a = this->vertices[index]; +// const auto b = this->vertices[index + 1]; +// const Eigen::Vector2d line_vec = b - a; + +// const Eigen::Vector2d ap = p - a; +// const Eigen::Vector2d bp = p - b; + +// // 计算投影参数 +// const auto t = std::clamp(line_vec.dot(ap) / line_vec.squaredNorm(), .0, 1.); +// const Eigen::Vector2d proj = a + t * line_vec; +// const auto dist = (p - proj).norm(); + +// // 计算GWN +// const auto dot = ap.dot(bp); +// const auto cross = ap.x() * bp.y() - ap.y() * bp.x(); +// const auto d_theta = std::atan2(cross, dot); + +// return std::array{dist, d_theta, proj.x(), proj.y()}; +// }; + +// auto circle_cpm = [&](auto index) { +// const auto bulge = this->bulges[index]; +// const auto a = this->vertices[index]; +// const auto b = this->vertices[index + 1]; +// const Eigen::Vector2d half_ab = 0.5 * (a + b); +// const Eigen::Vector2d bisec_vec = 0.5 * Eigen::Vector2d{a.y() - b.y(), b.x() - a.x()}; +// const Eigen::Vector2d center = half_ab + bisec_vec / bulge; +// Eigen::Vector2d ca = a - center; +// Eigen::Vector2d cb = b - center; +// Eigen::Vector2d cp = p - center; +// Eigen::Vector2d pa = a - p; +// Eigen::Vector2d pb = b - p; + +// const auto D = cp.norm(); // 查询点到圆心的距离 +// const auto r = ca.norm(); // 圆弧半径 +// const auto half_theta = std::atan(bulge); +// if (D < epsilon) return std::array{r, half_theta * 2, a.x(), a.y()}; + +// // 计算GWN +// const auto gwn_main = half_theta; +// const auto r_mul_D = r * D, r_plus_D = r + D, abs_r_sub_D = std::abs(r - D); +// const auto cross1 = cp.x() * ca.y() - ca.x() * cp.y(); +// const auto cross2 = cp.x() * cb.y() - cb.x() * cp.y(); +// auto gwn = gwn_main; +// if (abs_r_sub_D > epsilon) { +// auto dot = ca.dot(cp); +// const auto Ga = std::atan2(r_plus_D * cross1, abs_r_sub_D * (r_mul_D + dot)); +// dot = cb.dot(cp); +// const auto Gb = std::atan2(r_plus_D * cross2, abs_r_sub_D * (r_mul_D + dot)); +// auto gwn_fix = Gb - Ga; +// if (gwn_fix * half_theta < 0) gwn_fix += std::copysign(pi, half_theta); +// gwn += (r > D) ? gwn_fix : -gwn_fix; +// } + +// cp /= D; +// ca /= r; +// cb /= r; + +// // 原理：根据叉乘符号判断方向，不论何种情况cp都应该在ca和cb之间 +// if (cross1 * cross2 > 0) { +// // 在圆弧范围内 +// return std::array{D - r, gwn, center.x() + cp.x(), center.y() + cp.y()}; +// } else { +// // 最近点为圆弧端点 +// if (auto pa_norm = pa.norm(), pb_norm = pb.norm(); pa_norm < pb_norm) +// return std::array{pa_norm, gwn, a.x(), a.y()}; +// else +// return std::array{pa_norm, gwn, b.x(), b.y()}; +// } +// }; + +// // 单次遍历 +// std::vector indices(this->bulges.size()); +// std::iota(indices.begin(), indices.end(), 0); +// auto [min_dist, gwn, x, y] = std::transform_reduce( +// indices.begin(), +// indices.end(), +// std::array{std::numeric_limits::max(), .0, .0, .0}, +// [](auto&& lhs, auto&& rhs) { +// const auto gwn = lhs[1] + rhs[1]; +// if (lhs[0] < rhs[0]) +// return std::array{lhs[0], gwn, lhs[2], lhs[3]}; +// else +// return std::array{rhs[0], gwn, rhs[2], rhs[3]}; +// }, +// [&](size_t i) { return (std::abs(this->bulges[i]) <= epsilon) ? line_cpm(i) : circle_cpm(i); }); + +// return { +// // std::abs(gwn) < two_pi * epsilon, +// std::abs(gwn) == 0., +// Eigen::Vector2d{x, y} +// }; +// } } // namespace internal \ No newline at end of file diff --git a/primitive_process/src/subface/extrude/extrude_side_face.cpp b/primitive_process/src/subface/extrude/extrude_side_face.cpp index 57c6011..7ec1467 100644 --- a/primitive_process/src/subface/extrude/extrude_side_face.cpp +++ b/primitive_process/src/subface/extrude/extrude_side_face.cpp @@ -11,14 +11,10 @@ double extrude_side_function_impl::eval_sdf(pointer_wrapper object, Eig const auto v = axis_->get_closest_param(local_p.head<3>()); auto [pattern_p, local_TBN] = axis_->to_pattern_local(v, local_p.head<3>()); - pattern* pattern_ = reinterpret_cast(object->additional_structure[1]); - const auto [outside, closest_point] = pattern_->eval_sdf(pattern_p); + pattern* pattern_ = reinterpret_cast(object->additional_structure[1]); + const auto sdf = pattern_->eval_sdf(pattern_p); - auto local_q = local_TBN * Eigen::Vector3d{closest_point.x(), closest_point.y(), 0}; - auto vec_pq = local_p.head<3>() - local_q; - auto sdf = object->local_to_world_scale().cwiseProduct(vec_pq).norm(); - - return outside ? sdf : -sdf; + return sdf; } Eigen::Vector3d extrude_side_function_impl::eval_sdf_grad(pointer_wrapper object, Eigen::Vector3d p)