TriangulateMesh is almost the right tool for this problem. For example,

numb = 113;
region = BoundaryMeshRegion[{{0, 0}, {3, 0}, {3, 1}, {0, 1}}, Line[{1, 2, 3, 4, 1}]];
mesh = TriangulateMesh[region, MaxCellMeasure :> RegionMeasure[region]/numb]

meshed rectangle

Then I can do the following:

points = MeshCoordinates @ mesh;
Graphics @ Point @ points

points in rectangle

Length @ points

However, 113 points were desired, and this approach will always (or almost always) find slightly fewer points than what were desired. Is there an efficient way to get exactly the desired number of points? Perhaps we can spelunk the TriangulateMesh algorithm and use undocumented-internal functions to subdivide the largest cells to get the desired number of points.

  • $\begingroup$ Several methods are shown in this earlier MSE thread $\endgroup$ Commented Jan 29, 2017 at 15:27
  • $\begingroup$ This wolfram blog post by @silvia is definitely worth a look $\endgroup$
    – gpap
    Commented Jan 30, 2017 at 10:20
  • $\begingroup$ Do you desire points to be "stuck" on the boundary? (That is: "is the convex hull of the solution points EQUAL to the rectangle you are trying to fill?") Or may points be off-boundary? In the latter case, does the "distance from a point near the edge to its nearest edge" matter? $\endgroup$
    – Dan H
    Commented Jun 1, 2018 at 19:54

3 Answers 3


There is a neat (and widely known) trick to produce n approximately evenly distributed points in a disk, or on any surface of revolution. Points are placed on a spiral, each time turning by the "golden angle". The reference is:

We can adapt the same thing to a rectangle by "wrapping the $x$ coordinate around".

Let r be the aspect ratio of the rectangle (height/width) and n the number of points. Then use

r = 1.5;
n = 150;

 Point@Table[{Mod[k/GoldenRatio, 1], r k/n}, {k, 1., n}]

The nice thing about this arrangement is that the points will look approximately equally spaced regardless of the aspect ratio. You can compress or stretch the rectangle along either axis and they will still look equally spaced.

equidistributed points in a rectangle

  • 2
    $\begingroup$ The "phyllotactic" arrangement of points was also featured on the Wolfram Blog. Here's a slightly simplified version: With[{n = 150, w = 2, h = 3}, Graphics[Point[N @ Table[{Mod[w k (GoldenRatio - 1), w], h k/n}, {k, n}]]]] $\endgroup$ Commented Jan 29, 2017 at 11:55
  • $\begingroup$ Extra points for the "squishing" animation... and it's true that this operation preserves the equality of the edges connecting the nodes... HOWEVER it does fly in the face of "preserving the equidistance" because the "diagonal distance" ("across" a cell) soon becomes noticeably smaller than the original "edge" length. $\endgroup$
    – Dan H
    Commented Jun 1, 2018 at 19:27

Lloyd relaxation (previously featured here) gives a "nice" point distribution:

pts = With[{n = 150, w = 2, h = 3,
            maxit = 50, (* maximum iterations *)
            tol = 0.001 (* distance tolerance *)},
                               Block[{cells, ci, vm},
                                     vm = VoronoiMesh[pts, {{0, w}, {0, h}}];
                                     cells = MeshPrimitives[vm, "Faces"];
                                     ci = Region`Mesh`MeshMemberCellIndex[vm];
                                     RegionCentroid /@ cells[[ci[pts][[All, -1]]]]]], 
                      RescalingTransform[{{0, 1}, {0, 1}}, {{0, w}, {0, h}}] @
                      RandomReal[1, {n, 2}], maxit, 
                      SameTest -> (Max[MapThread[EuclideanDistance, {#1, #2}]] < tol &)]];


Lloyd-relaxed points


My implementation of the Poisson-disc distribution from this answer perhaps fits the bill:

findBestCandidate[samplePoints_, nrOfCandidates_, {w_, h_}] :=
Module[{nearestFunction, candidates},
  nearestFunction = Nearest[samplePoints];
  candidates = Transpose[{RandomInteger[{1, w}, nrOfCandidates],
                          RandomInteger[{1, h}, nrOfCandidates]}];
  Last @ SortBy[candidates, Norm[# - First @ nearestFunction @ #] &]

sample[nrOfSamplePoints_, nrOfCandidates_, {w_, h_}] := Nest[
  # ~Append~ findBestCandidate[#, nrOfCandidates, {w, h}] &,
  {{RandomInteger[w], RandomInteger[h]}},
  nrOfSamplePoints - 1

We can use it like this:

sample[113, 100, {3000, 1000}] // Point // Graphics

Mathematica graphics

113 is the number of points, and 100 is the number of iterations allowed to find the best position (so the higher number, the better), and {3000, 1000} is the size of the grid. The size of the grid controls the resolution (also in this case, the higher the better).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.