Random points inside a parallelogram

Question

I have a 4 side convex Polygon defined by 4 points in 2D, and I want to be able to generate random points inside it.

If it really simplifies the problem, I can limit the polygon to a parallelogram, but a more general answer is preferred.

Generating random points until one is inside the polygon wouldn't work because it's really unpredictable the time it takes.

Answer 1

The question by the OP is a bit ambiguous so the question I will answer is: How to generate a point from a uniform distribution within an arbitrary quadrilateral, which is actually a generalization of How to generate a point from a uniform distribution within an arbitrary (convex) polygon. The answer is based on the case of generating a sample from a uniform distribution in a triangle (see http://mathworld.wolfram.com/TrianglePointPicking.html, which has a very nice explanation).

In order to accomplish this we:

1. Triangulate the polygon (i.e. generate a collection of non-overlapping triangular regions which cover the polygon). For the case of a quadrilateral, create an edge across any two non-adjacent vertices. For other polygons, see http://en.wikipedia.org/wiki/Polygon_triangulation for a starting point, or http://www.cgal.org/ if you just need a library.

   

2. To pick one of the triangles randomly, let us assign an index to each triangle (i.e. 0,1,2,...). For the quadrilateral, they will be 0,1. For each triangle we assign a weight equal as follows:

   

3. Then generate a random index i from the finite distribution over indexes given their weights. For the quadrilateral, this is a Bernoulli distribution:

   

4. Let v0, v1, v2 be vertices of the triangle (represented by their point locations, so that v0 = (x0,y0), etc. Then we generate two random numbers a0 and a1, both drawn uniformly from the interval [0,1]. Then we calculate the random point x by x = a0 (v1-v0) + a1 (v2-v0).

   

5. Note that with probability 0.5, x lies outside outside the triangle, however if it does, it lies inside the parallelogram composed of the union of the triangle with it's image after a rotation of pi around the midpoint of (v1,v2) (dashed lines in the image). In that case, we can generate a new point x' = v0 + R(pi)(x - v3), where R(pi) is a rotation by pi (180 deg). The point x' will be inside the triangle.

6. Further note that, if the quadrilateral was already a parallelogram, then we do not have to pick a triangle at random, we can pick either one deterministically, and then choose the point x without testing that it is inside it's source triangle.

Answer 2

A. If you can restrict your input to parallelogram, this is really simple:

1. Take two random numbers between 0 and 1. We'll call then u and v.

2. If your parallelogram is defined by the points ABCD such that AB, BC, CD and DA are the sides, then take your point as being:

    p = A + (u * AB) + (v * AD)
   

Where AB is the vector from A to B and AD the vector from A to D.

B. Now, if you cannot, you can still use the barycentric coordinates. The barycentric coordinates correspond, for a quad, to 4 coordinates (a,b,c,d) such that a+b+c+d=1. Then, any point P within the quad can be described by a 4-uple such that:

P = a A + b B + c C + d D

In your case, you can draw 4 random numbers and normalize them so that they add up to 1. That will give you a point. Note that the distribution of points will NOT be uniform in that case.

C. You can also, as proposed elsewhere, decompose the quad into two triangles and use the half-parallelogram method (i.e., as the parallelogram but you add the condition u+v=1) or the barycentric coordinates for triangles. However, if you want uniform distribution, the probability of having a point in one of the triangle must be equal to the area of the triangle divided by the area of the quad.

Answer 3

Assuming you want a uniform distribution: Form two triangles from your polygon. Pick which triangle to generate the point in according to their area ratio.

Call the corners of the triangle A, B, C, the side vectors AB, BC, AC and generate two random numbers in [0,1] called u and v. Let p = u * AB + v * AC.

If A+p is inside the triangle, return A+p

If A+p is outside the triangle, return A + AB + AC - p

(This is basically PierreBdR's formula except for the preprocessing and the last step that folds the point back into a triangle, so it can handle other shapes than parallelograms).

Answer 4

Your polygon is two triangles, so why not randomly select one of those, then find a random point in the triangle.

Probably not the best solution, but it'd work.

Answer 5

A somewhat less "naïve" approach would be to use a polygon fill algorithm, and then select points from the fill lines randomly.

C Code Sample

//  public-domain code by Darel Rex Finley, 2007

int  nodes, nodeX[MAX_POLY_CORNERS], pixelX, pixelY, i, j, swap ;

//  Loop through the rows of the image.
for (pixelY=IMAGE_TOP; pixelY<IMAGE_BOT; pixelY++) {

  //  Build a list of nodes.
  nodes=0; j=polyCorners-1;
  for (i=0; i<polyCorners; i++) {
    if (polyY[i]<(double) pixelY && polyY[j]>=(double) pixelY
    ||  polyY[j]<(double) pixelY && polyY[i]>=(double) pixelY) {
      nodeX[nodes++]=(int) (polyX[i]+(pixelY-polyY[i])/(polyY[j]-polyY[i])
      *(polyX[j]-polyX[i])); }
    j=i; }

  //  Sort the nodes, via a simple “Bubble” sort.
  i=0;
  while (i<nodes-1) {
    if (nodeX[i]>nodeX[i+1]) {
      swap=nodeX[i]; nodeX[i]=nodeX[i+1]; nodeX[i+1]=swap; if (i) i--; }
    else {
      i++; }}

  //  Fill the pixels between node pairs.
  //  Code modified by SoloBold 27 Oct 2008
  //  The flagPixel method below will flag a pixel as a possible choice.
  for (i=0; i<nodes; i+=2) {
    if   (nodeX[i  ]>=IMAGE_RIGHT) break;
    if   (nodeX[i+1]> IMAGE_LEFT ) {
      if (nodeX[i  ]< IMAGE_LEFT ) nodeX[i  ]=IMAGE_LEFT ;
      if (nodeX[i+1]> IMAGE_RIGHT) nodeX[i+1]=IMAGE_RIGHT;
      for (j=nodeX[i]; j<nodeX[i+1]; j++) flagPixel(j,pixelY); }}}

   // TODO pick a flagged pixel randomly and fill it, then remove it from the list.
   // Repeat until no flagged pixels remain.

Answer 6

By "general" do you mean all non-parallelogram 4-side polygons in general or all possible polygons?

How about drawing a random line connecting the 4 sides e.g. If you have this:

.BBBB.
A    C
A    C
.DDDD.

Then generate a random point on a unit square, then mark the point on the line B and D at the percentage of distance on the X axis. Do the same on line A and C using value from the Y axis.

Then connect the point on line A to line C and line B to line D, the intersection point is then used as the random point.

It's not uniform because rounding errors will aid certain points but it should be close if you are working with floating points values.

Implementation should be rather easy, too, since you are already working with polygons. You should already have code that does those simple tasks.

Here's a quick pseudocode:

void GetRandomPoint(Polygon p, ref float x, ref float y) {

    float xrand = random();
    float yrand = random();

    float h0 = p.Vertices[0] + xrand * p.Vertices[1];
    float h1 = p.Vertices[2] + yrand * p.Vertices[3];

    float v0 = p.Vertices[0] + xrand * p.Vertices[2];
    float v1 = p.Vertices[1] + yrand * p.Vertices[3];

    GetLineIntersection(h0, h1, v0, v1, x, y);

}

Answer 7

This works for general, convex quadrilaterals:

You can borrow some concepts from the Finite Element Method, specifically for quadrilateral (4-sided) elements (refer to section 16.5 here). Basically, there is a bilinear parameterization that maps a square in u-v space (for u, v \in [-1, 1] in this case) to your quadrilateral that consists of points p_i (for i = 1,2,3,4). Note that In the provided reference, the parameters are called \eta and \xi.

Basic recipe:

1. Choose a suitable random number generator to generate well-distributed points in a square 2D domain
2. Generate random u-v pairs in the range [-1, 1]
3. For each u-v pair, the corresponding random point in your quad = 1/4 * ((1-u)(1-v) * p_1 + (1+u)(1-v) * p_2 + (1+u)(1+v) * p_3 + (1-u)(1+v) * p_4)

The only problem is that uniformly distributed points in the u-v space won't produce uniformly distributed points in your quad (in the Euclidean sense). If that is important, you can work directly in 2D within the bounding box of the quad and write a point-in-quad (maybe by splitting the problem into two point in tris) test to cull random points that are outside.

Answer 8

Do the points need to be uniformly distributed, or is any distribution ok?

Can the polygon be concave, or is it guarenteed to be convex?

If the answer to both the above is no, then pick any two of the vertexes and pick a random point on the line segment between them. This is limited to the line segements connecting the vertexes (ie, VERY non-uniform); you can do a bit better by picking a third vertex and then picking a point between that and the first point -- still non-uniform, but at least any point in the polygon is possible

Picking a random point on a line between two points is easy, just A + p(B-A), where A and B are the points and p is a random number between 0.0 and 1.0

Answer 9

What kind of distribution do you want the points to have? If you don't care, the above methods will work fine. If you want a uniform distribution, the following procedure will work: Divide the polygon into two triangles, a and b. Let A(a) and A(b) be their areas. Sample a point p from the uniform distribution on the interval between 0 and A(a)+A(b). If p < A(a), choose triangle a. Otherwise, choose triangle b. Choose a vertex v of the chosen triangle, and let c and d be the vectors corresponding to the sides of the triangle. Sample two numbers x and y from the exponential distribution with unit average. Then the point (xc+yd)/(x+y) is a sample from the uniform distribution on the polygon.

Answer 10

The MATLAB function cprnd generates points from the uniform distribution on a general convex polytope. For your question a more specialized algorithm based on decomposing the quadrilateral into triangles is more efficient.

Answer 11

For PostGIS, this is what I am using (you might want a ward for possible infinite loops). You might export the algorithm to your programming language:

CREATE or replace FUNCTION random_point(geometry)
RETURNS geometry
AS $$
DECLARE 
    env geometry;
    corner1 geometry;
    corner2 geometry;
    minx real;
    miny real;
    maxx real;
    maxy real;
    x real;
    y real;
    ret geometry;
begin

select ST_Envelope($1) into env;
select ST_PointN(ST_ExteriorRing(env),1) into corner1;
select ST_PointN(ST_ExteriorRing(env),3) into corner2;
select st_x(corner1) into minx;
select st_x(corner2) into maxx;
select st_y(corner1) into miny;
select st_y(corner2) into maxy;
loop
    select minx+random()*(maxx-minx) into x;
    select miny+random()*(maxy-miny) into y;
    select ST_SetSRID(st_point(x,y), st_srid($1)) into ret;
    if ST_Contains($1,ret) then
        return ret ;
    end if;
end loop;
end;
$$
LANGUAGE plpgsql
volatile
RETURNS NULL ON NULL INPUT;