👢 Remove LibGEOS dependency

Author: ZenanH

Project.toml

@@ -5,13 +5,11 @@ version = "0.1.0"
 
 [deps]
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
-LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 
 [compat]
 CondaPkg = "0.2"
-LibGEOS = "0.9"
 PythonCall = "0.9"
 julia = "1.11"
 

src/FastPointQuery.jl

@@ -1,6 +1,6 @@
 module FastPointQuery
 
-using CondaPkg, LibGEOS, Printf, PythonCall
+using CondaPkg, Printf, PythonCall
 
 const np      = PythonCall.pynew()
 const shapely = PythonCall.pynew()
@@ -19,11 +19,13 @@ function __init__()
     end
 end
 
+include(joinpath(@__DIR__, "STLInfo.jl"))
 include(joinpath(@__DIR__, "polygon.jl"))
 
+
 # export structs
-export QueryPolygon
+export QueryPolygon, STLInfo
 # export functions
-export get_polygon, pip_query
+export get_polygon, pip_query, loadmesh
 
 end

src/polygon.jl

@@ -21,37 +21,31 @@ struct QueryPolygon
 end
 
 """
-    get_polygon(pts; ratio=0.1)
+    get_polygon(points::AbstractArray; ratio::Bool=0.1)
 
 Description:
 ---
-Get the polygon from the given points. `pts` is a 2xN array, and `ratio` is the
+Get the polygon from the given points. `points` is a 2xN array, and `ratio` is the
 parameter for the concave hull. The default value is 0.1.
+
+Example:
+---
+```julia
+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)
     # inputs checking
     n, m = size(points)
     n == 2 || error("points must be a 2xN array")
     m >= 3 || error("at least 3 points are required")
     ratio > 0 || error("ratio must be positive")
-    allow_holes = 0
-    pts = Array{Float64}(points)
-
-    # convert to GEOS datatype
-    gpts = LibGEOS.MultiPoint([pts[:, i] for i in axes(pts, 2)])
-    ctx  = LibGEOS.get_context(gpts)
-    ptr  = LibGEOS.GEOSConcaveHullByLength_r(ctx, gpts, ratio, allow_holes)
-    conc = LibGEOS.geomFromGEOS(ptr)
-    conc == C_NULL && error("LibGEOS: Error in GEOSConvexHull") 
-
-    # convert back to julia datatype
-    ob_outer = LibGEOS.exteriorRing(conc)
-    cs_outer = LibGEOS.getCoordSeq(ob_outer)
-    xs_outer = LibGEOS.getX(cs_outer)
-    ys_outer = LibGEOS.getY(cs_outer)
-    xy       = vcat(permutedims(xs_outer), permutedims(ys_outer))
-
-    return QueryPolygon(xy, shapely.Polygon(xy'))
+    # use shapely to get the polygon
+    py_points = shapely.MultiPoint(np.array(points'))
+    poly = shapely.concave_hull(py_points, ratio=ratio, allow_holes=false)
+    xy = pyconvert(Array, np.array(poly.exterior.xy)) # get the exterior points of the polygon
+    return QueryPolygon(xy, poly)
 end
 
 """
@@ -68,8 +62,8 @@ Example:
 ```julia
 polypon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
 polygon = get_polygon(polygon_xy, ratio=1) # polygon.py_xy to visualize
-particle_in_polygon(polygon, 0, 0) # check points (0, 0) and return false
-particle_in_polygon(polygon, 0, 0, edge=true) # check points (0, 0) and return true
+pip_query(polygon, 0, 0) # check points (0, 0) and return false
+pip_query(polygon, 0, 0, edge=true) # check points (0, 0) and return true
 ```
 """
 function pip_query(
@@ -92,10 +86,12 @@ array, where each column represents a point (x, y). If `edge` is set to true, it
 the points are on the edge of the polygon as well.
 
 Example:
+---
 ```julia
 polygon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array
 polygon = get_polygon(polygon_xy, ratio=1) # polygon.py_xy to visualize
 pip_query(polygon, [0 0; 0.5 0.5; 1 1], edge=true) # returns [true, true, true]
+```
 """
 function pip_query(
     polygon::QueryPolygon, 
@@ -129,10 +125,12 @@ where each column represents a point (x, y). If `edge` is set to true, it checks
 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, 
@@ -164,3 +162,70 @@ function pip_query(
     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 
+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")
+points = [0 0; 0.5 0.5; 1 1]'
+pip_query(mesh, points, edge=true)
+```
+"""
+function pip_query(
+    mesh    ::Py,
+    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])
+            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])
+            mask = shapely.within(py_points, region)
+            return mask
+        """ => py_pip
+    end
+    tmp = py_pip(mesh, py_points, trimesh, 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)