Hello!
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:
The link includes some code based on a Python package to achieve this - I wondered if someone might have made a package like this for Julia?
Kind regards
See also this discussion: Mesh/grid over convex polytope (and also Is there a simple way to generate uniform points within a given region?). The former thread includes example code that produces a picture very similar to the one you linked.
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 (Processes · GeoStats.jl) 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.
Oh wauw thanks to both of you.
@stevengj I will be looking at that example now thanks.
@juliohm seems interesting and thank you very much for the correct terminology! Makes future searches much easier.
Kind regards
Okay, so now I have gotten it to work on my polygon:
And now my next step would be to make something like the grey particles in this picture below:
Would anyone happen to know if there is an easy way using LazySets to do this?
Kind regards
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]))
# superposition of two Binomial processes
p₁ = BinomialProcess(20)
p₂ = BinomialProcess(80)
p = p₁ ∪ p₂ # 100 points
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.
Hi again!
I am playing a bit around with the PoissonProcess. When I have generated “s”, how do I go about extracting the number of elements?
I think it is in the type, so that is why length/size does not work, could you explain it to me?
EDIT: Figured it out, one uses s[1].coords then size works as usual, my bad
Kind regards
Hi @Ahmed_Salih , there is a much simpler api now sitting in Meshes.jl:
points = sample(polygon, MinDistanceSampling(1.0))
This will handle non-convex polygons and holes. Sorry for the delay, I missed your reply because you didn’t tag my nick name.
@juliohm no problem at all, it should work without tagging you at all, but just tagging you know to say thanks for letting me know!
It is really cool to see that these kind of operations are being put into some great Julia packages.
Kind regards
@juliohm, spent some time reading Meshes.jl documentation but could not find easily the solution to OP’s donut regular gridding problem.
Could you please provide some tips on how to proceed and grid with say dx=20 and dy=10 using function:
sample(polygon_with_hole, RegularSampling(20, 10))
Where polygon_with_hole is the donut betwen (x1,y1) and (x2,y2):
using Meshes
t = LinRange(0, 2π, 100)
x1, y1 = 100*cos.(t), 50*sin.(t)
x1[end] = x1[1]
y1[end] = y1[1] # needed due to floating-point errors?
x2 = 0.7*x1
y2 = 0.7*y1
points1 = Point2[tuple.(x1,y1)...]
points2 = Point2[tuple.(x2,y2)...]
# how to define polygon_with_hole?
sample(polygon_with_hole, RegularSampling(20, 10))
Thanks in advance.
@rafael.guerra take a look at the test suite for the sampling methods: https://github.com/JuliaGeometry/Meshes.jl/blob/master/test/sampling.jl I am happy to review a PR with new examples for the sampling methods in the docs.
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 
Fyi, a regular grid on a donut can be produced using PolygonOps:
using PolygonOps, StaticArrays
t = LinRange(0, 2π, 100)
x1, y1 = 500*cos.(t), 250*sin.(t)
x1[end], y1[end] = x1[1], y1[1] # needed due to floating-point errors
x2, y2 = 0.7*x1, 0.7*y1
polygon1 = SVector.(x1, y1)
polygon2 = 0.7*polygon1
xr, yr = range(extrema(x1)..., step=20), range(extrema(y1)..., step=10)
points = Iterators.product(xr, yr)
mask1 = [inpolygon(p, polygon1; in=true, on=true, out=false) for p in points]
mask2 = [inpolygon(p, polygon2; in=false, on=true, out=true) for p in points]
finalgrid = collect(points)[mask1 .* mask2]
using Plots; gr(dpi=600)
plot(x1, y1, legend=false, lc=:gray, lw=0.3, ratio=1)
plot!(x2, y2, lc=:gray, lw=0.3)
scatter!(first.(finalgrid), last.(finalgrid), m=1, msw=0, mc=:red)
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.
LOL 
, you could have said so!
NB: many practical uses of this require doing it only once and then speed is irrelevant.
@juliohm Can Mehes.jl do the point-in-polyhedron test? Assume that the polyhedron is strictly convex.
@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?