How to create triangular lattice from two points?

Hey,

Let us say you have an m-by-n matrix, and in this matrix there are two points, A(x,y) and B(x,y). Can anyone help me figure out how to construct a triangular lattice of equilateral triangles where A and B are two of the vertices?

That is, if the line between A and B is one side of an equilatereal triangle, then point C will be located at equal distance from A and B, forming angles of 60 degrees. A point, D, will again form a triangle with points C and B and so on (see image).I want a code to find all vertices of such a triangular lattice within the boundaries of the matrix.

Any help is appreciated.

Thanks

Answer 1

[Guillaume]

I'm going to assume that given A and B you know how to find C (and C' the symmetric of C against AB). If not, use your favorite search engine, it's simple math.

Now, to fill your matrix, which I'll call area instead as it's not a matrix:

areabounds = [0 60;    ... x
              0 60];   %   y

You start with a list of points, an Nx2 matrix, which to starts with has only A and B:

points = [xa ya; 
          xb yb];  %will grow

and you start with a list of edges to process, a Mx2 matrix of point indices. To start with, it's just AB:

edges = [1 2]; %point 1 to point 2. will grow/shrink

Then your algorithm is simple. Take the first edge. Find the two points forming the triangles and if not present in the point list and within your area bounds, add them and add the 2 new edges from each point to the edge list. Remove the edge you've processed from the list and move on the next until the list is empty:

while ~isempty(edges)
   vertices = getvertices(points(edges(1, 1), :), points(edges(1, 2), :));  %your function for finding C and C'
   %vertices is a 2x2 matrix. row 1 is [x, y] of C, row2 is [x, y] of C'
   isnew = ismembertol(vertices, points, 'ByRows', true);  %are vertices new? Using ismembertol to avoid floating point errors. 
   isinbound = all(vertices(:, 1) >= areabounds(:, 1).', 2) & all(vertices(:, 2) <= areabounds(:, 2).', 2); %require R2016b or later. Use bsxfun in earlier versions
   for newpoint = vertices(isnew & isinbound).'  %transpose since for iterates over columns
     points(end+1, :) = newpoint.';  %add new point to point list
     edges(end+1:end+2, :) = [size(points, 1), edges(1, 1);  ... add CA to edge list
                              size(points, 1), edges(1, 2)]; %   add CB to edge list
   end
   edges = edges(2:end, :);  %remove processed edge
end  %move on to next edge

[John D'Errico]

+1. That is exactly how I would attack the problem, growing the triangulation from one edge, being careful not to create the same triangle twice. Work only on edges that are on the boundary of the triangulation. Each time you create a new triangle, you add two new edges to the list of those on the boundary. If one (or both) of those new edges were already on the boundary list, they you need to treat them properly. Use ismembertol to avoid duplication of vertices. If a new triangle would wander outside of the boundary, then do not accept it, and remove the edge that created it from the list of edges for future growth.

Basically a while loop, an exercise in careful bookkeeping.

Be careful in that the arrays generated will grow, so the whole mess will slow down quadratically with time. You can avoid this by preallocation, since you know how many triangles to expect in the end! That is, you know the final area of the domain. You know the area of one triangle. So you can make a very good estimate of the total number of triangles, and the total number of edges. Use an over-estimate to preallocate your arrays, then delete the unneeded space at the end.

[OH]

Thanks for the help! It seems good except that the ismembertol function does not see the vertices as new. With areabounds = [0 60;0 60]; points = [30 30;40 40]; and the output from getvertices: vertices = [26.3 43.6; 43.6 26.3], I get logical false from the ismembertol on the first iteration of the while loop. Might be something with the tolerance?

[John D'Errico]

It has nothing to do with the tolerance. Clearly these two points are not the same.

[26.3 43.6; 43.6 26.3]

They are mirror images across a line of symmetry.

It is true however that I would use a larger tolerance than the default. This is completely safe, since if properly generated, the new points will fall on a lattice. Any deviation from that will be in the least significant bits. But remember that the floating point trash will fluctuate as essentially a random walk in two dimensions. So the variance of the LSB trash will gradually grow in a predictable manner.

[Guillaume]

Well, it's untested code. Expect bugs! The result of ismembertol is correct, when a point is new it is not a member. I simply forgot to invert the result:

isnew = ~ismembertol(vertices, points, 'ByRows', true);

Unless your triangles are very small with regards to the size of your area, the default ismembertol tolerance should be fine.

[OH]

That works! Thanks!