ZCWang
/
ImplicitSurfaceNetwork
1
0
Fork 0
Code Issues 1 Pull Requests Projects Releases Wiki Activity
Browse Source

update tests

V2-integral
Zhicheng Wang 5 days ago
parent
7d79d2bc77
commit
57b015fdc2
4 changed files with 902 additions and 68 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 198
      application/test_circle_gwn.py
  2. 331
      application/test_extrude.cpp
  3. 324
      application/test_occ.py
  4. 117
      application/test_occ_draw.py

198
application/test_circle_gwn.py
View File

@ -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)

331
application/test_extrude.cpp
View File

@ -1,85 +1,280 @@
#include <iostream>
#include <array>
#include <fstream>
#include <cstdio>
#include <vector>
#include "io_primitive.h"
#include "math/math_defs.hpp"
#include "mimalloc.h"
#include "primitive_descriptor.h"
#include <solve.h>
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;
}

324
application/test_occ.py
View File

@ -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")

117
application/test_occ_draw.py
View File

@ -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()
Write Preview
Loading…
Cancel
Save