Skip to content

Commit 4c85fdf

Browse files
ZenanHZenanH
committed
✨ Add meshbuilder functions for 2D and 3D mesh generation and update utils.jl exports
1 parent bc88b0d commit 4c85fdf

2 files changed

Lines changed: 172 additions & 0 deletions

File tree

‎src/discretization/_polyhedron.jl‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,4 +27,76 @@ function _get_pts_voxel(stl_model::STLInfo3D, h::Real; fill::Bool=true)
2727
else
2828
return PyArray(pts_cen)
2929
end
30+
end
31+
32+
function _get_pts_ray(stl_model::STLInfo3D, h::Real; fill::Bool=true, ϵ::String="FP32")
33+
points, pts2 = prepareprojection(stl_model, h, ϵ)
34+
edges = projectionlist(stl_model, points)
35+
pts_cen = fill_particles(edges, pts2, h)
36+
if fill
37+
return populate_pts(pts_cen, h)
38+
else
39+
return pts_cen
40+
end
41+
end
42+
43+
function prepareprojection(stl_model::STLInfo3D, h::Real, ϵ)
44+
T = ϵ == "FP32" ? Float32 : Float64
45+
pts2 = meshbuilder(stl_model.vmin[1]-1.5h : h : stl_model.vmax[1]+1.5h,
46+
stl_model.vmin[2]-1.5h : h : stl_model.vmax[2]+1.5h, ϵ=ϵ)
47+
vzlimit = T(stl_model.vmin[3] - 10h)
48+
points = vcat(pts2, vzlimit .* ones(T, 1, size(pts2, 2)))
49+
return points, pts2
50+
end
51+
52+
function projectionlist(stl_model::STLInfo3D, points::AbstractArray{T}) where T
53+
# inputs check
54+
pts = points'; n, m = size(pts)
55+
m == 3 || error("points must be a 3xN array")
56+
57+
mesh = stl_model.mesh
58+
pts = np.asarray(pts, dtype=np.float32)
59+
n_rays = pts.shape[0]
60+
scene = o3d.t.geometry.RaycastingScene()
61+
tmesh = o3d.t.geometry.TriangleMesh.from_legacy(mesh, vertex_dtype=o3d.core.float32)
62+
scene.add_triangles(tmesh)
63+
64+
dirs = np.tile([0, 0, 1], (n_rays, 1))
65+
rays = o3d.core.Tensor(np.hstack((pts, dirs)), o3d.core.float32)
66+
hits = scene.list_intersections(rays)
67+
t_hit = hits["t_hit"].numpy()
68+
splits = hits["ray_splits"].numpy()
69+
70+
t_chunks = np.split(t_hit, splits[1:-1])
71+
@pyexec """
72+
def py_tmp(pts, t_chunks, n_rays, splits, t_hit):
73+
z_chunks = [pts[i, 2] + t_hit[splits[i] : splits[i + 1]] for i in range(n_rays)]
74+
return z_chunks
75+
""" => py_tmp
76+
result = py_tmp(pts, t_chunks, n, splits, t_hit)
77+
return py2ju(Vector{Vector{T}}, result)
78+
end
79+
80+
function fill_particles(edges::AbstractArray, pts::AbstractMatrix{T}, h::Real) where T
81+
x_all = Vector{T}(); sizehint!(x_all, 1024)
82+
y_all = Vector{T}(); sizehint!(y_all, 1024)
83+
z_all = Vector{T}(); sizehint!(z_all, 1024)
84+
@inbounds for i in eachindex(edges)
85+
v = edges[i]
86+
if !isempty(v) && iseven(length(v)) && any(!iszero, v)
87+
sort!(v); xi, yi = pts[1, i], pts[2, i]
88+
for j in 1:2:(length(v) - 1)
89+
for z in v[j]:h:v[j + 1]
90+
push!(x_all, xi)
91+
push!(y_all, yi)
92+
push!(z_all, z)
93+
end
94+
end
95+
end
96+
end
97+
xyz = Matrix{T}(undef, 3, length(x_all))
98+
xyz[1, :] .= x_all
99+
xyz[2, :] .= y_all
100+
xyz[3, :] .= z_all
101+
return xyz
30102
end

‎src/utils.jl‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,13 @@
22
| TABLE OF CONTENTS: |
33
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
44
| - get_polygon |
5+
| - meshbuilder |
6+
| - populate_pts |
57
+==========================================================================================#
68

79
export get_polygon
10+
export meshbuilder
11+
export populate_pts
812

913
"""
1014
get_polygon(points::AbstractMatrix; ratio::Bool=0.1, allow_holes::Bool=false)
@@ -32,4 +36,100 @@ function get_polygon(points::AbstractMatrix; ratio::Real=0.1, allow_holes::Bool=
3236
py_points = shapely.MultiPoint(np.array(points'))
3337
poly = shapely.concave_hull(py_points, ratio=ratio, allow_holes=allow_holes)
3438
return QueryPolygon(poly)
39+
end
40+
41+
"""
42+
meshbuilder(x::T, y::T; ϵ::String="FP64") where T <: AbstractRange
43+
44+
Description:
45+
---
46+
Generate structured mesh in 2D space.
47+
48+
Example:
49+
---
50+
```julia
51+
x = 0:0.1:1
52+
y = 0:0.1:1
53+
mesh = meshbuilder(x, y; ϵ="FP32") # returns a 2xN array of points
54+
```
55+
"""
56+
function meshbuilder(x::T, y::T; ϵ::String="FP64") where T <: AbstractRange
57+
x_tmp = repeat(x', length(y), 1) |> vec
58+
y_tmp = repeat(y , 1, length(x)) |> vec
59+
T1 = ϵ == "FP32" ? Float32 : Float64
60+
return Array{T1}(hcat(x_tmp, y_tmp)')
61+
end
62+
63+
"""
64+
meshbuilder(x::T, y::T, z::T; ϵ::String="FP64") where T <: AbstractRange
65+
66+
Description:
67+
---
68+
Generate structured mesh in 3D space.
69+
70+
Example:
71+
---
72+
```julia
73+
x = 0:0.1:1
74+
y = 0:0.1:1
75+
z = 0:0.1:1
76+
mesh = meshbuilder(x, y, z; ϵ="FP32") # returns a 3xN array of points
77+
"""
78+
function meshbuilder(x::T, y::T, z::T; ϵ::String="FP64") where T <: AbstractRange
79+
vx = x |> collect
80+
vy = y |> collect
81+
vz = z |> collect
82+
m, n, o = length(vy), length(vx), length(vz)
83+
vx = reshape(vx, 1, n, 1)
84+
vy = reshape(vy, m, 1, 1)
85+
vz = reshape(vz, 1, 1, o)
86+
om = ones(Int, m)
87+
on = ones(Int, n)
88+
oo = ones(Int, o)
89+
x_tmp = vec(vx[om, :, oo])
90+
y_tmp = vec(vy[:, on, oo])
91+
z_tmp = vec(vz[om, on, :])
92+
T1 = ϵ == "FP32" ? Float32 : Float64
93+
return Array{T1}(hcat(x_tmp, y_tmp, z_tmp)')
94+
end
95+
96+
"""
97+
populate_pts(pts_cen::AbstractMatrix{T}, h::Real) where T
98+
99+
Description:
100+
---
101+
Populate the points around the center points `pts_cen` with the spacing `h/4` (2/3D).
102+
"""
103+
@views function populate_pts(pts_cen::AbstractMatrix{T}, h::Real) where {T}
104+
h > 0 || error("h must be positive")
105+
D, N = size(pts_cen)
106+
(D == 2 || D == 3) ||
107+
throw(ArgumentError("pts_cen must have 2 or 3 rows (got $D)"))
108+
oft = T(h * 0.25)
109+
110+
if D == 3
111+
offsets = oft .* [
112+
-1 -1 -1 -1 1 1 1 1;
113+
-1 -1 1 1 -1 -1 1 1;
114+
-1 1 -1 1 -1 1 -1 1
115+
]
116+
K = 8
117+
else
118+
offsets = oft .* [
119+
-1 -1 1 1;
120+
-1 1 -1 1
121+
]
122+
K = 4
123+
end
124+
125+
pts = Matrix{T}(undef, D, K*N)
126+
@inbounds for i in 1:N
127+
base = (i-1)*K
128+
cen = pts_cen[:, i] # 取第 i 个中心点 (D×1)
129+
@simd for j in 1:K
130+
pts[:, base+j] .= cen .+ offsets[:, j]
131+
end
132+
end
133+
134+
return pts
35135
end

0 commit comments

Comments
 (0)