Matlab, generate and plot a point cloud distributed within a triangle

Question

I'm trying to generate a cloud of 2D points (uniformly) distributed within a triangle. So far, I've achieved the following:

The code I've used is this:

N = 1000;
X = -10:0.1:10;
for i=1:N
    j = ceil(rand() * length(X));
    x_i = X(j);
    y_i = (10 - abs(x_i)) * rand;

    E(:, i) = [x_i y_i];
end

However, the points are not uniformly distributed, as clearly seen in the left and right corners. How can I improve that result? I've been trying to search for the different shapes too, with no luck.

Defroster

Answer 1

You should first ask yourself what would make the points within a triangle distributed uniformly.

To make a long story short, given all three vertices of the triangle, you need to transform two uniformly distributed random values like so:

N = 1000;                    % # Number of points
V = [-10, 0; 0, 10; 10, 0];  % # Triangle vertices, pairs of (x, y)
t = sqrt(rand(N, 1));
s = rand(N, 1);
P = (1 - t) * V(1, :) + bsxfun(@times, ((1 - s) * V(2, :) + s * V(3, :)), t);

This will produce a set of points which are uniformly distributed inside the specified triangle:

scatter(P(:, 1), P(:, 2), '.')

Note that this solution does not involve repeated conditional manipulation of random numbers, so it cannot potentially fall into an endless loop.

For further reading, have a look at this article.

Noachian

Can you please explain how you calculated P? or just where can I read about the formula you used? Thanks! – Gingili

@CyberLingo Basically, it randomizes a point (x,y), represented by t and s in our case, and remaps it into a space confined by a triangle (represented by V). P is the product of the mapping matrix and (x, y). Note that the sqrt in the randomized t is just a way of simplifying the final calculation of the product. However, I think it would be easier to point you to a few links that explain it much better, such as this question, or this link... – Noachian

@CyberLingo Oh, and also note that this calculation randomizes an array of points (therefore yielding an array of points). Hence the use of bsxfun in the calculation of P. – Noachian

Answer 2

That concentration of points would be expected from the way you are building the points. Your points are equally distributed along the X axis. At the extremes of the triangle there is approximately the same amount of points present at the center of the triangle, but they are distributed along a much smaller region.

The first and best approach I can think of: brute force. Distribute the points equally around a bigger region, and then delete the ones that are outside the region you are interested in.

N = 1000;
points = zeros(N,2);
n = 0;
while (n < N)
    n = n + 1;
    x_i = 20*rand-10; % generate a number between -10 and 10
    y_i = 10*rand; % generate a number between 0 and 10
    if (y_i > 10 - abs(x_i)) % if the points are outside the triangle
       n = n - 1; % decrease the counter to try to generate one more point
    else % if the point is inside the triangle
       points(n,:) = [x_i y_i]; % add it to a list of points
    end
end

% plot the points generated
plot(points(:,1), points(:,2), '.');
title ('1000 points randomly distributed inside a triangle');

The result of the code I've posted:

one important disclaimer: Randomly distributed does not mean "uniformly" distributed! If you generate data randomly from an Uniform Distribution, that does not mean that it will be "evenly distributed" along the triangle. You will see, in fact, some clusters of points.

Natelson

Manually altering the iteration variable of a for loop in MATLAB is not a good idea. I don't think that this behaves like you're expecting it to. – Noachian

@Natelson - In matlab you cannot change the loop variable inside a for loop. Try the following for ii=1:2, ii=1, end; this loop will be executed only twice and will not result with an infinite loop. – Ezechiel

the code above works at Matlab R2011a, as proven by the plot. test it ;) – Natelson

Yes, the loop can be changed easily to use a while loop. For was the first approach because it is easier to write. From my experience, you can change the loop variable inside it, provided there is no nested loop. But I'll change the code to a while loop as it is indeed better. (p.s.: the code works for 5 points, 1000 points and 100000 points) – Natelson

@Natelson - (regarding the for version): when you say that "the code works" - are you sure it produces exactly N points? It seems like quite a few points falls at (0,0) due to the behavior of Matlab's for iterator. In your plot there is a point at (0,0) which I suspect is in fact quite a few points plotted one over the other... Please see my comment regarding for ii=1:2, ii=1, end; – Ezechiel

You are right, the loop does not get repeated. Thanks for pointing that out. Lesson learned :) blogs.mathworks.com/loren/2006/07/19/how-for-works – Natelson

Answer 3

You can imagine that the triangle is split vertically into two halves, and move one half so that together with the other it makes a rectangle. Now you sample uniformly in the rectangle, which is easy, and then move the half triangle back.

Also, it's easier to work with unit lengths (the rectangle becomes a square) and then stretch the triangle to the desired dimensions.

x = [-10 10]; % //triangle base
y = [0 10]; % //triangle height
N = 1000; %// number of points

points = rand(N,2); %// sample uniformly in unit square
ind = points(:,2)>points(:,1); %// points to be unfolded 
points(ind,:) = [2-points(ind,2) points(ind,1)]; %// unfold them
points(:,1) = x(1) + (x(2)-x(1))/2*points(:,1); %// stretch x as needed
points(:,2) = y(1) + (y(2)-y(1))*points(:,2); %// stretch y as needed
plot(points(:,1),points(:,2),'.')

Pt

Answer 4

We can generalize this case. If you want to sample points from some (n - 1)-dimensional simplex in Euclidean space UNIFORMLY (not necessarily a triangle - it can be any convex polytope), just sample a vector from a symmetric n-dimensional Dirichlet distribution with parameter 1 - these are the convex (or barycentric) coordinates relative to the vertices of the polytope.

Filipino