Generate random coordinates inside a convex polytope

I am trying to generate a random set of coordinates inside a randomly-shaped (convex) polytope defined by its bounding surfaces. My idea is to generate random trial coordinates in the smallest box containing the domain, and check if the point actually lies inside the polytope after. I am struggling with this last part: how do you check if a point is inside a given polytope?

[Matt J]

My idea is to generate random trial coordinates in the smallest box containing the domain, and check if the point actually lies inside the polytope after.

Is this a 3D polytope? In higher dimensions, the minimum bounding box can be quite a bit larger in volume than the polytope. It can take a huge number of randomizations to hit a point lying in the polytope depending on the dimension of the space.

[Alessandro]

Yes, it is a 3dimensional physical domain. Usually very regular so that the bounding box is comparable in volume with the polytope itself.

[John D'Errico]

No rejection needed. See my answer.

Answer 1

[John D'Errico]

This is actually easier than you think, IF you understand how to solve the necessary sub-problems. I even have code that does it, but I tend not to give it out, because people seem not to understand that toolbox. Interface is everything. Well, not everything, but important as hell. And since I've never had the ambition to create a complete new interface for those tools, they will languish on my list of things to do.

So, lets see if we can solve the problem. Is this a 2-d problem? 3-d? N-d? I'll show you how to do it in 2-d, because I can make pretty plots there that are easy to look at and verify we were successful. The extension to n-d is trivial.

First, create a convex tessellation.

xy = randn(10,2);
CH = convhulln(xy)
CH =
   1     8
   8    10
  10     6
   6     2
   2     7
   7     1

Those are edges, because I'm working in 2-d. n-d is the same idea.

plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(CH(:,1),1),xy(CH(:,2),1)]',[xy(CH(:,1),2),xy(CH(:,2),2)]','r-')

Compute a point in the interior of the polygon. Since it is known to be a convex polygon (polyhedron in n-d), mean will suffice.

xycent = mean(xy,1)
xycent =
   0.026363      0.42644

nxy = size(xy,1)
ncent = nxy+1;
xy(ncent,:) = xycent;

now convert the polygon into a triangulation of that convex domain. Each edge will become a triangle, or (n-1)-simplex in n-d. (A k-simplex is defined by k+1 vertices. It may be full volume, or it may live in a k-manifold, embedded in a higher dimensional space.) These are now triangles, living in the 2-d data plane.

tri = [CH,repmat(ncent,ntri,1)]
tri =
   1     8    11
   8    10    11
  10     6    11
   6     2    11
   2     7    11
   7     1    11

figure
plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(tri(:,1),1),xy(tri(:,2),1),xy(tri(:,3),1)]',[xy(tri(:,1),2),xy(tri(:,2),2),xy(tri(:,3),2)]','g-')

You can see one new point in that set, the one we added. It is a common vertex for every triangle.

Next, compute the volume of each simplex. In 2-d, it is just the area. This is mot easily done using determinants, with a loop. Yes, you might know that I hate determinants in general. But sometimes they can be useful, and I'm feeling too lazy.

Note that the loop below uses aR2016b feature. Earlier releases will need to use bsxfun.

V = zeros(1,ntri);
for ii = 1:ntri
  V(ii) = abs(det(xy(tri(ii,1:2),:) - xycent));
end

V
V =
       1.7526       1.2067      0.79587      0.94023       1.0889        2.993

So V is a list of areas (or volumes) for each simplex. Actually I needed to divide by factorial(ndim), where ndim is the dimension of the space we are living in to truly compute volume, but I need only something proportional to volume for the next computation. If that bothers you here, you could divide V by 2, but save the cpu cycles, because the next step is to normalize V to sum to 1.

V = V/sum(V)
V =
    0.19968      0.13747     0.090674      0.10712      0.12406      0.34099

All we need are the relative areas of each simplex.

We are almost done. Lets say we wish to generate M points in total.

M = 1000;

For each point, we first need to know which simplex it lives in. That will be proportional to the relative volume of the corresponding simplexes.

[~,~,simpind] = histcounts(rand(M,1),cumsum([0,V]));

So the bigger simplexes will have proportionally more points to fall in them.

FINALLY, generate the desired random points. Yes, I know this seems like a trip to get here, but most of it was me typing explanation and making pictures.

Again, the code that follows will require R2016b. Easy to fix though.

r1 = rand(M,1);
uv = xy(tri(simpind,1),:).*r1 + xy(tri(simpind,2),:).*(1-r1);
r2 = sqrt(rand(M,1));
uv = uv.*r2 + xy(tri(simpind,3),:).*(1-r2);

uv is the complete set of random points. Trust me? Yeah, I know, trust but verify.

figure
plot(xy(:,1),xy(:,2),'bo');
hold on
plot([xy(tri(:,1),1),xy(tri(:,2),1),xy(tri(:,3),1)]',[xy(tri(:,1),2),xy(tri(:,2),2),xy(tri(:,3),2)]','g-')
plot(uv(:,1),uv(:,2),'m.')

As you can see, the sampling is indeed uniform in that domain.

In 3 or more dimensions, just extend what I did as is appropriate. Ask if you have questions, but the extension really is easy. (Ok, there is one part that requires some thought by the programmer, and that is how to extend the very last sampling step into n-d. The trick is to understand why I used sqrt there.) Again, ask if you need the trick.

I suppose one day I should post code on the FEX that does this, but sampling from a convex n-dimensional polyhedron is not something for which I see people clamoring.

[Matt J]

Well, I like it, but I don't think the role of the square root and how to extend to 3D and higher is at all obvious.

An alternative that I've seen is to randomize over a parallelogram (or parallelepiped in 3D)

pt= rand*x+rand*y;  %2D
pt= rand*x+rand*y+rand*z; %3D

where x,y,z are rays defining the appropriate simplex. Then, you cut that region in half to obtain points in a triangle or tetrahedron. You either discard points in the wrong half or reflect them across the cutting plane to get them in correct region.

[John D'Errico]

I did say that the sqrt is the point of interest. A cube root will be of interest too. :)

So lets see what happens. Suppose you want to generate a random point inside an n-simplex? Lets see how you would do it for a triangle. A simple isosceles triangle is a great place to start, because it is easy to visualize what happens. So, consider the triangle with vertices at {(0,0), (0,1), (1,0)}. We want a uniformly random point in the triangle, so any point in the triangle will be equally probable. Start by slicing the triangle into infinitesimally thin slices, parallel to the y axis. Since we said the sampling will be uniform in the triangle, what is the marginal PDF, found by integrating over y? In fact, we would expect to see a triangular marginal PDF on x. That is, if P(x,y) is the PDF over this isosceles triangle, then if we integrate over y for any value of x, we get a marginal PDF.

Over the triangle:

P(x,y) = 2

Of course that 2-d PDF integrates to 1 over the area of the triangle.

The marginal PDF:

P_x(x) = 2*(1-x)

We can see that this marginal PDF also integrates to 1.

syms x
int(2*(1-x),[0 1])
ans =
1

Ok, so why have I bothered to go through all of this?

We will generate a random point in the triangle, by first generating a random point on AN edge of the triangle. Start with a random point along the edge between points {(0,0), (0,1)}. This is easy to do.

r1 = rand(M,1);
uv = [0 0].*r1 + [0 1].*(1-r1);

So the rows of uv are randomly sampled along the left edge of the triangle.

Now, generate a random point along the line segments between each of those points, and the third vertex of the triangle, at (1,0).

The distance along that line segment will be determined by the marginal PDF we generated before! So points near the left edge of the triangle will have higher probability of occurrence, that probability will decrease with the marginal PDF, so in this case, it will decrease linearly with distance between the edge and the third vertex.

So, if the marginal PDF in x is P_X(x)=2-2*x, then how do we sample? You compute the CDF.

CDF_X(x) = int(2*(1-u),[0,x])
ans =
-x*(x - 2)

So we have

CDF_x = 2*x - 2*x^2

Sampling from this marginal distribution is most easily done by inverting the CDF. Thus, generate a uniform random point (r2) in the interval [0,1]. Then invert the relation

CDF_X(x) = r2

You should be seeing where the sqrt arises now. Or maybe I've been too slow in my explanation and you saw it long ago. As it turns out, all you need to do is

r2 = sqrt(rand(M,1));
uv = uv.*r2 + [1 0].*(1-r2);

So lets try it out.

M = 1000;
r1 = rand(M,1);
uv = [0 0].*r1 + [0 1].*(1-r1);
r2 = sqrt(rand(M,1));
uv = uv.*r2 + [1 0].*(1-r2);
plot(uv(:,1),uv(:,2),'.')

Did it work? A picture is worth a thousand words. Wow! I must have gotten lucky. :)

As it turns out, the exact same procedure works to sample from ANY arbitrary triangle.

1. Generate points randomly along ONE edge. 2. Extend a segment from those points to the third vertex. Sample along that edge using the sqrt trick I showed.

Ok, so how do you extend this procedure to 3 dimensions? I did say it was easy, as long as you understand the sqrt, and where it came from.

A tetrahedron in 3-d has triangular facets. So, sample using the above scheme UNIFORMLY on that facet. Then sample along the line segment connecting those point to the 4TH vertex of the simplex. Use a cube root to sample along that line segment.

Lets try it, for a unit isosceles tetrahedron in 3-dimensions.

M = 10000;
r1 = rand(M,1);
uvw = [0 0 0].*r1 + [0 0 1].*(1-r1);
r2 = sqrt(rand(M,1));
uvw = uvw.*r2 + [0 1 0].*(1-r2);
r3 = nthroot(rand(M,1),3);
uvw = uvw.*r3 + [1 0 0].*(1-r3);

plot3(uvw(:,1),uvw(:,2),uvw(:,3),'.')
grid on
view(50,30)

Yes, those points are uniformly sampled over the simplex in question. You can convince yourself of that if you rotate the plot around in 3-d.

Easy, peasy. :) Ok, the cube root is not obvious, at least to someone who has never done this before. There is always a first time for everything. And now, easy peasy to all. In n-dimensions, you just extend the same trick to 4, 5, and higher dimensions. The beauty of this it is all vectorized, and requires no more than a loop over the dimensions of the space you are living in.

[fatima-zahra achbah]

Hi, can you please explain to me how can we do it the same to a parallelepiped. I am truly stuck

[Matt J]

You should be able to use the solution above unchanged. None of it was tied to a particular type of polygon.

[Carlos Ricci]

Thank you for this contribution. The explanation and development of the background were clear to follow (and thank you, not too slow!). It will enable the solving of an otherwise-vexing Monte-Carlo-style error-minimization problem. While your inhull() test is also excellent (and equally amazing), in high dimensions the test/reject method begins to take on astronomical time due to the small volume of the N-d convex hull relative to the bounding box (or even ellipsoid). This volume-generative method between the facet simplices and central point is nothing short of pure genius.

[Terrance Nearey]

@John D'Errico:

Here is one person clamoring (politely) for an N-d version on FEX.

[Mohammadali Kargar]

@John D'Errico

It was so genius! tnx

[John D'Errico]

Random points inside a polyhedron in n-dimensions? That is not that difficult, depending on how you define the polyhedron. Hmm. I can probably post that with not too much effort. Sorry I did not see the comment by Terrance sooner.

Answer 2

[Matt J]

A polytope is described by linear in/equalities

   A*x<=b
   Aeq*x=beq

To see if a point lies in the polytope, test whether it satisfies the inequalities to some tolerance.

You should be able to construct the system of inequalities directly, since you say you already know the bounding surfaces. If, on the other hand, you are starting with the polytope vertices, you can derive the inequalities, for example, using VERT2LCON ( Download ).

[Alessandro]

Thanks. I am trying to understand what your function does. The idea is that all the vertices of my polytope (n,3 size) will give me A,b,Aeq,beq and if my trial coordinate x satisfies both (in)equalities the point lies inside the polytope?

[Matt J]

Yes, that is the idea. Note, however, that you can test not just one, but many points using vectorized linear algebra operations, i.e.,

 isinpolytope =  all(  bsxfun(@le, A*X,b) ,1)

where X is a 3xN matrix of trial points.

Answer 3

[Walter Roberson]

https://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume-