|
| 1 | +#==========================================================================================++ |
| 2 | +| TABLE OF CONTENTS: | |
| 3 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 4 | +| struct list: | |
| 5 | +| - QueryPolygon | |
| 6 | ++------------------------------------------------------------------------------------------+ |
| 7 | +| function list: | |
| 8 | +| - get_polygon (from LibGEOS, need to improve) | |
| 9 | +| - pip_query | |
| 10 | ++==========================================================================================# |
| 11 | + |
| 12 | +struct QueryPolygon |
| 13 | + ju_xy::AbstractMatrix |
| 14 | + py_xy::Py |
| 15 | + function QueryPolygon(ju_xy::AbstractMatrix, py_xy::Py) |
| 16 | + n, m = size(ju_xy) |
| 17 | + m >= 3 || error("at least 3 points are required") |
| 18 | + n == 2 || error("points must be 2D (2×N array)") |
| 19 | + new(ju_xy, py_xy) |
| 20 | + end |
| 21 | +end |
| 22 | + |
| 23 | +""" |
| 24 | + get_polygon(pts; ratio=0.1) |
| 25 | +
|
| 26 | +Description: |
| 27 | +--- |
| 28 | +Get the polygon from the given points. `pts` is a 2xN array, and `ratio` is the |
| 29 | +parameter for the concave hull. The default value is 0.1. |
| 30 | +""" |
| 31 | +function get_polygon(points::AbstractArray; ratio::Real=0.1) |
| 32 | + # inputs checking |
| 33 | + n, m = size(points) |
| 34 | + n == 2 || error("points must be a 2xN array") |
| 35 | + m >= 3 || error("at least 3 points are required") |
| 36 | + ratio > 0 || error("ratio must be positive") |
| 37 | + allow_holes = 0 |
| 38 | + pts = Array{Float64}(points) |
| 39 | + |
| 40 | + # convert to GEOS datatype |
| 41 | + gpts = LibGEOS.MultiPoint([pts[:, i] for i in axes(pts, 2)]) |
| 42 | + ctx = LibGEOS.get_context(gpts) |
| 43 | + ptr = LibGEOS.GEOSConcaveHullByLength_r(ctx, gpts, ratio, allow_holes) |
| 44 | + conc = LibGEOS.geomFromGEOS(ptr) |
| 45 | + conc == C_NULL && error("LibGEOS: Error in GEOSConvexHull") |
| 46 | + |
| 47 | + # convert back to julia datatype |
| 48 | + ob_outer = LibGEOS.exteriorRing(conc) |
| 49 | + cs_outer = LibGEOS.getCoordSeq(ob_outer) |
| 50 | + xs_outer = LibGEOS.getX(cs_outer) |
| 51 | + ys_outer = LibGEOS.getY(cs_outer) |
| 52 | + xy = vcat(permutedims(xs_outer), permutedims(ys_outer)) |
| 53 | + |
| 54 | + return QueryPolygon(xy, shapely.Polygon(xy')) |
| 55 | +end |
| 56 | + |
| 57 | +""" |
| 58 | + pip_query(polygon::QueryPolygon, px::Real, py::Real; edge::Bool=false) |
| 59 | +
|
| 60 | +Description: |
| 61 | +--- |
| 62 | +Determine whether a point [px, py] is inside a polygon. Note to use the function |
| 63 | +`get_polygon` to get it; otherwise, it may lead to incorrect results. If the point is on the |
| 64 | +edge of the polygon, you can set `edge=true` to check if it is considered as inside. |
| 65 | +
|
| 66 | +Example: |
| 67 | +--- |
| 68 | +```julia |
| 69 | +polypon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array |
| 70 | +polygon = get_polygon(polygon_xy, ratio=1) # polygon.py_xy to visualize |
| 71 | +particle_in_polygon(polygon, 0, 0) # check points (0, 0) and return false |
| 72 | +particle_in_polygon(polygon, 0, 0, edge=true) # check points (0, 0) and return true |
| 73 | +``` |
| 74 | +""" |
| 75 | +function pip_query( |
| 76 | + polygon::QueryPolygon, |
| 77 | + px ::Real, |
| 78 | + py ::Real; |
| 79 | + edge ::Bool=false |
| 80 | +) |
| 81 | + func = edge ? shapely.intersects_xy : shapely.contains_xy |
| 82 | + return pyconvert(Bool, func(polygon.py_xy, px, py).item()) |
| 83 | +end |
| 84 | + |
| 85 | +""" |
| 86 | + pip_query(polygon::QueryPolygon, points::AbstractMatrix; edge::Bool=false) |
| 87 | +
|
| 88 | +Description: |
| 89 | +--- |
| 90 | +Determine whether a set of points is inside a polygon. The input `points` should be a 2xN |
| 91 | +array, where each column represents a point (x, y). If `edge` is set to true, it checks if |
| 92 | +the points are on the edge of the polygon as well. |
| 93 | +
|
| 94 | +Example: |
| 95 | +```julia |
| 96 | +polygon_xy = [0 0; 1 0; 1 1; 0 1]' # 2xN array |
| 97 | +polygon = get_polygon(polygon_xy, ratio=1) # polygon.py_xy to visualize |
| 98 | +pip_query(polygon, [0 0; 0.5 0.5; 1 1], edge=true) # returns [true, true, true] |
| 99 | +""" |
| 100 | +function pip_query( |
| 101 | + polygon::QueryPolygon, |
| 102 | + points ::AbstractMatrix; |
| 103 | + edge ::Bool=false |
| 104 | +) |
| 105 | + # inputs check for points |
| 106 | + n, m = size(points) |
| 107 | + m ≥ 1 || error("points should have at least 1 point") |
| 108 | + n == 2 || error("points should be a 2xN array") |
| 109 | + py_points = shapely.points(points') |
| 110 | + |
| 111 | + # check if the particle is inside the polygon |
| 112 | + if edge |
| 113 | + mask = shapely.covers(polygon.py_xy, py_points) |
| 114 | + else |
| 115 | + mask = shapely.within(py_points, polygon.py_xy) |
| 116 | + end |
| 117 | + |
| 118 | + return pyconvert(Vector{Bool}, mask) |
| 119 | +end |
| 120 | + |
| 121 | +""" |
| 122 | + pip_query(stl_file::String, points::AbstractMatrix; edge::Bool=false) |
| 123 | +
|
| 124 | +Description: |
| 125 | +--- |
| 126 | +Determine whether a set of points is inside a 3D mesh defined by an STL file. The input |
| 127 | +`stl_file` should be a valid file path to an STL file, and `points` should be a 2xN array |
| 128 | +where each column represents a point (x, y). If `edge` is set to true, it checks if the |
| 129 | +points are on the edge of the mesh as well. |
| 130 | +
|
| 131 | +Example: |
| 132 | +```julia |
| 133 | +stl_file = "path/to/your/file.stl" |
| 134 | +points = [0 0; 0.5 0.5; 1 1]' |
| 135 | +pip_query(stl_file, points, edge=true) |
| 136 | +""" |
| 137 | +function pip_query( |
| 138 | + stl_file::String, |
| 139 | + points ::AbstractMatrix; |
| 140 | + edge ::Bool=false |
| 141 | +) |
| 142 | + isfile(stl_file) || error("stl_file must be a valid file path") |
| 143 | + size(points, 1) == 2 || error("points must be a 2xN array") |
| 144 | + py_points = shapely.points(points') |
| 145 | + if edge |
| 146 | + @pyexec """ |
| 147 | + def py_pip(stl_file, py_points, trimesh, shapely): |
| 148 | + mesh = trimesh.load(stl_file, force="mesh") |
| 149 | + tris2d = mesh.triangles[:, :, :2] |
| 150 | + region = shapely.unary_union([shapely.Polygon(t) for t in tris2d]) |
| 151 | + mask = shapely.covers(region, py_points) |
| 152 | + return mask |
| 153 | + """ => py_pip |
| 154 | + else |
| 155 | + @pyexec """ |
| 156 | + def py_pip(stl_file, py_points, trimesh, shapely): |
| 157 | + mesh = trimesh.load(stl_file, force="mesh") |
| 158 | + tris2d = mesh.triangles[:, :, :2] |
| 159 | + region = shapely.unary_union([shapely.Polygon(t) for t in tris2d]) |
| 160 | + mask = shapely.within(py_points, region) |
| 161 | + return mask |
| 162 | + """ => py_pip |
| 163 | + end |
| 164 | + tmp = py_pip(stl_file, py_points, trimesh, shapely) |
| 165 | + return pyconvert(Vector{Bool}, tmp) |
| 166 | +end |
0 commit comments