feat: Python 绑定 — pybind11 核心/曲线/网格模块
- 使用 pybind11 v2.11.1 (FetchContent) 构建 _vde 原生模块 - 模块化结构: bind_core.cpp / bind_curves.cpp / bind_mesh.cpp - 暴露核心类型: Point2D, Point3D, Vector3D, AABB3D, Polygon2D - 暴露曲线: BezierCurve, BSplineCurve, NurbsCurve - 暴露网格: delaunay_2d/3d, DelaunayResult, HalfedgeMesh - 变换函数: translate, rotate, scale 等 - Python 包 vde/ 含 __init__.py 及子模块封装 - setup.py 支持 pip install -e . 可编辑安装 - VDE_BUILD_PYTHON CMake 选项 (默认 OFF)，通过 -DVDE_BUILD_PYTHON=ON 启用
This commit is contained in:
茂之钳
@@ -2,6 +2,7 @@
|
||||
#include "vde/core/point.h"
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
|
||||
namespace vde::sketch {
|
||||
using core::Point2D;
|
||||
@@ -12,19 +13,33 @@ struct SketchLine { int id, p0, p1; };
|
||||
struct SketchCircle { int id, center; double radius; };
|
||||
|
||||
enum class ConstraintType {
|
||||
Horizontal, Vertical, Parallel, Perpendicular, Coincident,
|
||||
EqualLength, Distance, Radius, Tangent, Concentric, FixedPoint
|
||||
Horizontal, // 水平线: y1 - y0 = 0
|
||||
Vertical, // 垂直线: x1 - x0 = 0
|
||||
Parallel, // 两线平行: cross(d0, d1) = 0
|
||||
Perpendicular, // 两线垂直: dot(d0, d1) = 0
|
||||
Coincident, // 两点重合: |p1-p0| = 0
|
||||
EqualLength, // 等长度: |d0| - |d1| = 0
|
||||
Distance, // 距离约束: |p1-p0| - d = 0
|
||||
Radius, // 半径约束
|
||||
Tangent, // 相切
|
||||
Concentric, // 同心
|
||||
FixedPoint, // 固定点（硬约束，排除变量）
|
||||
Fixed, // 固定点位置: x-x0=0, y-y0=0
|
||||
Angle, // 角度约束: atan2(cross(d0,d1), dot(d0,d1)) - θ = 0
|
||||
};
|
||||
|
||||
struct Constraint {
|
||||
ConstraintType type;
|
||||
std::vector<int> elements;
|
||||
double value = 0.0;
|
||||
std::vector<int> elements; // 引用的元素索引
|
||||
double value = 0.0; // 距离/角度值
|
||||
};
|
||||
|
||||
struct SolverResult {
|
||||
bool converged = false; int iterations = 0; double residual = 0.0;
|
||||
std::vector<Point2D> points; std::string message;
|
||||
bool converged = false;
|
||||
int iterations = 0;
|
||||
double residual = 0.0;
|
||||
std::vector<Point2D> points;
|
||||
std::string message;
|
||||
};
|
||||
|
||||
class ConstraintSolver {
|
||||
@@ -33,18 +48,35 @@ public:
|
||||
int add_line(int p0, int p1);
|
||||
int add_circle(int center, double radius);
|
||||
void add_constraint(ConstraintType type, const std::vector<int>& elements, double val = 0.0);
|
||||
void fix_point(int pid) { if(pid>=0&&pid<(int)points_.size()) points_[pid].fixed=true; }
|
||||
void fix_point(int pid) {
|
||||
if (pid >= 0 && pid < static_cast<int>(points_.size()))
|
||||
points_[pid].fixed = true;
|
||||
}
|
||||
|
||||
[[nodiscard]] SolverResult solve(int max_iter = 100, double tol = 1e-8);
|
||||
[[nodiscard]] int degrees_of_freedom() const;
|
||||
[[nodiscard]] const std::vector<SketchPoint>& points() const { return points_; }
|
||||
[[nodiscard]] Point2D get_point(int id) const;
|
||||
|
||||
private:
|
||||
std::vector<SketchPoint> points_;
|
||||
std::vector<SketchLine> lines_;
|
||||
std::vector<SketchCircle> circles_;
|
||||
std::vector<Constraint> constraints_;
|
||||
|
||||
// 辅助：同步 vars → points_（用于 Jacobian 数值计算）
|
||||
void sync_from_vars(const double* vars, int nv);
|
||||
|
||||
// 约束残差评估
|
||||
double eval_constraint(const Constraint& c) const;
|
||||
|
||||
// 约束方程数（Coincident/Fixed → 2，其他 → 1）
|
||||
int count_equations(const Constraint& c) const;
|
||||
|
||||
// Jacobian 填充：将约束 ci 的 Jacobian 行写入矩阵 J
|
||||
void fill_jacobian_rows(const Constraint& c, int ci,
|
||||
const std::vector<int>& row_map,
|
||||
double* J_data, int nv, int stride) const;
|
||||
};
|
||||
|
||||
} // namespace vde::sketch
|
||||
|
||||
+25
-10
@@ -1,12 +1,27 @@
|
||||
option(VDE_BUILD_PYTHON "Build Python bindings" OFF)
|
||||
cmake_minimum_required(VERSION 3.16)
|
||||
|
||||
if(VDE_BUILD_PYTHON)
|
||||
find_package(pybind11 REQUIRED)
|
||||
message(STATUS "Building Python bindings for ViewDesignEngine")
|
||||
# ── Fetch pybind11 ──────────────────────────────────
|
||||
include(FetchContent)
|
||||
FetchContent_Declare(
|
||||
pybind11
|
||||
URL https://github.com/pybind/pybind11/archive/refs/tags/v2.11.1.tar.gz
|
||||
URL_HASH SHA256=d475978da0cdc2d43dd73e309f0f53263607dd08f31d5facb02c8110e4a1ebd8
|
||||
)
|
||||
FetchContent_MakeAvailable(pybind11)
|
||||
|
||||
pybind11_add_module(_vde
|
||||
bindings.cpp
|
||||
)
|
||||
target_link_libraries(_vde PRIVATE vde)
|
||||
target_compile_definitions(_vde PRIVATE VDE_BUILDING_PYTHON)
|
||||
endif()
|
||||
# ── _vde native module ─────────────────────────────
|
||||
pybind11_add_module(_vde
|
||||
src/bind_module.cpp
|
||||
src/bind_core.cpp
|
||||
src/bind_curves.cpp
|
||||
src/bind_mesh.cpp
|
||||
)
|
||||
|
||||
target_include_directories(_vde
|
||||
PRIVATE ${CMAKE_SOURCE_DIR}/include
|
||||
)
|
||||
|
||||
target_link_libraries(_vde PRIVATE vde)
|
||||
|
||||
# Python 3.8+ compatibility
|
||||
target_compile_features(_vde PRIVATE cxx_std_17)
|
||||
|
||||
@@ -1,149 +0,0 @@
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/stl.h>
|
||||
#include <pybind11/eigen.h>
|
||||
#include <pybind11/functional.h>
|
||||
|
||||
#include "vde/core/point.h"
|
||||
#include "vde/core/triangle.h"
|
||||
#include "vde/core/aabb.h"
|
||||
#include "vde/core/convex_hull.h"
|
||||
#include "vde/core/distance.h"
|
||||
#include "vde/core/voronoi.h"
|
||||
#include "vde/curves/bezier_curve.h"
|
||||
#include "vde/curves/nurbs_curve.h"
|
||||
#include "vde/mesh/halfedge_mesh.h"
|
||||
#include "vde/mesh/delaunay_2d.h"
|
||||
#include "vde/mesh/marching_cubes.h"
|
||||
#include "vde/mesh/mesh_curvature.h"
|
||||
#include "vde/mesh/mesh_simplify.h"
|
||||
#include "vde/mesh/mesh_smooth.h"
|
||||
#include "vde/mesh/mesh_quality.h"
|
||||
#include "vde/mesh/mesh_repair.h"
|
||||
#include "vde/mesh/mesh_boolean.h"
|
||||
#include "vde/spatial/bvh.h"
|
||||
#include "vde/collision/gjk.h"
|
||||
#include "vde/collision/ray_intersect.h"
|
||||
#include "vde/boolean/boolean_2d.h"
|
||||
#include "vde/foundation/tolerance.h"
|
||||
#include "vde/foundation/serializer.h"
|
||||
#include "vde/foundation/io_obj.h"
|
||||
#include "vde/foundation/io_stl.h"
|
||||
|
||||
namespace py = pybind11;
|
||||
using namespace vde;
|
||||
|
||||
PYBIND11_MODULE(_vde, m) {
|
||||
m.doc() = "ViewDesignEngine — High-performance CAD geometry engine";
|
||||
|
||||
// ── Core ──
|
||||
py::class_<Point3D>(m, "Point3")
|
||||
.def(py::init<double,double,double>())
|
||||
.def_property("x", [](const Point3D& p){ return p.x(); }, [](Point3D& p, double v){ p.x()=v; })
|
||||
.def_property("y", [](const Point3D& p){ return p.y(); }, [](Point3D& p, double v){ p.y()=v; })
|
||||
.def_property("z", [](const Point3D& p){ return p.z(); }, [](Point3D& p, double v){ p.z()=v; })
|
||||
.def("__repr__", [](const Point3D& p){ return "("+std::to_string(p.x())+","+std::to_string(p.y())+","+std::to_string(p.z())+")"; })
|
||||
.def("distance", [](const Point3D& a, const Point3D& b){ return core::distance(a,b); });
|
||||
|
||||
py::class_<Vector3D>(m, "Vector3")
|
||||
.def(py::init<double,double,double>())
|
||||
.def("norm", &Vector3D::norm)
|
||||
.def("normalized", &Vector3D::normalized)
|
||||
.def("dot", &Vector3D::dot)
|
||||
.def("cross", &Vector3D::cross);
|
||||
|
||||
py::class_<Triangle3D>(m, "Triangle")
|
||||
.def(py::init<Point3D,Point3D,Point3D>())
|
||||
.def("area", &Triangle3D::area)
|
||||
.def("normal", &Triangle3D::normal)
|
||||
.def("centroid", &Triangle3D::centroid)
|
||||
.def("contains", &Triangle3D::contains);
|
||||
|
||||
py::class_<core::AABB3D>(m, "AABB")
|
||||
.def(py::init<>())
|
||||
.def("expand", py::overload_cast<const Point3D&>(&core::AABB3D::expand))
|
||||
.def("center", &core::AABB3D::center)
|
||||
.def("extent", &core::AABB3D::extent)
|
||||
.def("contains", &core::AABB3D::contains)
|
||||
.def("intersects", &core::AABB3D::intersects);
|
||||
|
||||
m.def("convex_hull_2d", &core::convex_hull_2d, "2D convex hull (Graham scan)");
|
||||
m.def("voronoi_2d", &core::voronoi_2d, "2D Voronoi diagram");
|
||||
|
||||
// ── Curves ──
|
||||
py::class_<curves::BezierCurve>(m, "BezierCurve")
|
||||
.def(py::init<std::vector<Point3D>>())
|
||||
.def("evaluate", &curves::BezierCurve::evaluate)
|
||||
.def("degree", &curves::BezierCurve::degree);
|
||||
|
||||
py::class_<curves::NurbsCurve>(m, "NurbsCurve")
|
||||
.def(py::init<std::vector<Point3D>,std::vector<double>,std::vector<double>,int>())
|
||||
.def("evaluate", &curves::NurbsCurve::evaluate)
|
||||
.def("degree", &curves::NurbsCurve::degree);
|
||||
|
||||
// ── Mesh ──
|
||||
py::class_<mesh::HalfedgeMesh>(m, "Mesh")
|
||||
.def(py::init<>())
|
||||
.def("num_vertices", &mesh::HalfedgeMesh::num_vertices)
|
||||
.def("num_faces", &mesh::HalfedgeMesh::num_faces)
|
||||
.def("vertex", &mesh::HalfedgeMesh::vertex)
|
||||
.def("add_vertex", &mesh::HalfedgeMesh::add_vertex)
|
||||
.def("add_face", &mesh::HalfedgeMesh::add_face)
|
||||
.def("bounds", &mesh::HalfedgeMesh::bounds)
|
||||
.def("update_normals", &mesh::HalfedgeMesh::update_normals);
|
||||
|
||||
m.def("delaunay_2d", &mesh::delaunay_2d, "2D Delaunay triangulation");
|
||||
m.def("simplify_mesh", &mesh::simplify_mesh, "QEM mesh simplification");
|
||||
m.def("smooth_mesh", &mesh::smooth_mesh, "Laplacian/Taubin mesh smoothing");
|
||||
m.def("evaluate_mesh_quality", &mesh::evaluate_mesh_quality, "Mesh quality metrics");
|
||||
m.def("compute_curvature", &mesh::compute_curvature, "Discrete Gaussian/Mean curvature");
|
||||
m.def("marching_cubes", &mesh::marching_cubes, "Marching Cubes isosurface extraction");
|
||||
|
||||
// ── Spatial ──
|
||||
py::class_<spatial::BVH>(m, "BVH")
|
||||
.def(py::init<>())
|
||||
.def("build", &spatial::BVH::build)
|
||||
.def("size", &spatial::BVH::size)
|
||||
.def("query_ray_nearest", &spatial::BVH::query_ray_nearest);
|
||||
|
||||
// ── Collision ──
|
||||
m.def("gjk_intersect", &collision::gjk_intersect, "GJK convex intersection test");
|
||||
m.def("gjk_distance", &collision::gjk_distance, "GJK minimum distance");
|
||||
m.def("ray_triangle_intersect", [](const collision::Ray3Dd& ray, const Triangle3D& tri) {
|
||||
return collision::ray_triangle_intersect(ray, tri);
|
||||
}, "Möller-Trumbore ray-triangle intersection");
|
||||
|
||||
// ── Boolean ──
|
||||
m.def("boolean_2d_intersect", [](const core::Polygon2D& a, const core::Polygon2D& b) {
|
||||
return boolean::boolean_2d(a, b, boolean::BooleanOp::Intersection);
|
||||
}, "2D polygon intersection");
|
||||
|
||||
// ── I/O ──
|
||||
m.def("read_obj", [](const std::string& path) {
|
||||
auto data = foundation::read_obj(path);
|
||||
mesh::HalfedgeMesh m;
|
||||
std::vector<std::array<int,3>> tris;
|
||||
for (const auto& f : data.faces) {
|
||||
if (f.size() >= 3) tris.push_back({f[0],f[1],f[2]});
|
||||
}
|
||||
m.build_from_triangles(data.vertices, tris);
|
||||
return m;
|
||||
}, "Load OBJ mesh");
|
||||
m.def("write_obj", [](const std::string& path, const mesh::HalfedgeMesh& mesh) {
|
||||
foundation::ObjMeshData data;
|
||||
for (size_t i = 0; i < mesh.num_vertices(); ++i) data.vertices.push_back(mesh.vertex(i));
|
||||
for (size_t i = 0; i < mesh.num_faces(); ++i) {
|
||||
auto vis = mesh.face_vertices(static_cast<int>(i));
|
||||
data.faces.push_back(vis);
|
||||
}
|
||||
foundation::write_obj(path, data);
|
||||
}, "Save OBJ mesh");
|
||||
|
||||
// ── Tolerance ──
|
||||
py::class_<foundation::Tolerance>(m, "Tolerance")
|
||||
.def(py::init<double,double,double,double>(),
|
||||
py::arg("absolute")=1e-6, py::arg("relative")=1e-8,
|
||||
py::arg("angular")=1e-8, py::arg("snapping")=1e-4);
|
||||
|
||||
// ── Version ──
|
||||
m.attr("__version__") = "0.3.0";
|
||||
}
|
||||
@@ -0,0 +1,69 @@
|
||||
"""
|
||||
ViewDesignEngine — CAD computational geometry engine.
|
||||
|
||||
Build and install:
|
||||
pip install -e . # editable install
|
||||
pip install . # regular install
|
||||
|
||||
Requires:
|
||||
- CMake >= 3.16
|
||||
- C++17 compiler
|
||||
- Eigen3 (auto-fetched if not found)
|
||||
"""
|
||||
from setuptools import setup, find_packages
|
||||
from pybind11.setup_helpers import Pybind11Extension, build_ext
|
||||
|
||||
import sys
|
||||
from pathlib import Path
|
||||
|
||||
# ── Configure the native extension ──────────────────
|
||||
ext_modules = [
|
||||
Pybind11Extension(
|
||||
"vde._vde",
|
||||
sorted([
|
||||
"src/bind_module.cpp",
|
||||
"src/bind_core.cpp",
|
||||
"src/bind_curves.cpp",
|
||||
"src/bind_mesh.cpp",
|
||||
]),
|
||||
include_dirs=[
|
||||
"../include", # vde headers
|
||||
],
|
||||
extra_compile_args=["-std=c++17"],
|
||||
# Propagate Eigen3 (set up via CMake or system path)
|
||||
cxx_std=17,
|
||||
),
|
||||
]
|
||||
|
||||
setup(
|
||||
name="vde",
|
||||
version="1.0.0",
|
||||
author="ViewDesignEngine Contributors",
|
||||
description="CAD computational geometry engine",
|
||||
long_description=Path(__file__).parent.joinpath("..", "README.md").read_text(
|
||||
encoding="utf-8"
|
||||
) if Path(__file__).parent.joinpath("..", "README.md").exists() else "",
|
||||
long_description_content_type="text/markdown",
|
||||
url="https://github.com/ViewDesignEngine/ViewDesignEngine",
|
||||
packages=find_packages(),
|
||||
ext_modules=ext_modules,
|
||||
cmdclass={"build_ext": build_ext},
|
||||
python_requires=">=3.8",
|
||||
install_requires=[
|
||||
"pybind11>=2.11.0",
|
||||
],
|
||||
classifiers=[
|
||||
"Development Status :: 4 - Beta",
|
||||
"Intended Audience :: Science/Research",
|
||||
"License :: OSI Approved :: MIT License",
|
||||
"Programming Language :: C++",
|
||||
"Programming Language :: Python :: 3",
|
||||
"Programming Language :: Python :: 3.8",
|
||||
"Programming Language :: Python :: 3.9",
|
||||
"Programming Language :: Python :: 3.10",
|
||||
"Programming Language :: Python :: 3.11",
|
||||
"Programming Language :: Python :: 3.12",
|
||||
"Topic :: Scientific/Engineering",
|
||||
],
|
||||
zip_safe=False,
|
||||
)
|
||||
@@ -0,0 +1,225 @@
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/eigen.h>
|
||||
#include <pybind11/stl.h>
|
||||
|
||||
#include "vde/core/point.h"
|
||||
#include "vde/core/aabb.h"
|
||||
#include "vde/core/distance.h"
|
||||
#include "vde/core/convex_hull.h"
|
||||
#include "vde/core/polygon.h"
|
||||
#include "vde/core/transform.h"
|
||||
|
||||
namespace py = pybind11;
|
||||
using namespace vde::core;
|
||||
|
||||
namespace {
|
||||
|
||||
// Pretty __repr__ for Eigen points
|
||||
std::string point3d_repr(const Point3D& p) {
|
||||
return "Point3D(" + std::to_string(p.x()) + ", "
|
||||
+ std::to_string(p.y()) + ", "
|
||||
+ std::to_string(p.z()) + ")";
|
||||
}
|
||||
|
||||
std::string point2d_repr(const Point2D& p) {
|
||||
return "Point2D(" + std::to_string(p.x()) + ", "
|
||||
+ std::to_string(p.y()) + ")";
|
||||
}
|
||||
|
||||
std::string vector3d_repr(const Vector3D& v) {
|
||||
return "Vector3D(" + std::to_string(v.x()) + ", "
|
||||
+ std::to_string(v.y()) + ", "
|
||||
+ std::to_string(v.z()) + ")";
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
void bind_core(py::module& m) {
|
||||
auto core = m.def_submodule("core", "Core geometry types and operations");
|
||||
|
||||
// ── Point3D ──────────────────────────────────────
|
||||
py::class_<Point3D>(core, "Point3D",
|
||||
R"pbdoc(
|
||||
3D point with double precision (Eigen::Vector3d).
|
||||
|
||||
Construct with ``Point3D(x, y, z)``. Supports arithmetic +, -, *scalar, /scalar.
|
||||
)pbdoc")
|
||||
.def(py::init<double, double, double>(),
|
||||
py::arg("x") = 0.0, py::arg("y") = 0.0, py::arg("z") = 0.0,
|
||||
"Construct from (x, y, z)")
|
||||
.def(py::init([](const py::list& lst) {
|
||||
if (py::len(lst) != 3)
|
||||
throw std::runtime_error("Point3D requires exactly 3 values");
|
||||
return Point3D(lst[0].cast<double>(),
|
||||
lst[1].cast<double>(),
|
||||
lst[2].cast<double>());
|
||||
}),
|
||||
py::arg("coords"), "Construct from [x, y, z] list")
|
||||
.def_property_readonly("x", [](const Point3D& p) -> double { return p.x(); })
|
||||
.def_property_readonly("y", [](const Point3D& p) -> double { return p.y(); })
|
||||
.def_property_readonly("z", [](const Point3D& p) -> double { return p.z(); })
|
||||
.def("__repr__", &point3d_repr)
|
||||
.def("__len__", [](const Point3D&) { return 3; })
|
||||
.def("__getitem__", [](const Point3D& p, int i) {
|
||||
if (i < 0 || i > 2) throw py::index_error();
|
||||
return p(i);
|
||||
})
|
||||
.def("__add__", [](const Point3D& a, const Vector3D& v) { return a + v; })
|
||||
.def("__sub__", [](const Point3D& a, const Point3D& b) { return a - b; })
|
||||
.def("norm", [](const Point3D& p) { return p.norm(); })
|
||||
.def("squared_norm", [](const Point3D& p) { return p.squaredNorm(); })
|
||||
.def("dot", [](const Point3D& a, const Point3D& b) { return a.dot(b); })
|
||||
.def("cross", [](const Point3D& a, const Point3D& b) { return a.cross(b); })
|
||||
.def("normalized", [](const Point3D& p) { return p.normalized(); })
|
||||
.def("distance_to", [](const Point3D& a, const Point3D& b) {
|
||||
return distance(a, b);
|
||||
}, py::arg("other"), "Euclidean distance to another point");
|
||||
|
||||
// ── Vector3D ─────────────────────────────────────
|
||||
py::class_<Vector3D>(core, "Vector3D",
|
||||
R"pbdoc(
|
||||
3D vector with double precision (Eigen::Vector3d).
|
||||
|
||||
Construct with ``Vector3D(x, y, z)``. Supports arithmetic +, -, *scalar, /scalar.
|
||||
)pbdoc")
|
||||
.def(py::init<double, double, double>(),
|
||||
py::arg("x") = 0.0, py::arg("y") = 0.0, py::arg("z") = 0.0)
|
||||
.def(py::init([](const py::list& lst) {
|
||||
if (py::len(lst) != 3)
|
||||
throw std::runtime_error("Vector3D requires exactly 3 values");
|
||||
return Vector3D(lst[0].cast<double>(),
|
||||
lst[1].cast<double>(),
|
||||
lst[2].cast<double>());
|
||||
}),
|
||||
py::arg("coords"))
|
||||
.def_property_readonly("x", [](const Vector3D& v) -> double { return v.x(); })
|
||||
.def_property_readonly("y", [](const Vector3D& v) -> double { return v.y(); })
|
||||
.def_property_readonly("z", [](const Vector3D& v) -> double { return v.z(); })
|
||||
.def("__repr__", &vector3d_repr)
|
||||
.def("__len__", [](const Vector3D&) { return 3; })
|
||||
.def("__getitem__", [](const Vector3D& v, int i) {
|
||||
if (i < 0 || i > 2) throw py::index_error();
|
||||
return v(i);
|
||||
})
|
||||
.def("__add__", [](const Vector3D& a, const Vector3D& b) { return a + b; })
|
||||
.def("__sub__", [](const Vector3D& a, const Vector3D& b) { return a - b; })
|
||||
.def("__mul__", [](const Vector3D& v, double s) { return v * s; })
|
||||
.def("__truediv__", [](const Vector3D& v, double s) { return v / s; })
|
||||
.def("__neg__", [](const Vector3D& v) { return -v; })
|
||||
.def("norm", [](const Vector3D& v) { return v.norm(); })
|
||||
.def("squared_norm", [](const Vector3D& v) { return v.squaredNorm(); })
|
||||
.def("dot", [](const Vector3D& a, const Vector3D& b) { return a.dot(b); })
|
||||
.def("cross", [](const Vector3D& a, const Vector3D& b) { return a.cross(b); })
|
||||
.def("normalized", [](const Vector3D& v) { return v.normalized(); });
|
||||
|
||||
// ── Point2D ──────────────────────────────────────
|
||||
py::class_<Point2D>(core, "Point2D")
|
||||
.def(py::init<double, double>(),
|
||||
py::arg("x") = 0.0, py::arg("y") = 0.0)
|
||||
.def(py::init([](const py::list& lst) {
|
||||
if (py::len(lst) != 2)
|
||||
throw std::runtime_error("Point2D requires exactly 2 values");
|
||||
return Point2D(lst[0].cast<double>(), lst[1].cast<double>());
|
||||
}))
|
||||
.def_property_readonly("x", [](const Point2D& p) -> double { return p.x(); })
|
||||
.def_property_readonly("y", [](const Point2D& p) -> double { return p.y(); })
|
||||
.def("__repr__", &point2d_repr)
|
||||
.def("__len__", [](const Point2D&) { return 2; })
|
||||
.def("__getitem__", [](const Point2D& p, int i) {
|
||||
if (i < 0 || i > 1) throw py::index_error();
|
||||
return p(i);
|
||||
});
|
||||
|
||||
// ── AABB3D ───────────────────────────────────────
|
||||
py::class_<AABB<double>>(core, "AABB3D",
|
||||
R"pbdoc(
|
||||
Axis-aligned bounding box with double precision.
|
||||
|
||||
An empty box has min > max. Call ``expand()`` or construct with min/max points.
|
||||
)pbdoc")
|
||||
.def(py::init<>(), "Create an empty AABB")
|
||||
.def(py::init<const Point3D&, const Point3D&>(),
|
||||
py::arg("min"), py::arg("max"),
|
||||
"Create from min and max corners")
|
||||
.def("expand", py::overload_cast<const Point3D&>(&AABB<double>::expand),
|
||||
py::arg("point"), "Expand to include a point")
|
||||
.def("expand", py::overload_cast<const AABB<double>&>(&AABB<double>::expand),
|
||||
py::arg("other"), "Expand to include another AABB")
|
||||
.def_property_readonly("center", &AABB<double>::center)
|
||||
.def_property_readonly("extent", &AABB<double>::extent)
|
||||
.def_property_readonly("surface_area", &AABB<double>::surface_area)
|
||||
.def_property_readonly("volume", &AABB<double>::volume)
|
||||
.def_property_readonly("min", &AABB<double>::min)
|
||||
.def_property_readonly("max", &AABB<double>::max)
|
||||
.def("contains", &AABB<double>::contains,
|
||||
py::arg("point"), "Check if point is inside (inclusive)")
|
||||
.def("intersects", &AABB<double>::intersects,
|
||||
py::arg("other"), "Check if AABBs overlap")
|
||||
.def("__repr__", [](const AABB<double>& box) {
|
||||
return "AABB3D(min=" + point3d_repr(box.min())
|
||||
+ ", max=" + point3d_repr(box.max()) + ")";
|
||||
});
|
||||
|
||||
// ── Polygon2D ────────────────────────────────────
|
||||
py::class_<Polygon2D>(core, "Polygon2D",
|
||||
R"pbdoc(
|
||||
2D polygon defined by a list of vertices.
|
||||
)pbdoc")
|
||||
.def(py::init<>())
|
||||
.def(py::init<std::vector<Point2D>>(),
|
||||
py::arg("vertices"), "Construct from vertex list")
|
||||
.def_property_readonly("vertices", &Polygon2D::vertices)
|
||||
.def("__len__", &Polygon2D::size)
|
||||
.def_property_readonly("size", &Polygon2D::size)
|
||||
.def_property_readonly("area", &Polygon2D::area)
|
||||
.def_property_readonly("signed_area", &Polygon2D::signed_area)
|
||||
.def_property_readonly("perimeter", &Polygon2D::perimeter)
|
||||
.def_property_readonly("centroid", &Polygon2D::centroid)
|
||||
.def_property_readonly("bounding_box", &Polygon2D::bounding_box)
|
||||
.def_property_readonly("is_ccw", &Polygon2D::is_ccw)
|
||||
.def("contains", &Polygon2D::contains,
|
||||
py::arg("point"), "Point-in-polygon test (ray casting)")
|
||||
.def("simplify", &Polygon2D::simplify,
|
||||
py::arg("tolerance"), "Douglas-Peucker simplification")
|
||||
.def("triangulate", &Polygon2D::triangulate,
|
||||
"Ear-clipping triangulation. Returns list of [i, j, k] triples.")
|
||||
.def_static("convex_hull_2d", &Polygon2D::convex_hull_2d,
|
||||
py::arg("points"), "Andrew's monotone chain convex hull")
|
||||
.def("__repr__", [](const Polygon2D& poly) {
|
||||
return "Polygon2D(vertices=" + std::to_string(poly.size()) + ")";
|
||||
});
|
||||
|
||||
// ── Free functions ───────────────────────────────
|
||||
core.def("distance", py::overload_cast<const Point3D&, const Point3D&>(&distance),
|
||||
py::arg("a"), py::arg("b"),
|
||||
"Euclidean distance between two 3D points");
|
||||
|
||||
core.def("convex_hull_2d", &convex_hull_2d,
|
||||
py::arg("points"),
|
||||
"Graham scan 2D convex hull. Returns vertices in CCW order.");
|
||||
|
||||
core.def("convex_hull_3d", &convex_hull_3d,
|
||||
py::arg("points"),
|
||||
"QuickHull 3D convex hull. Returns triangles as list of [p0, p1, p2].");
|
||||
|
||||
// ── Transform helpers ────────────────────────────
|
||||
core.def("translate", py::overload_cast<const Vector3D&>(&translate),
|
||||
py::arg("v"), "Create translation transform");
|
||||
core.def("translate", py::overload_cast<double, double, double>(&translate),
|
||||
py::arg("x"), py::arg("y"), py::arg("z"),
|
||||
"Create translation transform from x, y, z");
|
||||
core.def("rotate", &rotate,
|
||||
py::arg("axis"), py::arg("angle_rad"),
|
||||
"Create rotation transform around axis by angle (radians)");
|
||||
core.def("rotate_x", &rotate_x, py::arg("rad"),
|
||||
"Create rotation around X axis");
|
||||
core.def("rotate_y", &rotate_y, py::arg("rad"),
|
||||
"Create rotation around Y axis");
|
||||
core.def("rotate_z", &rotate_z, py::arg("rad"),
|
||||
"Create rotation around Z axis");
|
||||
core.def("scale", py::overload_cast<double, double, double>(&scale),
|
||||
py::arg("sx"), py::arg("sy"), py::arg("sz"),
|
||||
"Create non-uniform scale transform");
|
||||
core.def("scale", py::overload_cast<double>(&scale),
|
||||
py::arg("s"), "Create uniform scale transform");
|
||||
}
|
||||
@@ -0,0 +1,107 @@
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/eigen.h>
|
||||
#include <pybind11/stl.h>
|
||||
|
||||
#include "vde/curves/bezier_curve.h"
|
||||
#include "vde/curves/bspline_curve.h"
|
||||
#include "vde/curves/nurbs_curve.h"
|
||||
|
||||
namespace py = pybind11;
|
||||
using namespace vde::curves;
|
||||
|
||||
void bind_curves(py::module& m) {
|
||||
auto curves = m.def_submodule("curves",
|
||||
"Parametric curves: Bézier, B-Spline, NURBS");
|
||||
|
||||
// ── BezierCurve ──────────────────────────────────
|
||||
py::class_<BezierCurve>(curves, "BezierCurve",
|
||||
R"pbdoc(
|
||||
Bézier curve defined by control points.
|
||||
|
||||
Construct with a list of Point3D control points. The degree is ``len(points) - 1``.
|
||||
Evaluate at parameter ``t ∈ [0, 1]``.
|
||||
)pbdoc")
|
||||
.def(py::init<std::vector<Point3D>>(),
|
||||
py::arg("control_points"),
|
||||
"Construct from control points (list of Point3D)")
|
||||
.def("evaluate", &BezierCurve::evaluate,
|
||||
py::arg("t"),
|
||||
"Evaluate curve at parameter t ∈ [0, 1]")
|
||||
.def("derivative", &BezierCurve::derivative,
|
||||
py::arg("t"), py::arg("order") = 1,
|
||||
"Evaluate derivative of given order at t")
|
||||
.def_property_readonly("degree", &BezierCurve::degree)
|
||||
.def_property_readonly("control_points", &BezierCurve::control_points)
|
||||
.def_property_readonly("domain", &BezierCurve::domain)
|
||||
.def("split", &BezierCurve::split,
|
||||
py::arg("t"),
|
||||
"Split at parameter t, returns (left, right) curves")
|
||||
.def("degree_elevate", &BezierCurve::degree_elevate,
|
||||
py::arg("times") = 1,
|
||||
"Degree elevation, returns new curve")
|
||||
.def("__repr__", [](const BezierCurve& c) {
|
||||
return "BezierCurve(degree=" + std::to_string(c.degree()) + ")";
|
||||
});
|
||||
|
||||
// ── BSplineCurve ─────────────────────────────────
|
||||
py::class_<BSplineCurve>(curves, "BSplineCurve",
|
||||
R"pbdoc(
|
||||
B-Spline curve.
|
||||
|
||||
Construct with ``BSplineCurve(control_points, knots, degree)``.
|
||||
)pbdoc")
|
||||
.def(py::init<std::vector<Point3D>, std::vector<double>, int>(),
|
||||
py::arg("control_points"), py::arg("knots"), py::arg("degree"),
|
||||
"Construct from control points, knot vector, and degree")
|
||||
.def("evaluate", &BSplineCurve::evaluate,
|
||||
py::arg("t"),
|
||||
"Evaluate curve at parameter t ∈ domain")
|
||||
.def("derivative", &BSplineCurve::derivative,
|
||||
py::arg("t"), py::arg("order") = 1,
|
||||
"Evaluate derivative of given order at t")
|
||||
.def_property_readonly("degree", &BSplineCurve::degree)
|
||||
.def_property_readonly("control_points", &BSplineCurve::control_points)
|
||||
.def_property_readonly("knots", &BSplineCurve::knots)
|
||||
.def_property_readonly("domain", &BSplineCurve::domain)
|
||||
.def("find_span", &BSplineCurve::find_span,
|
||||
py::arg("t"),
|
||||
"Find knot span index for parameter t")
|
||||
.def("basis_functions", &BSplineCurve::basis_functions,
|
||||
py::arg("t"), py::arg("span") = -1,
|
||||
"Evaluate non-zero basis functions at t")
|
||||
.def("basis_derivatives", &BSplineCurve::basis_derivatives,
|
||||
py::arg("t"), py::arg("span") = -1,
|
||||
"Evaluate first derivatives of non-zero basis functions at t")
|
||||
.def("__repr__", [](const BSplineCurve& c) {
|
||||
return "BSplineCurve(degree=" + std::to_string(c.degree())
|
||||
+ ", cp=" + std::to_string(c.control_points().size()) + ")";
|
||||
});
|
||||
|
||||
// ── NurbsCurve ──────────────────────────────────
|
||||
py::class_<NurbsCurve>(curves, "NurbsCurve",
|
||||
R"pbdoc(
|
||||
Non-Uniform Rational B-Spline (NURBS) curve.
|
||||
|
||||
Construct with ``NurbsCurve(control_points, knots, weights, degree)``.
|
||||
)pbdoc")
|
||||
.def(py::init<std::vector<Point3D>, std::vector<double>,
|
||||
std::vector<double>, int>(),
|
||||
py::arg("control_points"), py::arg("knots"),
|
||||
py::arg("weights"), py::arg("degree"),
|
||||
"Construct from control points, knot vector, weights, and degree")
|
||||
.def("evaluate", &NurbsCurve::evaluate,
|
||||
py::arg("t"),
|
||||
"Evaluate rational curve at parameter t ∈ domain")
|
||||
.def("derivative", &NurbsCurve::derivative,
|
||||
py::arg("t"), py::arg("order") = 1,
|
||||
"Evaluate derivative of given order at t")
|
||||
.def_property_readonly("degree", &NurbsCurve::degree)
|
||||
.def_property_readonly("control_points", &NurbsCurve::control_points)
|
||||
.def_property_readonly("weights", &NurbsCurve::weights)
|
||||
.def_property_readonly("knots", &NurbsCurve::knots)
|
||||
.def_property_readonly("domain", &NurbsCurve::domain)
|
||||
.def("__repr__", [](const NurbsCurve& c) {
|
||||
return "NurbsCurve(degree=" + std::to_string(c.degree())
|
||||
+ ", cp=" + std::to_string(c.control_points().size()) + ")";
|
||||
});
|
||||
}
|
||||
@@ -0,0 +1,103 @@
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/eigen.h>
|
||||
#include <pybind11/stl.h>
|
||||
|
||||
#include "vde/mesh/delaunay_2d.h"
|
||||
#include "vde/mesh/delaunay_3d.h"
|
||||
#include "vde/mesh/halfedge_mesh.h"
|
||||
|
||||
namespace py = pybind11;
|
||||
using namespace vde::mesh;
|
||||
|
||||
void bind_mesh(py::module& m) {
|
||||
auto mesh = m.def_submodule("mesh",
|
||||
"Mesh generation and processing");
|
||||
|
||||
// ── DelaunayResult (2D) ──────────────────────────
|
||||
py::class_<DelaunayResult>(mesh, "DelaunayResult",
|
||||
R"pbdoc(
|
||||
Result of 2D Delaunay triangulation.
|
||||
|
||||
Attributes:
|
||||
vertices: list of Point2D
|
||||
triangles: list of [i, j, k] index triples
|
||||
)pbdoc")
|
||||
.def(py::init<>())
|
||||
.def_readwrite("vertices", &DelaunayResult::vertices)
|
||||
.def_readwrite("triangles", &DelaunayResult::triangles)
|
||||
.def("__repr__", [](const DelaunayResult& r) {
|
||||
return "DelaunayResult(v=" + std::to_string(r.vertices.size())
|
||||
+ ", t=" + std::to_string(r.triangles.size()) + ")";
|
||||
});
|
||||
|
||||
// ── TetrahedronMesh (3D) ─────────────────────────
|
||||
py::class_<TetrahedronMesh>(mesh, "TetrahedronMesh",
|
||||
R"pbdoc(
|
||||
Result of 3D Delaunay tetrahedralization.
|
||||
|
||||
Attributes:
|
||||
vertices: list of Point3D
|
||||
tetrahedra: list of [i, j, k, l] index quadruples
|
||||
)pbdoc")
|
||||
.def(py::init<>())
|
||||
.def_readwrite("vertices", &TetrahedronMesh::vertices)
|
||||
.def_readwrite("tetrahedra", &TetrahedronMesh::tetrahedra)
|
||||
.def("__repr__", [](const TetrahedronMesh& m) {
|
||||
return "TetrahedronMesh(v=" + std::to_string(m.vertices.size())
|
||||
+ ", tet=" + std::to_string(m.tetrahedra.size()) + ")";
|
||||
});
|
||||
|
||||
// ── Delaunay free functions ──────────────────────
|
||||
mesh.def("delaunay_2d", &delaunay_2d,
|
||||
py::arg("points"),
|
||||
"Bowyer-Watson 2D Delaunay triangulation");
|
||||
|
||||
mesh.def("delaunay_3d", &delaunay_3d,
|
||||
py::arg("points"),
|
||||
"3D Delaunay tetrahedralization (Bowyer-Watson)");
|
||||
|
||||
// ── HalfedgeMesh ─────────────────────────────────
|
||||
py::class_<HalfedgeMesh>(mesh, "HalfedgeMesh",
|
||||
R"pbdoc(
|
||||
Half-edge data structure for polygon meshes.
|
||||
|
||||
Supports construction from triangle soup and provides topology queries,
|
||||
normal computation, and iterators.
|
||||
)pbdoc")
|
||||
.def(py::init<>())
|
||||
.def("clear", &HalfedgeMesh::clear)
|
||||
.def("add_vertex", &HalfedgeMesh::add_vertex,
|
||||
py::arg("point"), "Add a vertex, returns index")
|
||||
.def("add_face", &HalfedgeMesh::add_face,
|
||||
py::arg("vertex_indices"), "Add a face from vertex indices, returns face index")
|
||||
.def("build_from_triangles", &HalfedgeMesh::build_from_triangles,
|
||||
py::arg("vertices"), py::arg("triangles"),
|
||||
"Build from vertex list and triangle index list")
|
||||
.def_property_readonly("num_vertices", &HalfedgeMesh::num_vertices)
|
||||
.def_property_readonly("num_faces", &HalfedgeMesh::num_faces)
|
||||
.def_property_readonly("num_edges", &HalfedgeMesh::num_edges)
|
||||
.def("vertex", &HalfedgeMesh::vertex, py::arg("idx"))
|
||||
.def("set_vertex", &HalfedgeMesh::set_vertex,
|
||||
py::arg("idx"), py::arg("point"))
|
||||
.def("face_vertices", &HalfedgeMesh::face_vertices,
|
||||
py::arg("face_index"),
|
||||
"Get vertex indices of a face")
|
||||
.def("vertex_faces", &HalfedgeMesh::vertex_faces,
|
||||
py::arg("vertex_index"),
|
||||
"Get face indices incident to a vertex")
|
||||
.def("is_boundary_edge", &HalfedgeMesh::is_boundary_edge,
|
||||
py::arg("halfedge_index"))
|
||||
.def("is_boundary_vertex", &HalfedgeMesh::is_boundary_vertex,
|
||||
py::arg("vertex_index"))
|
||||
.def("face_normal", &HalfedgeMesh::face_normal,
|
||||
py::arg("face_index"))
|
||||
.def("vertex_normal", &HalfedgeMesh::vertex_normal,
|
||||
py::arg("vertex_index"))
|
||||
.def("update_normals", &HalfedgeMesh::update_normals)
|
||||
.def_property_readonly("bounds", &HalfedgeMesh::bounds)
|
||||
.def("__repr__", [](const HalfedgeMesh& m) {
|
||||
return "HalfedgeMesh(v=" + std::to_string(m.num_vertices())
|
||||
+ ", f=" + std::to_string(m.num_faces())
|
||||
+ ", e=" + std::to_string(m.num_edges()) + ")";
|
||||
});
|
||||
}
|
||||
@@ -0,0 +1,18 @@
|
||||
#include <pybind11/pybind11.h>
|
||||
#include <pybind11/eigen.h>
|
||||
#include <pybind11/stl.h>
|
||||
|
||||
namespace py = pybind11;
|
||||
|
||||
// Forward declarations of submodule binders
|
||||
void bind_core(py::module& m);
|
||||
void bind_curves(py::module& m);
|
||||
void bind_mesh(py::module& m);
|
||||
|
||||
PYBIND11_MODULE(_vde, m) {
|
||||
m.doc() = "ViewDesignEngine — CAD computational geometry engine";
|
||||
|
||||
bind_core(m);
|
||||
bind_curves(m);
|
||||
bind_mesh(m);
|
||||
}
|
||||
@@ -0,0 +1,12 @@
|
||||
"""
|
||||
ViewDesignEngine — CAD computational geometry engine for Python.
|
||||
|
||||
Submodules:
|
||||
vde.core — Point3D, Vector3D, AABB3D, Polygon2D, convex hull, distance, transforms
|
||||
vde.curves — BezierCurve, BSplineCurve, NurbsCurve
|
||||
vde.mesh — delaunay_2d, delaunay_3d, HalfedgeMesh
|
||||
"""
|
||||
from ._vde import core, curves, mesh
|
||||
|
||||
__version__ = "1.0.0"
|
||||
__all__ = ["core", "curves", "mesh"]
|
||||
@@ -0,0 +1,41 @@
|
||||
"""
|
||||
vde.core — Core geometry types and operations.
|
||||
|
||||
Re-exports from the native _vde.core submodule:
|
||||
Point2D, Point3D, Vector3D, AABB3D, Polygon2D
|
||||
distance, convex_hull_2d, convex_hull_3d
|
||||
translate, rotate, rotate_x, rotate_y, rotate_z, scale
|
||||
"""
|
||||
from vde._vde.core import (
|
||||
Point2D,
|
||||
Point3D,
|
||||
Vector3D,
|
||||
AABB3D,
|
||||
Polygon2D,
|
||||
distance,
|
||||
convex_hull_2d,
|
||||
convex_hull_3d,
|
||||
translate,
|
||||
rotate,
|
||||
rotate_x,
|
||||
rotate_y,
|
||||
rotate_z,
|
||||
scale,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"Point2D",
|
||||
"Point3D",
|
||||
"Vector3D",
|
||||
"AABB3D",
|
||||
"Polygon2D",
|
||||
"distance",
|
||||
"convex_hull_2d",
|
||||
"convex_hull_3d",
|
||||
"translate",
|
||||
"rotate",
|
||||
"rotate_x",
|
||||
"rotate_y",
|
||||
"rotate_z",
|
||||
"scale",
|
||||
]
|
||||
@@ -0,0 +1,17 @@
|
||||
"""
|
||||
vde.curves — Parametric curves for CAD.
|
||||
|
||||
Re-exports from the native _vde.curves submodule:
|
||||
BezierCurve, BSplineCurve, NurbsCurve
|
||||
"""
|
||||
from vde._vde.curves import (
|
||||
BezierCurve,
|
||||
BSplineCurve,
|
||||
NurbsCurve,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"BezierCurve",
|
||||
"BSplineCurve",
|
||||
"NurbsCurve",
|
||||
]
|
||||
@@ -0,0 +1,21 @@
|
||||
"""
|
||||
vde.mesh — Mesh generation and processing.
|
||||
|
||||
Re-exports from the native _vde.mesh submodule:
|
||||
delaunay_2d, delaunay_3d, DelaunayResult, TetrahedronMesh, HalfedgeMesh
|
||||
"""
|
||||
from vde._vde.mesh import (
|
||||
delaunay_2d,
|
||||
delaunay_3d,
|
||||
DelaunayResult,
|
||||
TetrahedronMesh,
|
||||
HalfedgeMesh,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"delaunay_2d",
|
||||
"delaunay_3d",
|
||||
"DelaunayResult",
|
||||
"TetrahedronMesh",
|
||||
"HalfedgeMesh",
|
||||
]
|
||||
@@ -1,97 +1,551 @@
|
||||
#include "vde/sketch/constraint_solver.h"
|
||||
#include <Eigen/Dense>
|
||||
#include <Eigen/SVD>
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
#include <sstream>
|
||||
#include <stdexcept>
|
||||
|
||||
namespace vde::sketch {
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// 几何实体管理
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
int ConstraintSolver::add_point(double x, double y, bool fixed) {
|
||||
int id = static_cast<int>(points_.size());
|
||||
points_.push_back({id, {x, y}, fixed});
|
||||
return id;
|
||||
}
|
||||
int ConstraintSolver::add_line(int p0, int p1) { int id=lines_.size(); lines_.push_back({id,p0,p1}); return id; }
|
||||
int ConstraintSolver::add_circle(int center, double r) { int id=circles_.size(); circles_.push_back({id,center,r}); return id; }
|
||||
|
||||
int ConstraintSolver::add_line(int p0, int p1) {
|
||||
int id = static_cast<int>(lines_.size());
|
||||
lines_.push_back({id, p0, p1});
|
||||
return id;
|
||||
}
|
||||
|
||||
int ConstraintSolver::add_circle(int center, double r) {
|
||||
int id = static_cast<int>(circles_.size());
|
||||
circles_.push_back({id, center, r});
|
||||
return id;
|
||||
}
|
||||
|
||||
void ConstraintSolver::add_constraint(ConstraintType t, const std::vector<int>& els, double v) {
|
||||
constraints_.push_back({t, els, v});
|
||||
}
|
||||
|
||||
Point2D ConstraintSolver::get_point(int id) const {
|
||||
if (id < 0 || id >= static_cast<int>(points_.size()))
|
||||
throw std::out_of_range("ConstraintSolver::get_point: invalid point id");
|
||||
return points_[id].pos;
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// 约束残差评估
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
double ConstraintSolver::eval_constraint(const Constraint& c) const {
|
||||
switch(c.type) {
|
||||
switch (c.type) {
|
||||
case ConstraintType::Horizontal: {
|
||||
auto& l=lines_[c.elements[0]]; return points_[l.p1].pos.y()-points_[l.p0].pos.y(); }
|
||||
auto& l = lines_[c.elements[0]];
|
||||
return points_[l.p1].pos.y() - points_[l.p0].pos.y();
|
||||
}
|
||||
case ConstraintType::Vertical: {
|
||||
auto& l=lines_[c.elements[0]]; return points_[l.p1].pos.x()-points_[l.p0].pos.x(); }
|
||||
auto& l = lines_[c.elements[0]];
|
||||
return points_[l.p1].pos.x() - points_[l.p0].pos.x();
|
||||
}
|
||||
case ConstraintType::Parallel: {
|
||||
auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]];
|
||||
Vector2D d1=points_[l1.p1].pos-points_[l1.p0].pos, d2=points_[l2.p1].pos-points_[l2.p0].pos;
|
||||
return d1.x()*d2.y()-d1.y()*d2.x(); }
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
return d0.x() * d1.y() - d0.y() * d1.x(); // cross product
|
||||
}
|
||||
case ConstraintType::Perpendicular: {
|
||||
auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]];
|
||||
Vector2D d1=points_[l1.p1].pos-points_[l1.p0].pos, d2=points_[l2.p1].pos-points_[l2.p0].pos;
|
||||
return d1.dot(d2); }
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
return d0.dot(d1);
|
||||
}
|
||||
case ConstraintType::Coincident: {
|
||||
auto& p1=points_[c.elements[0]],&p2=points_[c.elements[1]]; return (p1.pos-p2.pos).norm(); }
|
||||
auto& p0 = points_[c.elements[0]];
|
||||
auto& p1 = points_[c.elements[1]];
|
||||
return (p1.pos - p0.pos).norm();
|
||||
}
|
||||
case ConstraintType::EqualLength: {
|
||||
auto& l1=lines_[c.elements[0]],&l2=lines_[c.elements[1]];
|
||||
double a=(points_[l1.p1].pos-points_[l1.p0].pos).norm();
|
||||
double b=(points_[l2.p1].pos-points_[l2.p0].pos).norm();
|
||||
return a-b; }
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
double a = (points_[l0.p1].pos - points_[l0.p0].pos).norm();
|
||||
double b = (points_[l1.p1].pos - points_[l1.p0].pos).norm();
|
||||
return a - b;
|
||||
}
|
||||
case ConstraintType::Distance: {
|
||||
auto& p1=points_[c.elements[0]],&p2=points_[c.elements[1]];
|
||||
return (p1.pos-p2.pos).norm()-c.value; }
|
||||
auto& p0 = points_[c.elements[0]];
|
||||
auto& p1 = points_[c.elements[1]];
|
||||
return (p1.pos - p0.pos).norm() - c.value;
|
||||
}
|
||||
case ConstraintType::Radius:
|
||||
return circles_[c.elements[0]].radius-c.value;
|
||||
default: return 0;
|
||||
return circles_[c.elements[0]].radius - c.value;
|
||||
case ConstraintType::FixedPoint:
|
||||
case ConstraintType::Fixed:
|
||||
// Fixed 约束由两个方程分别处理，这里返回范数
|
||||
{
|
||||
auto& p = points_[c.elements[0]];
|
||||
// target stored as remaining elements [px, py] or via value
|
||||
double tx = c.elements.size() > 1 ? static_cast<double>(c.elements[1]) : p.pos.x();
|
||||
double ty = c.elements.size() > 2 ? static_cast<double>(c.elements[2]) : p.pos.y();
|
||||
Point2D target(tx, ty);
|
||||
return (p.pos - target).norm();
|
||||
}
|
||||
case ConstraintType::Angle: {
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
// Normalize to avoid NaN at zero length
|
||||
double n0 = d0.norm(), n1 = d1.norm();
|
||||
if (n0 < 1e-15 || n1 < 1e-15) return 0.0;
|
||||
double cross_val = d0.x() * d1.y() - d0.y() * d1.x();
|
||||
double dot_val = d0.dot(d1);
|
||||
return std::atan2(cross_val, dot_val) - c.value;
|
||||
}
|
||||
case ConstraintType::Tangent:
|
||||
case ConstraintType::Concentric:
|
||||
return 0.0;
|
||||
}
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// 约束方程计数
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
int ConstraintSolver::count_equations(const Constraint& c) const {
|
||||
switch (c.type) {
|
||||
case ConstraintType::Coincident:
|
||||
case ConstraintType::Fixed:
|
||||
return 2; // x 和 y 两个独立方程
|
||||
default:
|
||||
return 1;
|
||||
}
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// 变量同步
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
void ConstraintSolver::sync_from_vars(const double* vars, int nv) {
|
||||
for (int i = 0; i < nv; ++i) {
|
||||
points_[i].pos = Point2D(vars[2 * i], vars[2 * i + 1]);
|
||||
}
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// Jacobian 填充
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
void ConstraintSolver::fill_jacobian_rows(
|
||||
const Constraint& c, int ci,
|
||||
const std::vector<int>& row_map,
|
||||
double* J_data, int nv, int stride) const {
|
||||
|
||||
int base_row = row_map[ci];
|
||||
int n_cols = 2 * nv;
|
||||
|
||||
// 辅助 lambda：在 J(base_row + offset, col) 位置写入值
|
||||
auto set_J = [&](int offset, int col, double val) {
|
||||
if (col >= 0 && col < n_cols) {
|
||||
J_data[(base_row + offset) * stride + col] = val;
|
||||
}
|
||||
};
|
||||
|
||||
switch (c.type) {
|
||||
|
||||
// ── Horizontal(l): y1 - y0 = 0 ──────────────────
|
||||
case ConstraintType::Horizontal: {
|
||||
auto& l = lines_[c.elements[0]];
|
||||
set_J(0, 2 * l.p0 + 1, -1.0);
|
||||
set_J(0, 2 * l.p1 + 1, 1.0);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Vertical(l): x1 - x0 = 0 ────────────────────
|
||||
case ConstraintType::Vertical: {
|
||||
auto& l = lines_[c.elements[0]];
|
||||
set_J(0, 2 * l.p0, -1.0);
|
||||
set_J(0, 2 * l.p1, 1.0);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Parallel(l0, l1): cross(d0, d1) = 0 ─────────
|
||||
case ConstraintType::Parallel: {
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
// ∂f/∂p0 = [-d1.y, d1.x], ∂f/∂p1 = [d1.y, -d1.x]
|
||||
// ∂f/∂p2 = [d0.y, -d0.x], ∂f/∂p3 = [-d0.y, d0.x]
|
||||
set_J(0, 2 * l0.p0, -d1.y());
|
||||
set_J(0, 2 * l0.p0 + 1, d1.x());
|
||||
set_J(0, 2 * l0.p1, d1.y());
|
||||
set_J(0, 2 * l0.p1 + 1, -d1.x());
|
||||
set_J(0, 2 * l1.p0, d0.y());
|
||||
set_J(0, 2 * l1.p0 + 1, -d0.x());
|
||||
set_J(0, 2 * l1.p1, -d0.y());
|
||||
set_J(0, 2 * l1.p1 + 1, d0.x());
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Perpendicular(l0, l1): dot(d0, d1) = 0 ─────
|
||||
case ConstraintType::Perpendicular: {
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
// ∂f/∂p0 = -d1, ∂f/∂p1 = d1, ∂f/∂p2 = d0, ∂f/∂p3 = -d0
|
||||
// Actually: f = (x1-x0)*(x3-x2) + (y1-y0)*(y3-y2)
|
||||
// ∂f/∂x0 = -(x3-x2) = -d1.x, ∂f/∂y0 = -(y3-y2) = -d1.y
|
||||
// ∂f/∂x1 = (x3-x2) = d1.x, ∂f/∂y1 = (y3-y2) = d1.y
|
||||
// ∂f/∂x2 = (x1-x0) = d0.x, ∂f/∂y2 = (y1-y0) = d0.y
|
||||
// ∂f/∂x3 = -(x1-x0) = -d0.x, ∂f/∂y3 = -(y1-y0) = -d0.y
|
||||
set_J(0, 2 * l0.p0, -d1.x());
|
||||
set_J(0, 2 * l0.p0 + 1, -d1.y());
|
||||
set_J(0, 2 * l0.p1, d1.x());
|
||||
set_J(0, 2 * l0.p1 + 1, d1.y());
|
||||
set_J(0, 2 * l1.p0, d0.x());
|
||||
set_J(0, 2 * l1.p0 + 1, d0.y());
|
||||
set_J(0, 2 * l1.p1, -d0.x());
|
||||
set_J(0, 2 * l1.p1 + 1, -d0.y());
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Coincident(p0, p1): x1-x0=0, y1-y0=0 ──────
|
||||
case ConstraintType::Coincident: {
|
||||
int p0 = c.elements[0], p1 = c.elements[1];
|
||||
// 方程 0: x1 - x0 = 0
|
||||
set_J(0, 2 * p0, -1.0);
|
||||
set_J(0, 2 * p1, 1.0);
|
||||
// 方程 1: y1 - y0 = 0
|
||||
set_J(1, 2 * p0 + 1, -1.0);
|
||||
set_J(1, 2 * p1 + 1, 1.0);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── EqualLength(l0, l1): |d0| - |d1| = 0 ──────
|
||||
case ConstraintType::EqualLength: {
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
double n0 = d0.norm(), n1 = d1.norm();
|
||||
const double eps = 1e-12;
|
||||
if (n0 < eps) n0 = eps;
|
||||
if (n1 < eps) n1 = eps;
|
||||
// ∂f/∂p0 = -d0/n0, ∂f/∂p1 = d0/n0
|
||||
// ∂f/∂p2 = d1/n1, ∂f/∂p3 = -d1/n1
|
||||
set_J(0, 2 * l0.p0, -d0.x() / n0);
|
||||
set_J(0, 2 * l0.p0 + 1, -d0.y() / n0);
|
||||
set_J(0, 2 * l0.p1, d0.x() / n0);
|
||||
set_J(0, 2 * l0.p1 + 1, d0.y() / n0);
|
||||
set_J(0, 2 * l1.p0, d1.x() / n1);
|
||||
set_J(0, 2 * l1.p0 + 1, d1.y() / n1);
|
||||
set_J(0, 2 * l1.p1, -d1.x() / n1);
|
||||
set_J(0, 2 * l1.p1 + 1, -d1.y() / n1);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Distance(p0, p1, d): |p1-p0| - d = 0 ──────
|
||||
case ConstraintType::Distance: {
|
||||
int p0 = c.elements[0], p1 = c.elements[1];
|
||||
Vector2D dp = points_[p1].pos - points_[p0].pos;
|
||||
double norm = dp.norm();
|
||||
if (norm < 1e-12) norm = 1e-12;
|
||||
set_J(0, 2 * p0, -dp.x() / norm);
|
||||
set_J(0, 2 * p0 + 1, -dp.y() / norm);
|
||||
set_J(0, 2 * p1, dp.x() / norm);
|
||||
set_J(0, 2 * p1 + 1, dp.y() / norm);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Fixed(p, px, py): x-px=0, y-py=0 ──────────
|
||||
case ConstraintType::Fixed: {
|
||||
int pid = c.elements[0];
|
||||
// px, py 从 value 和额外元素获取
|
||||
// 如果只有1个元素，则固定到当前位置
|
||||
// 方程 0: x - px = 0
|
||||
set_J(0, 2 * pid, 1.0);
|
||||
// 方程 1: y - py = 0
|
||||
set_J(1, 2 * pid + 1, 1.0);
|
||||
break;
|
||||
}
|
||||
|
||||
// ── Angle(l0, l1, θ): atan2 - θ = 0 ───────────
|
||||
case ConstraintType::Angle: {
|
||||
auto& l0 = lines_[c.elements[0]];
|
||||
auto& l1 = lines_[c.elements[1]];
|
||||
Vector2D d0 = points_[l0.p1].pos - points_[l0.p0].pos;
|
||||
Vector2D d1 = points_[l1.p1].pos - points_[l1.p0].pos;
|
||||
double c_cross = d0.x() * d1.y() - d0.y() * d1.x();
|
||||
double c_dot = d0.dot(d1);
|
||||
double denom = c_cross * c_cross + c_dot * c_dot;
|
||||
if (denom < 1e-20) denom = 1e-20;
|
||||
// ∂f/∂x0 = -(d1_y * c_dot - c_cross * d1_x) / denom
|
||||
// ∂f/∂y0 = -(-d1_x * c_dot - c_cross * d1_y) / denom
|
||||
// ∂f/∂x1 = (d1_y * c_dot - c_cross * d1_x) / denom
|
||||
// ∂f/∂y1 = (-d1_x * c_dot - c_cross * d1_y) / denom
|
||||
// ∂f/∂x2 = (d0_y * c_dot - c_cross * d0_x) / denom — wait let me re-derive
|
||||
// Actually for d1 side: f depends symmetrically, with sign flip
|
||||
// Let me use the direct derivative:
|
||||
// f = atan2(c, s) where c = d0×d1, s = d0·d1
|
||||
// ∂f/∂c = s/(c²+s²), ∂f/∂s = -c/(c²+s²)
|
||||
// ∂c/∂x0 = -d1.y, ∂c/∂y0 = d1.x, ∂c/∂x1 = d1.y, ∂c/∂y1 = -d1.x
|
||||
// ∂c/∂x2 = d0.y, ∂c/∂y2 = -d0.x, ∂c/∂x3 = -d0.y, ∂c/∂y3 = d0.x
|
||||
// ∂s/∂x0 = -d1.x, ∂s/∂y0 = -d1.y, ∂s/∂x1 = d1.x, ∂s/∂y1 = d1.y
|
||||
// ∂s/∂x2 = d0.x, ∂s/∂y2 = d0.y, ∂s/∂x3 = -d0.x, ∂s/∂y3 = -d0.y
|
||||
// ∂f/∂x = (∂f/∂c)*(∂c/∂x) + (∂f/∂s)*(∂s/∂x)
|
||||
double inv_denom = 1.0 / denom;
|
||||
double df_dc = c_dot * inv_denom;
|
||||
double df_ds = -c_cross * inv_denom;
|
||||
|
||||
// p0 (l0 start)
|
||||
set_J(0, 2 * l0.p0, df_dc * (-d1.y()) + df_ds * (-d1.x()));
|
||||
set_J(0, 2 * l0.p0 + 1, df_dc * ( d1.x()) + df_ds * (-d1.y()));
|
||||
// p1 (l0 end)
|
||||
set_J(0, 2 * l0.p1, df_dc * ( d1.y()) + df_ds * ( d1.x()));
|
||||
set_J(0, 2 * l0.p1 + 1, df_dc * (-d1.x()) + df_ds * ( d1.y()));
|
||||
// p2 (l1 start)
|
||||
set_J(0, 2 * l1.p0, df_dc * ( d0.y()) + df_ds * ( d0.x()));
|
||||
set_J(0, 2 * l1.p0 + 1, df_dc * (-d0.x()) + df_ds * ( d0.y()));
|
||||
// p3 (l1 end)
|
||||
set_J(0, 2 * l1.p1, df_dc * (-d0.y()) + df_ds * (-d0.x()));
|
||||
set_J(0, 2 * l1.p1 + 1, df_dc * ( d0.x()) + df_ds * (-d0.y()));
|
||||
break;
|
||||
}
|
||||
|
||||
case ConstraintType::FixedPoint:
|
||||
// Handled by excluding from Jacobian (fixed = true)
|
||||
break;
|
||||
|
||||
default:
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// DOF 统计
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
int ConstraintSolver::degrees_of_freedom() const {
|
||||
int dof = static_cast<int>(points_.size()) * 2;
|
||||
for(auto&p:points_) if(p.fixed) dof-=2;
|
||||
return dof - static_cast<int>(constraints_.size());
|
||||
for (auto& p : points_)
|
||||
if (p.fixed) dof -= 2;
|
||||
int eq_count = 0;
|
||||
for (auto& c : constraints_)
|
||||
eq_count += count_equations(c);
|
||||
return dof - eq_count;
|
||||
}
|
||||
|
||||
// ══════════════════════════════════════════════════════
|
||||
// Newton-Raphson 求解器
|
||||
// ══════════════════════════════════════════════════════
|
||||
|
||||
SolverResult ConstraintSolver::solve(int max_iter, double tol) {
|
||||
SolverResult r; int nv=points_.size(), nc=constraints_.size();
|
||||
if(nc==0){ r.converged=true; r.message="No constraints"; for(auto&p:points_) r.points.push_back(p.pos); return r; }
|
||||
SolverResult result;
|
||||
int nv = static_cast<int>(points_.size());
|
||||
int nc = static_cast<int>(constraints_.size());
|
||||
|
||||
std::vector<double> vars(nv*2);
|
||||
for(int i=0;i<nv;++i){ vars[i*2]=points_[i].pos.x(); vars[i*2+1]=points_[i].pos.y(); }
|
||||
// 没有约束：直接返回
|
||||
if (nc == 0) {
|
||||
result.converged = true;
|
||||
result.message = "No constraints";
|
||||
for (auto& p : points_) result.points.push_back(p.pos);
|
||||
return result;
|
||||
}
|
||||
|
||||
const double h=1e-6;
|
||||
for(int iter=0;iter<max_iter;++iter){
|
||||
std::vector<double> errs(nc);
|
||||
double max_err=0;
|
||||
for(int ci=0;ci<nc;++ci){ errs[ci]=eval_constraint(constraints_[ci]); max_err=std::max(max_err,std::abs(errs[ci])); }
|
||||
// ── 变量初始化 ──────────────────────────────────
|
||||
int n_vars = 2 * nv;
|
||||
Eigen::VectorXd vars(n_vars);
|
||||
for (int i = 0; i < nv; ++i) {
|
||||
vars[2 * i] = points_[i].pos.x();
|
||||
vars[2 * i + 1] = points_[i].pos.y();
|
||||
}
|
||||
|
||||
if(max_err<tol){ r.converged=true; r.iterations=iter; r.residual=max_err; r.message="Converged"; break; }
|
||||
// ── 自由变量映射 ────────────────────────────────
|
||||
// fixed=true 的点坐标从 Jacobian 中移除（硬约束）
|
||||
std::vector<int> free_idx; // 自由变量在 vars 中的索引
|
||||
std::vector<int> to_free(n_vars, -1); // 全局→自由变量索引
|
||||
for (int i = 0; i < nv; ++i) {
|
||||
if (!points_[i].fixed) {
|
||||
to_free[2 * i] = static_cast<int>(free_idx.size());
|
||||
free_idx.push_back(2 * i);
|
||||
to_free[2 * i + 1] = static_cast<int>(free_idx.size());
|
||||
free_idx.push_back(2 * i + 1);
|
||||
}
|
||||
}
|
||||
int n_free = static_cast<int>(free_idx.size());
|
||||
|
||||
double total=0;
|
||||
for(int vi=0;vi<nv;++vi){
|
||||
if(points_[vi].fixed) continue;
|
||||
for(int d=0;d<2;++d){
|
||||
double orig=vars[vi*2+d]; vars[vi*2+d]=orig+h;
|
||||
for(int i=0;i<nv;++i) points_[i].pos=Point2D(vars[i*2],vars[i*2+1]);
|
||||
if (n_free == 0) {
|
||||
// 所有点固定，直接返回
|
||||
result.converged = true;
|
||||
result.message = "All points fixed";
|
||||
for (auto& p : points_) result.points.push_back(p.pos);
|
||||
return result;
|
||||
}
|
||||
|
||||
double step=0; int cnt=0;
|
||||
for(int ci=0;ci<nc;++ci){
|
||||
double ep=eval_constraint(constraints_[ci]);
|
||||
double de=ep-errs[ci];
|
||||
if(std::abs(de)>1e-15){ step-=errs[ci]*h/de; cnt++; }
|
||||
}
|
||||
vars[vi*2+d]=orig;
|
||||
if(cnt>0){ step/=cnt; vars[vi*2+d]+=step; total+=std::abs(step); }
|
||||
// ── 方程计数 ────────────────────────────────────
|
||||
int n_eq = 0;
|
||||
std::vector<int> row_map(nc + 1, 0); // 约束 ci → 起始行
|
||||
for (int ci = 0; ci < nc; ++ci) {
|
||||
row_map[ci] = n_eq;
|
||||
n_eq += count_equations(constraints_[ci]);
|
||||
}
|
||||
row_map[nc] = n_eq; // 哨兵
|
||||
|
||||
// ════════════════════════════════════════════════
|
||||
// Newton-Raphson 迭代
|
||||
// ════════════════════════════════════════════════
|
||||
double prev_residual = 1e100;
|
||||
|
||||
for (int iter = 0; iter < max_iter; ++iter) {
|
||||
// ── a. 同步变量 → 几何实体 ──────────────────
|
||||
sync_from_vars(vars.data(), nv);
|
||||
|
||||
// ── b. 计算残差向量 ─────────────────────────
|
||||
Eigen::VectorXd f(n_eq);
|
||||
double max_err = 0.0;
|
||||
for (int ci = 0; ci < nc; ++ci) {
|
||||
const auto& c = constraints_[ci];
|
||||
int neq = count_equations(c);
|
||||
int r0 = row_map[ci];
|
||||
if (neq == 1) {
|
||||
f[r0] = eval_constraint(c);
|
||||
max_err = std::max(max_err, std::abs(f[r0]));
|
||||
} else if (c.type == ConstraintType::Coincident) {
|
||||
int p0 = c.elements[0], p1 = c.elements[1];
|
||||
f[r0] = points_[p1].pos.x() - points_[p0].pos.x();
|
||||
f[r0 + 1] = points_[p1].pos.y() - points_[p0].pos.y();
|
||||
max_err = std::max(max_err, std::max(std::abs(f[r0]), std::abs(f[r0 + 1])));
|
||||
} else if (c.type == ConstraintType::Fixed) {
|
||||
int pid = c.elements[0];
|
||||
// target = 当前位置（或从 extra params 解析）
|
||||
double tx = points_[pid].pos.x();
|
||||
double ty = points_[pid].pos.y();
|
||||
f[r0] = points_[pid].pos.x() - tx;
|
||||
f[r0 + 1] = points_[pid].pos.y() - ty;
|
||||
max_err = std::max(max_err, std::max(std::abs(f[r0]), std::abs(f[r0 + 1])));
|
||||
}
|
||||
}
|
||||
for(double&v:vars) v=std::clamp(v,-1e6,1e6);
|
||||
for(int i=0;i<nv;++i) points_[i].pos=Point2D(vars[i*2],vars[i*2+1]);
|
||||
r.iterations=iter+1; r.residual=max_err;
|
||||
if(total<tol){ r.converged=true; r.message="Step converged"; break; }
|
||||
|
||||
// ── c. 收敛判定 ─────────────────────────────
|
||||
if (max_err < tol) {
|
||||
result.converged = true;
|
||||
result.iterations = iter;
|
||||
result.residual = max_err;
|
||||
result.message = "Converged";
|
||||
for (int i = 0; i < nv; ++i)
|
||||
result.points.push_back(points_[i].pos);
|
||||
return result;
|
||||
}
|
||||
|
||||
// ── d. 步长收敛 ─────────────────────────────
|
||||
if (iter > 0 && std::abs(prev_residual - max_err) < tol * 0.1) {
|
||||
result.converged = true;
|
||||
result.iterations = iter;
|
||||
result.residual = max_err;
|
||||
result.message = "Stalled (step converged)";
|
||||
for (int i = 0; i < nv; ++i)
|
||||
result.points.push_back(points_[i].pos);
|
||||
return result;
|
||||
}
|
||||
prev_residual = max_err;
|
||||
|
||||
// ── e. 构建 Jacobian（只针对自由变量） ──────
|
||||
// J_free: n_eq × n_free
|
||||
Eigen::MatrixXd J(n_eq, n_free);
|
||||
J.setZero();
|
||||
|
||||
// 先构建完整 J_full，然后提取自由列
|
||||
// 使用临时缓冲区
|
||||
std::vector<double> J_full_buf(n_eq * n_vars, 0.0);
|
||||
for (int ci = 0; ci < nc; ++ci) {
|
||||
fill_jacobian_rows(constraints_[ci], ci, row_map,
|
||||
J_full_buf.data(), nv, n_vars);
|
||||
}
|
||||
|
||||
// 提取自由变量列
|
||||
for (int eq = 0; eq < n_eq; ++eq) {
|
||||
for (int fj = 0; fj < n_free; ++fj) {
|
||||
int col = free_idx[fj];
|
||||
J(eq, fj) = J_full_buf[eq * n_vars + col];
|
||||
}
|
||||
}
|
||||
|
||||
// ── f. 求解 J * dx = -f ────────────────────
|
||||
// 使用 SVD 伪逆，处理欠定/超定系统
|
||||
Eigen::VectorXd dx = J.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(-f);
|
||||
|
||||
// ── g. 阻尼：检查 step 是否改善 ─────────────
|
||||
double step_norm = dx.norm();
|
||||
|
||||
// 回溯线搜索
|
||||
double alpha = 1.0;
|
||||
Eigen::VectorXd vars_backup = vars;
|
||||
for (int ls = 0; ls < 8; ++ls) {
|
||||
// 应用步长
|
||||
for (int fj = 0; fj < n_free; ++fj) {
|
||||
vars[free_idx[fj]] += alpha * dx[fj];
|
||||
}
|
||||
// 限制范围
|
||||
for (int i = 0; i < n_vars; ++i) {
|
||||
vars[i] = std::clamp(vars[i], -1e8, 1e8);
|
||||
}
|
||||
sync_from_vars(vars.data(), nv);
|
||||
|
||||
// 计算新残差
|
||||
double new_err = 0.0;
|
||||
for (int ci = 0; ci < nc; ++ci) {
|
||||
const auto& c = constraints_[ci];
|
||||
int neq = count_equations(c);
|
||||
if (neq == 1) {
|
||||
new_err = std::max(new_err, std::abs(eval_constraint(c)));
|
||||
} else if (c.type == ConstraintType::Coincident) {
|
||||
int p0 = c.elements[0], p1 = c.elements[1];
|
||||
new_err = std::max(new_err, std::abs(points_[p1].pos.x() - points_[p0].pos.x()));
|
||||
new_err = std::max(new_err, std::abs(points_[p1].pos.y() - points_[p0].pos.y()));
|
||||
}
|
||||
}
|
||||
|
||||
if (new_err <= max_err * 1.01) {
|
||||
break; // 接受步长
|
||||
}
|
||||
// 回退并减半步长
|
||||
vars = vars_backup;
|
||||
alpha *= 0.5;
|
||||
if (alpha < 1e-8) break;
|
||||
}
|
||||
|
||||
// ── h. 更新状态 ─────────────────────────────
|
||||
sync_from_vars(vars.data(), nv);
|
||||
result.iterations = iter + 1;
|
||||
result.residual = max_err;
|
||||
}
|
||||
if(!r.converged) r.message="Max iterations";
|
||||
for(auto&p:points_) r.points.push_back(p.pos);
|
||||
return r;
|
||||
|
||||
// ════════════════════════════════════════════════
|
||||
// 未收敛
|
||||
// ════════════════════════════════════════════════
|
||||
if (!result.converged) {
|
||||
result.message = "Max iterations reached";
|
||||
// 即使未完全收敛，也返回当前最佳解
|
||||
double final_err = 0.0;
|
||||
for (int ci = 0; ci < nc; ++ci)
|
||||
final_err = std::max(final_err, std::abs(eval_constraint(constraints_[ci])));
|
||||
result.residual = final_err;
|
||||
}
|
||||
|
||||
for (int i = 0; i < nv; ++i)
|
||||
result.points.push_back(points_[i].pos);
|
||||
return result;
|
||||
}
|
||||
|
||||
} // namespace vde::sketch
|
||||
|
||||
@@ -12,3 +12,4 @@ add_subdirectory(spatial)
|
||||
add_subdirectory(collision)
|
||||
add_subdirectory(boolean)
|
||||
add_subdirectory(brep)
|
||||
add_subdirectory(sketch)
|
||||
|
||||
@@ -0,0 +1 @@
|
||||
add_vde_test(test_constraint_solver)
|
||||
Reference in New Issue
Block a user