# Computational Geometry: Regular Distribution of Points inside of a Polygon

**URL:** <https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523>\
**Category:** Numerics\
**Tags:** polygons\
**Created:** [May 15, 2020, 1:43pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523 "2020-05-15T13:43:06Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [May 15, 2020, 1:43pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/1 "2020-05-15T13:43:06Z")

</div>

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:

 ![image](https://global.discourse-cdn.com/julialang/original/2X/7/7a7477f4543cf26b6e440f9f22e175e810fe4722.png)

> **[Generative Pseudo-Random Polygon Fill Patterns in QGIS](https://impermanent.io/2017/05/05/generative-pseudo-random-polygon-fill-patterns-in-qgis/)**
>
> QGIS doesn’t support pseudo-random fill patterns out-of-the-box. However, using the Geometry Generator we can achieve the same effect. Random fill pattern? Here’s a polygon with such a …

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 15, 2020, 2:00pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/2 "2020-05-15T14:00:28Z")

</div>

See also this discussion: [Mesh/grid over convex polytope](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047) (and also [Is there a simple way to generate uniform points within a given region?](https://discourse.julialang.org/t/is-there-a-simple-way-to-generate-uniform-points-within-a-given-region/29800)). The former thread includes [example code](https://discourse.julialang.org/t/mesh-grid-over-convex-polytope/36047/11) that produces a picture very similar to the one you linked.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [May 15, 2020, 2:00pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/3 "2020-05-15T14:00:45Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [May 15, 2020, 2:10pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/4 "2020-05-15T14:10:28Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [May 15, 2020, 4:00pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/5 "2020-05-15T16:00:38Z")

</div>

Okay, so now I have gotten it to work on my polygon:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/2/1269ef5419f781b358e4a48e1af82392086c1476.png)

And now my next step would be to make something like the grey particles in this picture below:

![image](https://global.discourse-cdn.com/julialang/original/3X/4/9/49bf2f7d7ffdbc2e740b768ea1059140ab6b3e31.png)  
[https://www.researchgate.net/publication/256688486\_Particle\_packing\_algorithm\_for\_SPH\_schemes?enrichId=rgreq-b428b49aae0168986c147685bcf97faf-XXX&enrichSource=Y292ZXJQYWdlOzI1NjY4ODQ4NjtBUzo1ODQ0Njg5MTgxOTgyNzJAMTUxNjM1OTY1NzA4NA%3D%3D&el=1\_x\_2&\_esc=publicationCoverPdf](https://www.researchgate.net/publication/256688486_Particle_packing_algorithm_for_SPH_schemes?enrichId=rgreq-b428b49aae0168986c147685bcf97faf-XXX&enrichSource=Y292ZXJQYWdlOzI1NjY4ODQ4NjtBUzo1ODQ0Njg5MTgxOTgyNzJAMTUxNjM1OTY1NzA4NA%3D%3D&el=1_x_2&_esc=publicationCoverPdf)

Would anyone happen to know if there is an easy way using LazySets to do this?

Kind regards

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [May 15, 2020, 6:10pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/6 "2020-05-15T18:10:26Z")

</div>

> [@Ahmed\_Salih](#):
>
> Would anyone happen to know if there is an easy way using LazySets to do this?

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](https://github.com/saschatimme/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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [May 16, 2020, 12:37am UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/7 "2020-05-16T00:37:13Z")

</div>

For future reference, this is what is currently possible:

```julia
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]))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/1/71da8454e6573e6e44659fea647cb5806c4a4a1a.png)

```julia
# 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]))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/5/a51afe68032eb4439cb2269b45a1e20548e3f218.png)

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.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [August 1, 2020, 9:48pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/8 "2020-08-01T21:48:18Z")

</div>

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?

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/b/7b56bef0219f16217eed5fb8de3696d60af72a36.png)

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 1, 2021, 4:45pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/9 "2021-11-01T16:45:06Z")

</div>

Hi @Ahmed_Salih , there is a much simpler api now sitting in Meshes.jl:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [November 1, 2021, 5:20pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/10 "2021-11-01T17:20:51Z")

</div>

@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

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 1, 2021, 10:05pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/11 "2021-11-01T22:05:04Z")

</div>

@juliohm, spent some time reading [Meshes.jl](https://juliageometry.github.io/Meshes.jl/stable/index.html) 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:

```julia
sample(polygon_with_hole, RegularSampling(20, 10))

```

Where `polygon_with_hole` is the donut betwen `(x1,y1)` and `(x2,y2)`:

```julia
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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 1, 2021, 10:32pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/12 "2021-11-01T22:32:13Z")

</div>

@rafael.guerra take a look at the test suite for the sampling methods: [https://github.com/JuliaGeometry/Meshes.jl/blob/master/test/sampling.jl](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 👍🏽

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 2, 2021, 12:25am UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/13 "2021-11-02T00:25:16Z")

</div>

Fyi, a regular grid on a donut can be produced using `PolygonOps`:

 ![PolygonOps_regular_grid_on_donut](https://global.discourse-cdn.com/julialang/original/3X/2/4/24e219f7af68ef73a645b61aba8955e19c3665ef.png)

> **MWE**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 2, 2021, 12:29am UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/14 "2021-11-02T00:29:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 2, 2021, 12:32am UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/15 "2021-11-02T00:32:16Z")

</div>

LOL 🙂, you could have said so!

_ **NB:** many practical uses of this require doing it only once and then speed is irrelevant._

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [December 12, 2021, 2:14am UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/16 "2021-12-12T02:14:03Z")

</div>

@juliohm Can Mehes.jl do the point-in-polyhedron test? Assume that the polyhedron is strictly convex.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 12, 2021, 12:01pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/17 "2021-12-12T12:01:38Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![jacobusmmsmit](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacobusmmsmit/32/217669_2.png) [@jacobusmmsmit](https://discourse.julialang.org/u/jacobusmmsmit)\
**Post date:** [December 12, 2021, 1:41pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/18 "2021-12-12T13:41:36Z")

</div>

@juliohm @liuyxpp I stole some code from [a python user](https://github.com/ranjeethmahankali/SampleCode/blob/master/PointInPolyhedron/PointInMesh.py) who implemented the technique from [this paper](https://www.sciencedirect.com/science/article/pii/0734189X84901336). 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)`

```julia
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

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 12, 2021, 1:46pm UTC](https://discourse.julialang.org/t/computational-geometry-regular-distribution-of-points-inside-of-a-polygon/39523/19 "2021-12-12T13:46:57Z")

</div>

@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?