# Generate random coordinates inside a convex polytope

105 views (last 30 days)
Alessandro on 3 Mar 2017
Commented: Terrance Nearey on 23 Jan 2020
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 on 3 Mar 2017
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 on 3 Mar 2017
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 on 4 Mar 2017
No rejection needed. See my answer.

John D'Errico on 4 Mar 2017
Edited: John D'Errico on 7 Mar 2017
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 on 6 Oct 2017
You should be able to use the solution above unchanged. None of it was tied to a particular type of polygon.
Carlos Ricci on 2 Mar 2019
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 on 23 Jan 2020
@John D'Errico:
Here is one person clamoring (politely) for an N-d version on FEX.

Matt J on 3 Mar 2017
Edited: Matt J on 3 Mar 2017
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 on 3 Mar 2017
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 on 4 Mar 2017
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.