diff --git a/application/test_circle_gwn.py b/application/test_circle_gwn.py new file mode 100644 index 0000000..3059d94 --- /dev/null +++ b/application/test_circle_gwn.py @@ -0,0 +1,198 @@ +import math +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + + +def arc_gwn_batch(a, b, bulge, P): + """ + 向量化批处理接口：对 (N, 2) 的点阵 P 同时计算单段圆弧的广义绕数 W。 + """ + a = np.asarray(a, dtype=float) + b = np.asarray(b, dtype=float) + P = np.asarray(P, dtype=float) + + half_ab = 0.5 * (a + b) + bisec_vec = 0.5 * np.array([a[1] - b[1], b[0] - a[0]]) + center = half_ab + bisec_vec / bulge + + ca = a - center + cb = b - center + cp = P - center + + D = np.linalg.norm(cp, axis=1) + r = np.linalg.norm(ca) + theta = 2.0 * math.atan(bulge) + + W = np.empty(len(P)) + + mask_center = D < 1e-12 + mask_circle = np.abs(D - r) < 1e-12 + mask_normal = ~(mask_circle) + + # W[mask_center] = theta / (2.0 * math.pi) + W[mask_circle] = theta / (4.0 * math.pi) + + cp_n = cp[mask_normal] + D_n = D[mask_normal] + + sigma = np.where(r > D_n, 1.0, -1.0) + Rd = r * D_n + r_p_D = r + D_n + abs_r_D = np.abs(r - D_n) + + dot_a = ca.dot(cp_n.T) + cross1 = cp_n[:, 0] * ca[1] - ca[0] * cp_n[:, 1] + Ga = np.arctan2(r_p_D * cross1, abs_r_D * (Rd + dot_a)) + + dot_b = cb.dot(cp_n.T) + cross2 = cp_n[:, 0] * cb[1] - cb[0] * cp_n[:, 1] + Gb = np.arctan2(r_p_D * cross2, abs_r_D * (Rd + dot_b)) + + delta = Gb - Ga + fix = delta * theta < 0.0 + delta[fix] += math.copysign(math.pi, theta) + + two_pi_W = 0.5 * theta + sigma * delta + W[mask_normal] = two_pi_W / (2.0 * math.pi) + + return W + + +def plot_multi_arc_gwn_field(arcs, margin=0.5, nx=400, ny=400): + """ + 绘制多段圆弧的广义绕数之和。 + + 参数 + ---- + arcs : list of tuple + 每段圆弧表示为 (a, b, bulge)，其中 a, b 为长度为 2 的序列。 + margin : float + AABB 外扩边界。 + nx, ny : int + 采样分辨率。 + """ + arcs = [(np.asarray(a, dtype=float), np.asarray(b, dtype=float), float(bulge)) + for a, b, bulge in arcs] + + # ---- 1. 计算所有圆弧的几何信息，确定全局 AABB ---- + all_arc_x = [] + all_arc_y = [] + all_centers = [] + all_rs = [] + + for a, b, bulge in arcs: + half_ab = 0.5 * (a + b) + bisec_vec = 0.5 * np.array([a[1] - b[1], b[0] - a[0]]) + center = half_ab + bisec_vec / bulge + r = np.linalg.norm(a - center) + theta = 2.0 * math.atan(bulge) + + start_angle = math.atan2(a[1] - center[1], a[0] - center[0]) + end_angle = start_angle + theta + t_arc = np.linspace(start_angle, end_angle, 200) + arc_x = center[0] + r * np.cos(t_arc) + arc_y = center[1] + r * np.sin(t_arc) + + all_arc_x.append(arc_x) + all_arc_y.append(arc_y) + all_centers.append(center) + all_rs.append(r) + + # 全局 AABB + x_min = min(ax.min() for ax in all_arc_x) - margin + x_max = max(ax.max() for ax in all_arc_x) + margin + y_min = min(ay.min() for ay in all_arc_y) - margin + y_max = max(ay.max() for ay in all_arc_y) + margin + + # ---- 2. 均匀采样并累加所有圆弧的 GWN ---- + xs = np.linspace(x_min, x_max, nx) + ys = np.linspace(y_min, y_max, ny) + X, Y = np.meshgrid(xs, ys) + + points = np.stack([X.ravel(), Y.ravel()], axis=1) + W = np.zeros(len(points)) + + for a, b, bulge in arcs: + W += arc_gwn_batch(a, b, bulge, points) + + W = W.reshape(X.shape) + + # ---- 3. 三分区独立归一化 ---- + eps = np.finfo(float).eps * 1e6 + bound_low = -eps + bound_high = eps + + # 负区 + W_neg = W.copy() + W_neg[W >= bound_low] = np.nan + W_min = np.nanmin(W_neg) + if not np.isnan(W_min): + W_neg = (bound_low - W_neg) / (bound_low - W_min) + + # 零区 + W_zero = W.copy() + W_zero[(W < bound_low) | (W > bound_high)] = np.nan + W_zero = (W_zero - bound_low) / (2.0 * eps) + + # 正区 + W_pos = W.copy() + W_pos[W <= bound_high] = np.nan + W_max = np.nanmax(W_pos) + if not np.isnan(W_max): + W_pos = (W_pos - bound_high) / (W_max - bound_high) + + # ---- 4. 绘图 ---- + fig, ax = plt.subplots(figsize=(10, 8)) + + im_neg = ax.pcolormesh(X, Y, W_neg, cmap='Blues', vmin=0, vmax=1, shading='auto') + im_zero = ax.pcolormesh(X, Y, W_zero, cmap='Greys', vmin=0, vmax=1, shading='auto') + im_pos = ax.pcolormesh(X, Y, W_pos, cmap='Reds', vmin=0, vmax=1, shading='auto') + + # 叠加所有圆弧 + for idx, (a, b, bulge) in enumerate(arcs): + ax.plot(all_arc_x[idx], all_arc_y[idx], 'k-', linewidth=2.5, zorder=4) + ax.plot(a[0], a[1], 'go', markersize=7, zorder=5) + ax.plot(b[0], b[1], 'mo', markersize=7, zorder=5) + ax.plot(all_centers[idx][0], all_centers[idx][1], 'k+', markersize=8, mew=1.5, zorder=5) + + ax.set_aspect('equal') + ax.set_xlim(x_min, x_max) + ax.set_ylim(y_min, y_max) + ax.set_title(f'Multi-Arc GWN Sum ({len(arcs)} arcs)') + + # ---- 5. 三个独立 colorbar ---- + divider = make_axes_locatable(ax) + + cax_neg = divider.append_axes("left", size="5%", pad=0.4) + cbar_neg = fig.colorbar(im_neg, cax=cax_neg, orientation='vertical') + cbar_neg.ax.yaxis.set_ticks_position('left') + cbar_neg.ax.yaxis.set_label_position('left') + cbar_neg.set_label('Negative $W$ (normalized)', rotation=90, labelpad=10) + cbar_neg.set_ticks([0, 1]) + cbar_neg.set_ticklabels([f'{bound_low:.2e}', f'{W_min:.3f}' if not np.isnan(W_min) else 'N/A']) + + cax_pos = divider.append_axes("right", size="5%", pad=0.4) + cbar_pos = fig.colorbar(im_pos, cax=cax_pos, orientation='vertical') + cbar_pos.set_label('Positive $W$ (normalized)', rotation=270, labelpad=15) + cbar_pos.set_ticks([0, 1]) + cbar_pos.set_ticklabels([f'{bound_high:.2e}', f'{W_max:.3f}' if not np.isnan(W_max) else 'N/A']) + + cax_zero = divider.append_axes("bottom", size="5%", pad=0.6) + cbar_zero = fig.colorbar(im_zero, cax=cax_zero, orientation='horizontal') + cbar_zero.set_label(f'Near-zero $W$ (normalized) [{bound_low:.2e}, {bound_high:.2e}]', labelpad=5) + cbar_zero.set_ticks([0, 1]) + cbar_zero.set_ticklabels([f'{bound_low:.2e}', f'{bound_high:.2e}']) + + plt.tight_layout() + plt.show() + + +# ==================== 主程序 ==================== +if __name__ == '__main__': + # 示例 1：两段圆弧拼成近似整圆 + arcs = [ + ([2.0, 0.0], [0.0, 0.0], float("inf")), # 上半圆 + ([0.0, 0.0], [2.0, 0.0], float("inf")), # 下半圆（注意 bulge>0 表示逆时针） + ] + plot_multi_arc_gwn_field(arcs, margin=0.5, nx=1024, ny=1024) \ No newline at end of file diff --git a/application/test_extrude.cpp b/application/test_extrude.cpp index e4816f3..05a87e9 100644 --- a/application/test_extrude.cpp +++ b/application/test_extrude.cpp @@ -1,85 +1,280 @@ #include -#include +#include +#include #include #include "io_primitive.h" #include "math/math_defs.hpp" #include "mimalloc.h" #include "primitive_descriptor.h" - #include -int main() +void export_obj(const char* filename, const polymesh_t& mesh) { + std::ofstream f(filename); + if (!f) { std::cerr << "Cannot open " << filename << std::endl; return; } + for (uint32_t i = 0; i < mesh.num_vertices; ++i) { + auto& v = mesh.vertices[i]; + f << "v " << v.x << " " << v.y << " " << v.z << "\n"; + } + uint32_t offset = 0; + for (uint32_t fi = 0; fi < mesh.num_faces; ++fi) { + uint32_t count = mesh.vertex_counts[fi]; + f << "f"; + for (uint32_t j = 0; j < count; ++j) + f << " " << (mesh.faces[offset + j] + 1); + f << "\n"; + offset += count; + } + std::cout << " -> " << filename << ": " << mesh.num_vertices << " verts, " << mesh.num_faces << " faces" << std::endl; +} + +void run_test_case(const char* name, + polyline_pattern_descriptor_t& pattern_desc, + polyline_axis_descriptor_t& axis_desc) { - helixline_axis_descriptor_t axis_desc{}; - axis_desc.axis_start = {0, 0, 0}; - axis_desc.axis_end = {0, 0, 13}; - axis_desc.is_righthanded = true; - axis_desc.radius = 1.; - axis_desc.advance_per_round = 1; - axis_desc.start_direction = {1, 0, 0}; - - // polyline_axis_descriptor_t axis_desc{}; - // std::vector points_{ - // vector3d{0, 0, 0}, - // vector3d{0, 1, 1} - // }; - // axis_desc.point_number = 2; - // axis_desc.points = points_.data(); - // std::vector bulges_{std::tan(pi / 8)}; - // axis_desc.bulge_number = 1; - // axis_desc.bulges = bulges_.data(); - // axis_desc.reference_normal = {1, 0, 0}; - // axis_desc.is_closed = false; - - polyline_pattern_descriptor_t pattern_desc{}; - pattern_desc.anchor = {0, 0}; - std::vector points{ - vector2d{-.1, 0}, - vector2d{.1, 0}, - vector2d{-.1, 0} - }; - pattern_desc.point_number = 3; - pattern_desc.points = points.data(); - std::vector bulges{1., 1.}; - pattern_desc.bulge_number = 2; - pattern_desc.bulges = bulges.data(); - - auto primitive_data_center = create_primitive_data_center(); - auto pattern = create_pattern(primitive_data_center, &pattern_desc, PATTERN_TYPE_POLYLINE); - auto axis = create_extrude_axis(primitive_data_center, &axis_desc, AXIS_TYPE_EXTRUDE_HELIXLINE); - auto extrude = create_param_primitive(&axis, &pattern, primitive_data_center, PRIMITIVE_TYPE_EXTRUDE); - std::cout << "primitive created..." << std::endl; - - auto runtime_blobtree = create_blobtree(); - auto node_iter1 = blobtree_add_primitive_node(runtime_blobtree, extrude); - auto baked_blobtree = bake_blobtree(runtime_blobtree); - destroy_blobtree(runtime_blobtree); - - std::cout << "blobtree created..." << std::endl; - - s_settings settings{}; - settings.resolution = 32; - settings.scene_aabb_margin = 1e-5; - settings.restricted_primitive_bounding_test = true; - settings.integrand_density = 32; - auto solver = create_solver(baked_blobtree, settings); + auto dc = create_primitive_data_center(); + auto pat = create_pattern(dc, &pattern_desc, PATTERN_TYPE_POLYLINE); + auto ax = create_extrude_axis(dc, &axis_desc, AXIS_TYPE_EXTRUDE_POLYLINE); + auto ext = create_param_primitive(&ax, &pat, dc, PRIMITIVE_TYPE_EXTRUDE); + auto rt = create_blobtree(); + blobtree_add_primitive_node(rt, ext); + auto baked = bake_blobtree(rt); + destroy_blobtree(rt); + s_settings s{}; + s.resolution = 32; + s.scene_aabb_margin = 1e-5; + s.restricted_primitive_bounding_test = true; + s.integrand_density = 32; + auto solver = create_solver(baked, s); auto result = generate_polymesh(solver); - // auto integrator = generate_integral_scheme(solver); - print_statistics(solver); + std::cout << name << ": " << (result.success ? "OK" : "FAIL") + << " verts=" << result.mesh.num_vertices + << " faces=" << result.mesh.num_faces << std::endl; - // integrand_handle_t handle{1., 0, 0, 0}; - // std::cout << "surface area: " << calculate_surface_integrand(integrator, handle) << std::endl; - // std::cout << "volume: " << calculate_volume_integrand(integrator, handle) << std::endl; - // std::cout << "baseline (by polymesh): \n"; - // auto baseline_integral = polymesh_area_and_volume(solver); - // std::cout << "surface area: " << baseline_integral.area << std::endl; - // std::cout << "volume: " << baseline_integral.volume << std::endl; + char filename[256]; + snprintf(filename, sizeof(filename), "test_extrude_%s.obj", name); + export_obj(filename, result.mesh); destroy_solver(solver); - destroy_baked_blobtree(baked_blobtree); - destroy_primitive(extrude); - destroy_primitive_data_center(primitive_data_center); + destroy_baked_blobtree(baked); + destroy_primitive(ext); + destroy_primitive_data_center(dc); +} + +int main() +{ + const double b_quarter = std::tan(pi / 8); // tan(theta/4), theta=pi/2 -> quarter circle + + // --- Curved axis: quarter-circle arc from (0,0,0) to (0,1,1) --- + polyline_axis_descriptor_t axis_curved{}; + std::vector axis_pts_curved{ vector3d{0, 0, 0}, vector3d{0, 1, 1} }; + axis_curved.point_number = 2; + axis_curved.points = axis_pts_curved.data(); + std::vector axis_blg_curved{ b_quarter }; + axis_curved.bulge_number = 1; + axis_curved.bulges = axis_blg_curved.data(); + axis_curved.reference_normal = {1, 0, 0}; + axis_curved.is_closed = false; + + // --- Straight axis: from (0,0,0) to (0,0,1) --- + polyline_axis_descriptor_t axis_straight{}; + std::vector axis_pts_straight{ vector3d{0, 0, 0}, vector3d{0, 0, 1} }; + axis_straight.point_number = 2; + axis_straight.points = axis_pts_straight.data(); + std::vector axis_blg_straight{ 0. }; + axis_straight.bulge_number = 1; + axis_straight.bulges = axis_blg_straight.data(); + axis_straight.reference_normal = {1, 0, 0}; + axis_straight.is_closed = false; + + // ================================================================== + // Test 1: Original semicircle (reference, all arcs, legacy bulge=1) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ vector2d{-.1, 0}, vector2d{.1, 0}, vector2d{-.1, 0} }; + pat.point_number = 3; + pat.points = pts.data(); + std::vector blg{ 1., 1. }; + pat.bulge_number = 2; + pat.bulges = blg.data(); + run_test_case("01_original_semicircle", pat, axis_curved); + } + + // ================================================================== + // Test 2: Rectangle (all lines, bulge=0) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.1, -.1}, vector2d{.1, -.1}, + vector2d{.1, .1}, vector2d{-.1, .1}, + vector2d{-.1, -.1} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ 0., 0., 0., 0. }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("02_rectangle", pat, axis_curved); + } + + // ================================================================== + // Test 3: Curved diamond (4 quarter-circle arcs, total theta=2pi) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.12, 0}, vector2d{0, .12}, + vector2d{.12, 0}, vector2d{0, -.12}, + vector2d{-.12, 0} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ b_quarter, b_quarter, b_quarter, b_quarter }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("03_curved_diamond", pat, axis_curved); + } + + // ================================================================== + // Test 4: Rectangle with one rounded side (3 lines + 1 arc) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.12, -.08}, vector2d{.12, -.08}, + vector2d{.12, .08}, vector2d{-.12, .08}, + vector2d{-.12, -.08} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ 0., 0., 0., b_quarter }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("04_rect_one_arc", pat, axis_curved); + } + + // ================================================================== + // Test 5: Alternating line-arc hexagon (3 lines + 3 quarter arcs) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.10, -.05}, vector2d{0., -.10}, + vector2d{.10, -.05}, vector2d{.10, .05}, + vector2d{0., .10}, vector2d{-.10, .05}, + vector2d{-.10, -.05} + }; + pat.point_number = 7; + pat.points = pts.data(); + std::vector blg{ 0., b_quarter, 0., b_quarter, 0., b_quarter }; + pat.bulge_number = 6; + pat.bulges = blg.data(); + run_test_case("05_alternating_hex", pat, axis_curved); + } + + // ================================================================== + // Test 6: Rectangle on straight axis (baseline, all lines) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.1, -.1}, vector2d{.1, -.1}, + vector2d{.1, .1}, vector2d{-.1, .1}, + vector2d{-.1, -.1} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ 0., 0., 0., 0. }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("06_rect_straight_axis", pat, axis_straight); + } + + // ================================================================== + // Test 7: Curved diamond on straight axis + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.12, 0}, vector2d{0, .12}, + vector2d{.12, 0}, vector2d{0, -.12}, + vector2d{-.12, 0} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ b_quarter, b_quarter, b_quarter, b_quarter }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("07_curved_diamond_straight", pat, axis_straight); + } + + // ================================================================== + // Test 8: Pill/Stadium shape (2 lines + 2 quarter arcs, straight axis) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.06, -.12}, vector2d{-.06, .12}, + vector2d{ .06, .12}, vector2d{ .06, -.12}, + vector2d{-.06, -.12} + }; + pat.point_number = 5; + pat.points = pts.data(); + std::vector blg{ 0., b_quarter, 0., b_quarter }; + pat.bulge_number = 4; + pat.bulges = blg.data(); + run_test_case("08_pill_shape", pat, axis_straight); + } + + // ================================================================== + // Test 9: Mixed octagon (4 lines + 4 arcs, straight axis) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.12, -.08}, vector2d{-.06, -.12}, + vector2d{ .06, -.12}, vector2d{ .12, -.08}, + vector2d{ .12, .08}, vector2d{ .06, .12}, + vector2d{-.06, .12}, vector2d{-.12, .08}, + vector2d{-.12, -.08} + }; + pat.point_number = 9; + pat.points = pts.data(); + std::vector blg{ b_quarter, 0., b_quarter, 0., b_quarter, 0., b_quarter, 0. }; + pat.bulge_number = 8; + pat.bulges = blg.data(); + run_test_case("09_mixed_octagon", pat, axis_straight); + } + + // ================================================================== + // Test 10: Rounded triangle (2 lines + 1 arc, straight axis) + // ================================================================== + { + polyline_pattern_descriptor_t pat{}; + pat.anchor = {0, 0}; + std::vector pts{ + vector2d{-.10, -.10}, vector2d{.10, -.10}, + vector2d{-.10, .10}, vector2d{-.10, -.10} + }; + pat.point_number = 4; + pat.points = pts.data(); + std::vector blg{ 0., b_quarter, 0. }; + pat.bulge_number = 3; + pat.bulges = blg.data(); + run_test_case("10_rounded_triangle", pat, axis_straight); + } + std::cout << "\nAll tests completed." << std::endl; return 0; } \ No newline at end of file diff --git a/application/test_occ.py b/application/test_occ.py new file mode 100644 index 0000000..1263978 --- /dev/null +++ b/application/test_occ.py @@ -0,0 +1,324 @@ +import numpy as np +import time +from math import cos, sin, pi, sqrt +from typing import List, Optional + +# --- OCC Core imports --- +from OCC.Core.gp import ( + gp_Pnt, gp_Dir, gp_Ax2, gp_Ax3, gp_Vec, gp_Trsf, gp_Circ +) +from OCC.Core.BRepPrimAPI import ( + BRepPrimAPI_MakeBox, + BRepPrimAPI_MakeCylinder, + BRepPrimAPI_MakeSphere +) +from OCC.Core.BRepBuilderAPI import ( + BRepBuilderAPI_MakeEdge, BRepBuilderAPI_MakeWire, + BRepBuilderAPI_MakeFace, BRepBuilderAPI_Transform +) +from OCC.Core.BRepOffsetAPI import BRepOffsetAPI_MakePipe +from OCC.Core.BRepAlgoAPI import ( + BRepAlgoAPI_Fuse, BRepAlgoAPI_Cut, BRepAlgoAPI_Common +) +from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh +from OCC.Core.BRep import BRep_Tool +from OCC.Core.TopExp import TopExp_Explorer +from OCC.Core.TopAbs import TopAbs_FACE, TopAbs_REVERSED +from OCC.Core.TopoDS import topods_Face +from OCC.Core.TopLoc import TopLoc_Location +from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepBndLib import brepbndlib_Add +from OCC.Core.TColgp import TColgp_Array1OfPnt +from OCC.Core.GeomAPI import GeomAPI_PointsToBSpline +from OCC.Core.BRepCheck import BRepCheck_Analyzer + + +# ============================================================ +# 1. 基础几何体构造 +# ============================================================ +def make_box(x: float, y: float, z: float, + dx: float, dy: float, dz: float): + axes = gp_Ax2(gp_Pnt(x, y, z), gp_Dir(0, 0, 1)) + return BRepPrimAPI_MakeBox(axes, dx, dy, dz).Shape() + + +def make_cylinder(cx: float, cy: float, cz: float, + radius: float, height: float, + axis_dir=None): + if axis_dir is None: + axis_dir = gp_Dir(0, 0, 1) + axes = gp_Ax2(gp_Pnt(cx, cy, cz), axis_dir) + return BRepPrimAPI_MakeCylinder(axes, radius, height).Shape() + + +def make_sphere(cx: float, cy: float, cz: float, radius: float): + center = gp_Pnt(cx, cy, cz) + return BRepPrimAPI_MakeSphere(center, radius).Shape() + + +# ============================================================ +# 2. 螺旋线扫掠 +# ============================================================ +def make_bspline_helix(radius: float, pitch: float, + t0: float, t1: float, n_points: int = 128): + points = TColgp_Array1OfPnt(1, n_points) + for i in range(n_points): + t = t0 + (t1 - t0) * i / (n_points - 1) + points.SetValue(i + 1, gp_Pnt( + radius * cos(t), + radius * sin(t), + pitch * t / (2.0 * pi) + )) + approx = GeomAPI_PointsToBSpline(points) + return approx.Curve() + + +def make_helix_sweep(helix_r: float, pitch: float, turns: float, + section_r: float, n_points: int = 128): + t0, t1 = 0.0, 2.0 * pi * turns + helix_curve = make_bspline_helix(helix_r, pitch, t0, t1, n_points) + + helix_edge = BRepBuilderAPI_MakeEdge(helix_curve).Edge() + helix_wire = BRepBuilderAPI_MakeWire(helix_edge).Wire() + + sec_axes = gp_Ax2(gp_Pnt(helix_r, 0.0, 0.0), gp_Dir(0, 0, 1)) + sec_circ = gp_Circ(sec_axes, section_r) + sec_edge = BRepBuilderAPI_MakeEdge(sec_circ).Edge() + sec_wire = BRepBuilderAPI_MakeWire(sec_edge).Wire() + sec_face = BRepBuilderAPI_MakeFace(sec_wire).Face() + + pipe = BRepOffsetAPI_MakePipe(helix_wire, sec_face) + pipe.Build() + if not pipe.IsDone(): + raise RuntimeError("Helix sweep failed") + return pipe.Shape() + + +# ============================================================ +# 3. 布尔运算 +# ============================================================ +def boolean_op(a, b, op: str = "fuse"): + if op == "fuse": + api = BRepAlgoAPI_Fuse(a, b) + elif op == "cut": + api = BRepAlgoAPI_Cut(a, b) + elif op == "common": + api = BRepAlgoAPI_Common(a, b) + else: + raise ValueError(f"Unknown op: {op}") + + if not api.IsDone(): + raise RuntimeError(f"Boolean {op} failed") + + result = api.Shape() + analyzer = BRepCheck_Analyzer(result) + if not analyzer.IsValid(): + print("[Warning] Boolean result has invalid topology") + return result + + +def cascade_boolean(shapes: List, op: str = "fuse") -> Optional: + if not shapes: + return None + if len(shapes) == 1: + return shapes[0] + + current = shapes[:] + while len(current) > 1: + nxt = [] + for i in range(0, len(current) - 1, 2): + nxt.append(boolean_op(current[i], current[i + 1], op)) + if len(current) % 2 == 1: + nxt.append(current[-1]) + current = nxt + return current[0] + + +# ============================================================ +# 4. 网格提取 +# ============================================================ +def extract_triangulation(shape, merge_tol: float = 1e-6): + vertices = [] + faces = [] + vmap = {} + + def quantize(v: float) -> int: + return int(round(v / merge_tol)) + + explorer = TopExp_Explorer(shape, TopAbs_FACE) + while explorer.More(): + face = topods_Face(explorer.Current()) + + loc = TopLoc_Location() + tri = BRep_Tool.Triangulation(face, loc) + if tri is None: + explorer.Next() + continue + + trsf = loc.Transformation() + n_nodes = tri.NbNodes() + n_tri = tri.NbTriangles() + + local2global = [0] * (n_nodes + 1) + + for i in range(1, n_nodes + 1): + p = tri.Node(i).Transformed(trsf) + key = (quantize(p.X()), quantize(p.Y()), quantize(p.Z())) + if key not in vmap: + idx = len(vertices) + vmap[key] = idx + vertices.append((p.X(), p.Y(), p.Z())) + local2global[i] = idx + else: + local2global[i] = vmap[key] + + for i in range(1, n_tri + 1): + t = tri.Triangle(i) + n1, n2, n3 = t.Get() + if face.Orientation() == TopAbs_REVERSED: + n1, n2 = n2, n1 + faces.append(( + local2global[n1], + local2global[n2], + local2global[n3] + )) + + explorer.Next() + + return ( + np.array(vertices, dtype=np.float64), + np.array(faces, dtype=np.int32) + ) + + +# ============================================================ +# 5. 自适应目标面数网格化 +# ============================================================ +def bounding_diagonal(shape) -> float: + bbox = Bnd_Box() + brepbndlib_Add(shape, bbox) + x1, y1, z1, x2, y2, z2 = bbox.Get() + return sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2) + + +def mesh_to_target(shape, target_faces: int, + ang_deflection: float = 0.5, + tol_percent: float = 0.05, + max_iter: int = 25): + diag = bounding_diagonal(shape) + lo = 1e-4 * diag + hi = diag + + best = None + best_diff = float('inf') + + for _ in range(max_iter): + mid = (lo + hi) * 0.5 + + BRepMesh_IncrementalMesh(shape, mid, + False, + ang_deflection, + True) + + verts, faces = extract_triangulation(shape) + diff = faces.shape[0] - target_faces + abs_diff = abs(diff) + + if abs_diff < best_diff: + best_diff = abs_diff + best = (verts.copy(), faces.copy()) + + if abs_diff <= target_faces * tol_percent: + break + + if diff > 0: + lo = mid + else: + hi = mid + + return best + + +# ============================================================ +# 6. 导出为 Wavefront OBJ +# ============================================================ +def write_obj(path: str, vertices: np.ndarray, faces: np.ndarray): + """ + 导出 Wavefront OBJ 格式。 + 面索引使用 1-based（OBJ 标准）。 + """ + with open(path, 'w') as f: + f.write("# Generated by pythonOCC\n") + for v in vertices: + f.write(f"v {v[0]:.17g} {v[1]:.17g} {v[2]:.17g}\n") + for tri in faces: + # OBJ 为 1-based 索引 + f.write(f"f {tri[0]+1} {tri[1]+1} {tri[2]+1}\n") + + +# ============================================================ +# 7. 主流程（含计时，仅输出 OBJ） +# ============================================================ +if __name__ == "__main__": + t_total_0 = time.perf_counter() + + # -------------------------------------------------------- + # 7.1 几何体构造 + # -------------------------------------------------------- + t_geom_0 = time.perf_counter() + + # base = make_box(0, 0, 0, 100, 80, 20) + # boss = make_cylinder(50, 40, 20, 15, 30) + # hole = make_cylinder(50, 40, 0, 8, 60) + # sphere = make_sphere(50, 40, 55, 12) + + # spring = make_helix_sweep(20, 10, 3, 3) + # trsf = gp_Trsf() + # trsf.SetTranslation(gp_Vec(20, 20, 20)) + # spring = BRepBuilderAPI_Transform(spring, trsf).Shape() + + cylinder = make_cylinder(0, 0, 0, 1, 1) + sphere = make_sphere(1.9, 0, 0, 1) + + t_geom_1 = time.perf_counter() + print(f"[Timing] Geometry construction: {t_geom_1 - t_geom_0:.4f} s") + + # -------------------------------------------------------- + # 7.2 布尔运算 + # -------------------------------------------------------- + t_bool_0 = time.perf_counter() + + # body = cascade_boolean([base, boss, sphere, spring], "fuse") + part = boolean_op(cylinder, sphere, "common") + + t_bool_1 = time.perf_counter() + print(f"[Timing] Boolean operations: {t_bool_1 - t_bool_0:.4f} s") + + # -------------------------------------------------------- + # 7.3 网格化 + # -------------------------------------------------------- + t_mesh_0 = time.perf_counter() + + verts, faces = mesh_to_target(part, target_faces=10000, + ang_deflection=0.5, + tol_percent=0.03) + + t_mesh_1 = time.perf_counter() + print(f"[Timing] Meshing (target 50k faces): {t_mesh_1 - t_mesh_0:.4f} s") + print(f"[Result] {verts.shape[0]} vertices, {faces.shape[0]} faces") + + # -------------------------------------------------------- + # 7.4 导出 OBJ + # -------------------------------------------------------- + t_export_0 = time.perf_counter() + + write_obj("common.obj", verts, faces) + + t_export_1 = time.perf_counter() + print(f"[Timing] OBJ export: {t_export_1 - t_export_0:.4f} s") + + # -------------------------------------------------------- + # 总计时 + # -------------------------------------------------------- + t_total_1 = time.perf_counter() + print(f"[Timing] Total elapsed: {t_total_1 - t_total_0:.4f} s") \ No newline at end of file diff --git a/application/test_occ_draw.py b/application/test_occ_draw.py new file mode 100644 index 0000000..7d8dc3d --- /dev/null +++ b/application/test_occ_draw.py @@ -0,0 +1,117 @@ +import vtk +import datetime + +if __name__ == "__main__": + # Create renderer, render window, and interactor + renderer = vtk.vtkRenderer() + render_window = vtk.vtkRenderWindow() + render_window.AddRenderer(renderer) + render_window.SetSize(1200, 800) + render_window_interactor = vtk.vtkRenderWindowInteractor() + render_window_interactor.SetRenderWindow(render_window) + + # 添加坐标轴 + axes = vtk.vtkAxesActor() + axes_widget = vtk.vtkOrientationMarkerWidget() + axes_widget.SetOutlineColor(0.93, 0.57, 0.13) + axes_widget.SetOrientationMarker(axes) + axes_widget.SetInteractor(render_window_interactor) + axes_widget.SetViewport(0.0, 0.0, 0.2, 0.2) + axes_widget.EnabledOn() + axes_widget.InteractiveOn() + + # Initialize OBJ importer + obj_importer = vtk.vtkOBJImporter() + obj_importer.SetFileName("common.obj") # Replace with your OBJ file path + # obj_importer.SetFileNameMTL("halfpatch_final.mtl") # Replace with your MTL file path + obj_importer.SetTexturePath("./") # Directory containing texture files + # Set render window for the importer + obj_importer.SetRenderWindow(render_window) + obj_importer.Update() + + reader = vtk.vtkOBJReader() + reader.SetFileName("halfpatch_final.obj") # Replace with your OBJ file path + reader.Update() + polydata = reader.GetOutput() + triangle_filter = vtk.vtkTriangleFilter() + triangle_filter.SetInputData(polydata) + triangle_filter.Update() + polydata_tri = triangle_filter.GetOutput() + + mass = vtk.vtkMassProperties() + mass.SetInputData(polydata_tri) + print("Surface Area:", mass.GetSurfaceArea()) + print("Volume:", mass.GetVolume()) + + # 获取导入模型的所有actor + actors = renderer.GetActors() + actors.InitTraversal() + + # 创建线框actor列表 + wireframe_actors = [] + + # 为每个模型actor创建对应的线框actor + for i in range(actors.GetNumberOfItems()): + actor = actors.GetNextActor() + + # 创建线框表示 + wireframe_mapper = vtk.vtkPolyDataMapper() + wireframe_mapper.SetInputData(actor.GetMapper().GetInput()) + + # 创建线框actor + wireframe_actor = vtk.vtkActor() + wireframe_actor.SetMapper(wireframe_mapper) + + # 设置线框属性 + wireframe_actor.GetProperty().SetRepresentationToWireframe() # 设置为线框模式 + wireframe_actor.GetProperty().SetColor(1.0, 1.0, 1.0) # 设置线框颜色为白色 + wireframe_actor.GetProperty().SetLineWidth(2.0) # 设置线框线条粗细，可以调整这个值 + + # 使线框显示在模型上方（通过调整渲染顺序） + wireframe_actor.GetProperty().SetOpacity(0.5) # 设置不透明度 + wireframe_actor.GetProperty().LightingOff() # 关闭光照，使线框更清晰 + + # 设置线框actor的位置和方向与原始actor一致 + wireframe_actor.SetPosition(actor.GetPosition()) + wireframe_actor.SetOrientation(actor.GetOrientation()) + wireframe_actor.SetScale(actor.GetScale()) + + # 添加到渲染器和列表 + renderer.AddActor(wireframe_actor) + wireframe_actors.append(wireframe_actor) + + # Set background color + renderer.SetBackground(0.1, 0.1, 0.1) # Dark gray background + + # 设置默认交互器样式 + style = vtk.vtkInteractorStyleTrackballCamera() + render_window_interactor.SetInteractorStyle(style) + style.SetMotionFactor(5.0) # 默认是10.0，较小的值会减慢移动速度 + + # 定义按键回调函数，处理F10截图 + def on_key_press(obj, event): + key = obj.GetKeySym() + if key == 'F10': + # 创建窗口到图像的过滤器，捕获整个窗口内容 + window_to_image = vtk.vtkWindowToImageFilter() + window_to_image.SetInput(obj.GetRenderWindow()) + window_to_image.Update() + + # 生成带时间戳的文件名，避免多次截图覆盖旧文件 + timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") + filename = f"screenshot_{timestamp}.png" + + # 将图像写入PNG文件 + writer = vtk.vtkPNGWriter() + writer.SetFileName(filename) + writer.SetInputConnection(window_to_image.GetOutputPort()) + writer.Write() + + print(f"截图已成功保存到当前目录：{filename}") + + # 绑定按键事件观察者，监听键盘输入 + render_window_interactor.AddObserver('KeyPressEvent', on_key_press) + + # Render and start interaction + render_window.Render() + render_window_interactor.Start() \ No newline at end of file