6

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

Matlab plot

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.

Eitan T
  • 32,660
  • 14
  • 72
  • 109
José Tomás Tocino
  • 9,873
  • 5
  • 44
  • 78

4 Answers4

9

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), '.')

enter image description here

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.

Community
  • 1
  • 1
Eitan T
  • 32,660
  • 14
  • 72
  • 109
  • Can you please explain how you calculated P? or just where can I read about the formula you used? Thanks! – Ayushi Jha Mar 23 '17 at 06:35
  • 1
    @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](http://stackoverflow.com/questions/4778147/sample-random-point-in-triangle), or [this link](http://mathworld.wolfram.com/TrianglePointPicking.html)... – Eitan T Mar 23 '17 at 10:34
  • 1
    @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`. – Eitan T Mar 23 '17 at 10:38
  • 1
    Thank you for the explanation! – Ayushi Jha Mar 23 '17 at 11:38
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: Plot generated by the code above

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.

Castilho
  • 3,147
  • 16
  • 15
  • 1
    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. – Eitan T Dec 24 '12 at 12:56
  • @Castilho - 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. – Shai Dec 24 '12 at 12:57
  • the code above works at Matlab R2011a, as proven by the plot. test it ;) – Castilho Dec 24 '12 at 12:59
  • 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) – Castilho Dec 24 '12 at 13:05
  • @Castilho - (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;` – Shai Dec 24 '12 at 13:13
  • 1
    You are right, the loop does not get repeated. Thanks for pointing that out. Lesson learned :) http://blogs.mathworks.com/loren/2006/07/19/how-for-works/ – Castilho Dec 24 '12 at 13:20
1

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),'.')

enter image description here

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
0

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.

Nik Ved
  • 168
  • 2
  • 8