MATLAB: Generate random coordinates inside a convex polytope

polytoperandom number generatorrandomize over nonrectangular region

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?

Best Answer

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.

Related Solutions

MATLAB: Pretty simple question regarding convhulln. Now that I have K…

(Somebody should post a thing on the FEX explaining all the basics about convex hulls, tessellations, etc. One more thing to add to the round tuit list.)

Think in 2-d first. A convex hull of some data points is a simple polygon, composed of linear segments. Any vertex of the convex hull is also one of the data points. A convex hull does not introduce new points. So essentially a convex hull is just a set of references into the original set of points. Since these things are best done visually, lets generate an example.

>> xy = rand(8,2)
xy =
      0.80028      0.84913
      0.14189      0.93399
      0.42176      0.67874
      0.91574      0.75774
      0.79221      0.74313
      0.95949      0.39223
      0.65574      0.65548
     0.035712      0.17119
>> plot(xy(:,1),xy(:,2),'o')

See that there appear to be 5 of those points that lie on the boundary of the convex hull in this case, and 3 points lie clearly in the interior.

We can find the hull using either convhull, or convhulln.

>> convhull(xy(:,1),xy(:,2))
ans =
     1
     2
     8
     6
     4
     1
>> convhulln(xy)
ans =
     8     6
     6     4
     4     1
     1     2
     2     8

What do these numbers mean? Each index returned corresponds to ONE vertex of the hull polygon, but also to one of our original data points. Thus, the convex hull goes between the points [1 2 8 6 4 1] in the set IN THAT ORDER. The points that are inside the hull are of course then numbered [3 5 7].

See that each row of the array returned by convhulln is the same information, just in a slightly different form. It tells us that one edge of the hull is the edge [8 6], connecting point 8 and point 6. Then we connect points 6 and 4, with another edge, etc. These edges in the array returned from convhulln are given in no specific order, but they will be the same information overall as that returned by convhull.

The difference is that convhulln will work in any number of dimensions. (Ok, it can get pretty nasty to visualize in high dimensions. But 3 or 4 or 5 dimensions is entirely adequate. And we can't even draw pictures past 3-d anyway.)

So think about the 2-d case. The convex hull in 2-d is composed of "facets" composed of edges, line segments. In 3-d, the convex hull will always be composed of triangles. A 3-d convex hull will have planar facets. They will always break down into triangular facets though. Yes, one could imagine the set of points that comprise the unit 3-cube. (A cube in R^3.) Clearly its convex hull is a set of squares.

So here are the 8 vertices of the unit 3-cube.

V =
     0     0     0
     1     0     0
     0     1     0
     1     1     0
     0     0     1
     1     0     1
     0     1     1
     1     1     1

ndgrid could essentially help you to generate that list of points if you wanted.

Let us generate the convex hull using matlab. It is a set of triangles.

>> convhulln(V)
ans =
     4     3     1
     2     4     1
     6     2     1
     5     6     1
     6     4     2
     4     6     8
     3     7     1
     7     5     1
     4     7     3
     7     4     8
     7     6     5
     6     7     8

Each square has been split, into a pair of triangles. Since there are 6 square faces of a cube, there will be 12 triangles.

>> size(ans)
ans =
    12     3

If you prefer to visualize them this way, we can sort the triangles in increasing order of the vertices, and then use sortrows. This does not change the information in any serious way.

>> sortrows(sort(convhulln(V),2))
ans =
     1     2     4
     1     2     6
     1     3     4
     1     3     7
     1     5     6
     1     5     7
     2     4     6
     3     4     7
     4     6     8
     4     7     8
     5     6     7
     6     7     8

Still the same information. Each row is ONE triangle. A set of references into the original set of data points. In the case of the unit 3-cube of course, EVERY point lies one the surface, in at least a few of the triangles of that hull. In fact, depending on how those surface squares are split into triangles each point may lie in from 3 to 6 different triangles.

What you need to see is that a convex hull in 3-d is composed of triangles. If we go to higher dimensions, a convex hull in n-dimensions will be composed of (n-1)-simplexes. (A k-simplex is composed of k+1 points. So an (n-1)-simplex is defined by n vertices.) For example...

>> size(convhulln(rand(100,4)))
ans =
   290     4
>> size(convhulln(rand(100,5)))
ans =
        1068           5
>> size(convhulln(rand(100,6)))
ans =
        4047           6
>> size(convhulln(rand(100,7)))
ans =
       16094           7

See that each row of the result is a set of indexes into the original array of points, and that in n-dimensions, each row will be a facet, an (n-1) simplex. We can't plot these things on a 2-d monitor or piece of paper, but who cares?

MATLAB: Generate Random Number inside a closed area

Not hard. But it takes a little effort. (I should post a tool for this on the FEX. SIGH. Actually, I have a tool that does this, but it is more complicated to use than I want it to be for me to post it.)

The general idea is to...

Triangulate the region. If you have a convex region, bounded by a polygon, then there are many ways to do this.
Determine the area of each triangle.
Choose a random triangle inside the region, with probability based on the relative fraction of the area of the differfent triangles.
Once you have the triangle, then generate a random point that lies uniformly inside the triangle.

Each of the above steps is easy enough in theory, though many users might find them complex. I can't say. In fact, all of the above steps are doable in a way that is fast, efficient, and even vectorized.

So, given a list of points that form a polygon, if the polygon is convex, a simple way to triangulate it is a delaunay triangulation. I'll want a list of vertices for later, so do it like this:

XY = [x(:),y(:)];
tri = delaunay(XY);

Or, you could do it using the delaunayTriangulation function. Or, you can do it by starting with the convex hull of the points, and then add on point at the centroid of the polygon. Then connect each edge of the convex hull to the centroidal point. Or, given a polygon, you could use ear clipping to create a triangulation. As I said, lots of ways to do it.

Next, we can compute the area of each triangle simply enough. In fact, this is easily done as a vectorized computation.

% use a 2x2 determinant, in a vectorized form
v1 = XY(tri(:,1),:);
v2 = XY(tri(:,2),:);
v3 = XY(tri(:,3),:);
    
% translate
v1 = v1-v3;
v2 = v2-v3;
    
% vectorized determinant
% divide by factorial(2) for the area
areas = (v1(:,1).*v2(:,2) - v1(:,2).*v2(:,1))/2;
% normalize the areas to sum to 1
areas = areas/sum(areas);

Now, pick a random triangle, proportional to the relative area of each triangle. You can do this for many points at once, using the discretize tool. I'll assume you want to generate nsample points. I'll do these computations for the data provided, here, for 1000 points.

nsample = 1000;
R = rand(nsample,1);
tind = discretize(R,cumsum([0;areas]));

Next, for each triangle chosen, we will find ONE point that lives inside that triangle. Again, a vectorized computation.

v1 = XY(tri(tind,1),:);
v2 = XY(tri(tind,2),:);
v3 = XY(tri(tind,3),:);
R1 = rand(nsample,1);
xyrand = v1.*repmat(R1,[1 2]) + v2.*repmat(1-R1,[1 2]);
R2 = sqrt(rand(nsample,1));
xyrand = xyrand.*repmat(R2,[1 2]) + v3.*repmat(1-R2,[1 2]);

I used repmat above, but it is also easily done using scalar dimension expansion as found in R2016b or later, or using bsxfun.

Did it work?

plot(xyrand(:,1),xyrand(:,2),'.')
hold on
plot(x,y,'r-')

Of course. Uniformly generated samples in a polygonal region. Fully vectorized too.

Best Answer

Related Solutions

MATLAB: Pretty simple question regarding convhulln. Now that I have K…

MATLAB: Generate Random Number inside a closed area

Related Question