So I am basically asking since I struggle to figure out if there is an efficient way to do this. Basically I want to be able to take any arbitrary polygon and fill it with particles as shown in this figure:

This is an upcoming feature of GeoStats.jl, particularly the PointPatterns.jl submodule. The problem of sampling points with certain characteristics inside a geometry is known in the literature as point pattern analysis and synthesis. Currently, we have some methods to sample inside hyperrectangles (https://juliaearth.github.io/GeoStats.jl/stable/pointpatterns/pointprocs) but to generalize these methods to arbitrary polygons we need to first address a major milestone in the project regarding general finite element meshes. It is the next item in my TODO list. Meanwhile, I think you could achieve a similar result by generating random numbers in the bounding box of the polygon, and then taking an intersection. @visr and @evetion may have thoughts on this as well.

There is no internal representation for the red curve (“solid boundary profile”) in LazySets. Not sure that i understood what you are trying to achieve, but to draw implicit plots like these you can use ImplicitPlots.jl. If you know how to test if the point is inside/outside the given region, then you can use the same idea of rejection sampling from the other discussion and make such plot.

For future reference, this is what is currently possible:

using GeoStats
using Plots
# rectangular region of interest
r = RectangleRegion((0.,0.), (100.,100.))
# Poisson process with λ=0.1 intensity
p = PoissonProcess(0.1)
# two samples from the process
s = rand(p, r, 2)
plot(plot(s[1]), plot(s[2]))

Notice that this method works in N-dimensional regions. The plan is to handle other types of regions (e.g. polygons, polyhedrons) in the next minor release.

Also, notice that the sampling method I suggested produces points that are well-spaced given a minimum distance between any two points. If the OP needs points that are perfectly aligned in a grid, then we need to add a new method for RegularSampling with Polygon. Contributions are welcome

We have the equivalent method in Meshes.jl. Instead of PolygonOps.inpolygon you can use point \in polygon. The problem with this solution is that it is not fast enough for large polygonal areas. The MinDistanceSampling method I suggested above is.

@liuyxpp we have very few algorithms implemented for polyhedron at this point, but the building blocks are there. If you have the skills, please consider submitting a PR. We are continuously adding new features and that would be super helpful to the community as a whole.

@juliohm@liuyxpp I stole some code from a python user who implemented the technique from this paper. I have no familiarity with Meshes.jl or the geometry ecosystem so the only function I didn’t know how to implement (or whether it would be necessary) is the badly named convert_polyhedron_to_vector_of_triangles(polyhedron)

using LinearAlgebra
# calculates the solid angle between three vectors
# reference: https://math.stackexchange.com/questions/202261/solid-angle-between-vectors-in-n-dimensional-space
function get_solid_angle(a,b,c)
a, b, c = normalize.((a, b, c))
numer = (a × b)⋅ c
denom = 1 + (a ⋅ b) + (b ⋅ c) + (c ⋅ a)
angle = 2 * atan(numer, denom)
return abs(angle)
end
# calculating the normal for a triangular face
normal(p, q, r) = (q - p) × (r - q)
# calculating the centroid of the triangular face
centroid(p,q,r) = (p + q + r)/3
# the logic of the point in mesh algorithm,
# assumes a vector of triangles has structure like
# Vector{Tuple{Point, Point, Point}
# where Point = SVector{3, Float64}
function inpolyhedron(polyhedron, pt, tolerance = 1e-6)
# triangle_mesh = convert_polyhedron_to_vector_of_triangles(polyhedron)
total_angle = 0
for triangle in triangle_mesh
ptA, ptB, ptC = triangle
a = ptA .- pt
b = ptB .- pt
c = ptC .- pt
angle = get_solid_angle(a,b,c)
normalvector = normal(ptA, ptB, ptC)
centre = centroid(ptA, ptB, ptC)
facevec = pt - centre
dot = normalvector ⋅ facevec
factor = dot > 0 ? 1 : -1
total_angle += angle * factor
end
abs_total = abs(total_angle)
return abs(abs_total - (4*π)) < tolerance
end

@jacobusmmsmit a PR is welcome. We already have the triangle functions. We only need need to define the boundary(polyhedron) (which returns a mesh of n-gon) and use the existing functions in the project. Would you like to work on that?