🙆🏻 Refactor STL handling

Author: ZenanH

src/STLInfo.jl

@@ -1,27 +1,12 @@
-#==========================================================================================++
+#==========================================================================================+
 | TABLE OF CONTENTS:                                                                       |
 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-| struct list:                                                                             |
-|  - QueryPolygon                                                                          |
-+------------------------------------------------------------------------------------------+
-| function list:                                                                           |
-|  - get_polygon (from LibGEOS, need to improve)                                           |
+|  - get_polygon                                                                           |
 |  - pip_query                                                                             |
 +==========================================================================================#
 
-struct QueryPolygon
-    ju_xy::AbstractMatrix
-    py_xy::Py
-    function QueryPolygon(ju_xy::AbstractMatrix, py_xy::Py)
-        n, m = size(ju_xy)
-        m >= 3 || error("at least 3 points are required")
-        n == 2 || error("points must be 2D (2×N array)")
-        new(ju_xy, py_xy)
-    end
-end
-
 """
-    get_polygon(points::AbstractArray; ratio::Bool=0.1)
+    get_polygon(points::AbstractMatrix; ratio::Bool=0.1)
 
 Description:
 ---
@@ -35,7 +20,7 @@ polygon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
 polygon = get_polygon(polygon_xy, ratio=1) # polygon.py_xy to visualize
 ```
 """
-function get_polygon(points::AbstractArray; ratio::Real=0.1)
+function get_polygon(points::AbstractMatrix; ratio::Real=0.1)
     # inputs checking
     n, m = size(points)
     n == 2 || error("points must be a 2xN array")
@@ -115,117 +100,49 @@ function pip_query(
 end
 
 """
-    pip_query(stl_file::String, points::AbstractMatrix; edge::Bool=false)
+    pip_query(stl_model::STLInfo2D, points::AbstractMatrix; edge::Bool=false)
 
 Description:
 ---
-Determine whether a set of points is inside a 3D mesh defined by an STL file. The input
-`stl_file` should be a valid file path to an STL file, and `points` should be a 2xN array
-where each column represents a point (x, y). If `edge` is set to true, it checks if the 
-points are on the edge of the mesh as well.
-
-Example:
----
-```julia
-stl_file = "path/to/your/file.stl"
-points = [0 0; 0.5 0.5; 1 1]'
-pip_query(stl_file, points, edge=true)
-```
-"""
-function pip_query(
-    stl_file::String, 
-    points  ::AbstractMatrix;
-    edge    ::Bool=false
-)
-    isfile(stl_file) || error("stl_file must be a valid file path")
-    size(points, 1) == 2 || error("points must be a 2xN array")
-    py_points = shapely.points(points')
-    if edge
-        @pyexec """
-        def py_pip(stl_file, py_points, trimesh, shapely):
-            mesh = trimesh.load(stl_file, force="mesh")
-            tris2d = mesh.triangles[:, :, :2]
-            region = shapely.unary_union([shapely.Polygon(t) for t in tris2d])
-            mask = shapely.covers(region, py_points)
-            return mask
-        """ => py_pip
-    else
-        @pyexec """
-        def py_pip(stl_file, py_points, trimesh, shapely):
-            mesh = trimesh.load(stl_file, force="mesh")
-            tris2d = mesh.triangles[:, :, :2]
-            region = shapely.unary_union([shapely.Polygon(t) for t in tris2d])
-            mask = shapely.within(py_points, region)
-            return mask
-        """ => py_pip
-    end
-    tmp = py_pip(stl_file, py_points, trimesh, shapely)
-    return pyconvert(Vector{Bool}, tmp)
-end
-
-"""
-    pip_query(mesh::Py, points::AbstractMatrix; edge::Bool=false)
-
-Description:
----
-Determine whether a set of points is inside a 3D mesh defined by a Python object. The input 
-`mesh` should be a Python object representing the mesh (`trimesh`), and `points` should be 
+Determine whether a set of points is inside a 3D mesh defined by a 'STLInfo'. The input 
+`stl_model` should be initiated by function `readSTL2D(stl_file)`, and `points` should be 
 a 2xN array where each column represents a point (x, y). If `edge` is set to true, it checks 
 if the points are on the edge of the mesh as well.
 
 Example:
 ---
 ```julia
-mesh = trimesh.load("path/to/your/file.stl", force="mesh")
+stl_model = readSTL2D("path/to/your/file.stl")
 points = [0 0; 0.5 0.5; 1 1]'
-pip_query(mesh, points, edge=true)
+pip_query(stl_model, points, edge=true)
 ```
 """
 function pip_query(
-    mesh    ::Py,
-    points  ::AbstractMatrix;
-    edge    ::Bool=false
+    stl_model::STLInfo2D,
+    points   ::AbstractMatrix;
+    edge     ::Bool=false
 )
     size(points, 1) == 2 || error("points must be a 2xN array")
     py_points = shapely.points(points')
     if edge
         @pyexec """
-        def py_pip(mesh, py_points, trimesh, shapely):
-            tris2d = mesh.triangles[:, :, :2]
-            region = shapely.unary_union([shapely.Polygon(t) for t in tris2d])
+        def py_pip(vertices, triangles, py_points, shapely):
+            triangle_coords_2d = vertices[triangles][:, :, :2]
+            tris2d = shapely.polygons(triangle_coords_2d)
+            region = shapely.unary_union(tris2d)
             mask = shapely.covers(region, py_points)
             return mask
         """ => py_pip
     else
         @pyexec """
-        def py_pip(mesh, py_points, trimesh, shapely):
-            tris2d = mesh.triangles[:, :, :2]
-            region = shapely.unary_union([shapely.Polygon(t) for t in tris2d])
+        def py_pip(vertices, triangles, py_points, shapely):
+            triangle_coords_2d = vertices[triangles][:, :, :2]
+            tris2d = shapely.polygons(triangle_coords_2d)
+            region = shapely.unary_union(tris2d)
             mask = shapely.within(py_points, region)
             return mask
         """ => py_pip
     end
-    tmp = py_pip(mesh, py_points, trimesh, shapely)
+    tmp = py_pip(stl_model.py_vertices, stl_model.py_triangles, py_points, shapely)
     return pyconvert(Vector{Bool}, tmp)
-end
-
-"""
-    pip_query(stl_model::STLInfo, points::AbstractMatrix; edge::Bool=false)
-
-Description:
----
-Determine whether a set of points is inside a 3D mesh defined by a 'STLInfo'. The input 
-`stl_model` should be initiated by function `loadmesh(stl_file)`, and `points` should be 
-a 2xN array where each column represents a point (x, y). If `edge` is set to true, it checks 
-if the points are on the edge of the mesh as well.
-
-Example:
----
-```julia
-stl_model = loadmesh("path/to/your/file.stl")
-points = [0 0; 0.5 0.5; 1 1]'
-pip_query(stl_model, points, edge=true)
-```
-"""
-pip_query(stl_model::STLInfo, points::AbstractMatrix; edge::Bool=false) = pip_query(
-    stl_model.mesh, points; edge=edge)
+end

src/polygon.jl

@@ -0,0 +1,35 @@
+#==========================================================================================+
+| TABLE OF CONTENTS:                                                                       |
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+|  - pip_query                                                                             |
++==========================================================================================#
+
+function pip_query(
+    stl_model::STLInfo3D, 
+    points   ::AbstractMatrix; 
+    nsamples ::Int=3, 
+    dev      ::String="CPU:0"
+)
+    # inputs check for points
+    n, m = size(points)
+    n == 3 || error("points must be a 3xN array")
+    m >= 1 || error("at least 1 points are required")
+    py_points = np.array(points')
+
+    # query by using open3d
+    vertices = stl_model.py_vertices
+    faces = stl_model.py_triangles
+    _dev = o3d.core.Device(dev)
+    verts_o3d = o3d.core.Tensor(vertices, dtype=o3d.core.Dtype.Float32, device=_dev)
+    faces_o3d = o3d.core.Tensor(faces, dtype=o3d.core.Dtype.Int32, device=_dev)
+    mesh_t = o3d.t.geometry.TriangleMesh(_dev)
+    mesh_t.vertex.positions = verts_o3d
+    mesh_t.triangle.indices = faces_o3d
+    scene = o3d.t.geometry.RaycastingScene(device=_dev)
+    _ = scene.add_triangles(mesh_t)
+    pts_t = o3d.core.Tensor(py_points.astype(np.float32), device=_dev)
+    occ = scene.compute_occupancy(pts_t, nsamples=nsamples)
+    tmp = occ.numpy().astype(pybuiltins.bool)
+
+    return pyconvert(Vector{Bool}, tmp)
+end

src/polyhedron.jl

@@ -0,0 +1,55 @@
+#==========================================================================================+
+| TABLE OF CONTENTS:                                                                       |
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+|  - readSTL2D                                                                             |
+|  - readSTL3D                                                                             |
++==========================================================================================#
+
+"""
+    readSTL2D(stl_file::String)
+
+Description:
+---
+Read a STL file as a 2D object (z=0) containing the mesh, bounding box, vertices, and 
+triangles.
+
+Example:
+---
+```julia
+stl_info = readSTL2D("path/to/your/file.stl")
+```
+"""
+function readSTL2D(stl_file)
+    isfile(stl_file) || error("stl_file should be a valid file path")
+    mesh = o3d.io.read_triangle_mesh(stl_file)
+    vertices = np.asarray(mesh.vertices)           # (N, 3)
+    triangles = np.asarray(mesh.triangles)
+    aabb = mesh.get_axis_aligned_bounding_box()
+    vmin = pyconvert(Vector, aabb.get_min_bound())
+    vmax = pyconvert(Vector, aabb.get_max_bound())
+    return STLInfo2D(mesh, vmin, vmax, vertices, triangles) 
+end
+
+"""
+    readSTL3D(stl_file::String)
+
+Description:
+---
+Read a STL file as a 3D object containing the mesh, bounding box, vertices, and triangles.
+
+Example:
+---
+```julia
+stl_info = readSTL3D("path/to/your/file.stl")
+```
+"""
+function readSTL3D(stl_file)
+    isfile(stl_file) || error("stl_file should be a valid file path")
+    mesh = o3d.io.read_triangle_mesh(stl_file)
+    vertices = np.asarray(mesh.vertices)
+    triangles = np.asarray(mesh.triangles)
+    aabb = mesh.get_axis_aligned_bounding_box()
+    vmin = pyconvert(Vector, aabb.get_min_bound())
+    vmax = pyconvert(Vector, aabb.get_max_bound())
+    return STLInfo3D(mesh, vmin, vmax, vertices, triangles) 
+end