MATLAB: Delaunay Triangulation only for edges below certain length

delaunaymaximum edge lengthtriangulation

Hello, having tried to find out if that works for quite some time, I thought I'll try here.

So basically I am performing this triangulation on a network and at the edges of the image, there are points connected by the triangulation that should not be, and I would like to sort them out by giving a constraint on the maximum edge length of the triangles. Is that possible?

Thanks for any help!

Best Answer

Is it possible? Of course. Compute the edge lengths. Decide which triangle have a longer edge length than what you want.

What will you do with that information? If the anser is you want the Delaunay triangulation to NOT create those triangles, and have it be smart enough to do that, think again. In context of a Delaunay triangulation, that is not an option you have. In fact, there is relatively little in context of a Delaunay triangulation you can do.

A Delaunay triangulation forms a triangulation of a region, such that the boundary of the region will be the convex hull of the data.

Given that, IF the region you want to see triangulated is not a convex region, then expect problems! You will get those long thin triangles. It will fill the gaps around the perimeter of your region where the data is NOT in fact convex.

Is your real question, can you create a triangulation that does NOT fill those gaps, thus not creating long, thin triangles around the perimeter? To some extent, yes, you do have an alternative - an alpha shape. An alpha shape starts with a Delaunay triangulation, then effectively rolls a ball (or sphere in 3-d) of radius alpha around the perimeter. Where the ball can penetrate a facet of the triangulation far enough to touch an interior node, it erodes those simplexes away. Depending on the options chosen, you can allow the code to erode interior holes or not.

xy = rand(100,2);
S = alphaShape(xy)
S = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: 0.0949990743504186
      HoleThreshold: 0
    RegionThreshold: 0
    
plot(S)

The code chose Alpha as 0.094999..., which is pretty small in general. I'd guess an alpha radius of vaguely 1/2 to 2 would be not be unreasonable, since the data is generated from uniform data in the unit square.

S.Alpha = .5
S = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: 0.5
      HoleThreshold: 0
    RegionThreshold: 0
    
plot(S)

So where the alpha ball can penetrate a perimeter facet of the triangulation far enough to TOUCH an internal node, the triangle will be eroded, thus deleted from what you get.

If you know in advance the value of alpha you might want to use, then pass it in directly when you create the shape. Thus use it as:

S = alphaShape(xy,0.5);

I generally find an alpha shape to be moderately insensitive to alpha radius as long as you don't go too small. Make the radius a little bit larger, and you will add a few moderately thin triangles. Go too large, and it starts to turn into the convex hull. Note that an alpha radius of inf will generally be equivalent to the Delaunay triangulation itself.

Finally, you do get two parameters to control hole and region thresholds. Without making this a complete tutorial on the idea, just read the docs for alphaShape. Play with it. Learn how to use the tool.

Related Solutions

MATLAB: Alpha shape vs convex hull volume value

A problem with your data is it has holes in it. The inside looks to be empty. But worse, it has a void on one corner, where the alpha ball can penetrate.

S0 = alphaShape(data(:,1:3))
S0 = 
  alphaShape with properties:
             Points: [76370×3 double]
              Alpha: 246.2
      HoleThreshold: 0
    RegionThreshold: 0
    
volume(S0)
ans =
   1.5087e+09

A volume of 1.5e9. But look at the alpha radius choen by alphaShape: 246.2. When I looked at the plot, the hole was roughly 5 times as large. Let me plot the result, and rotate it to show the interior:

plot(S0)
view(-2.5,11.5)

Do you see the hole? The result is completely HOLLOW.

Instead, if I use a larger value for alpha, I get a larger volume, closer to that of the convex hull.

S = alphaShape(data(:,1:3),1500);
volume(S)
ans =
   1.7524e+10

Admittedly, it does not resolve the details of that shape as well. But if you don't want the holes to allow the alpha ball to erode the entire interior of the volume, then you cannot have large holes in your data!!!

Perhaps you don't understand the idea of an alpha shape. You essentially start with the convex hull, as a delaunay triangulation. Then you roll a ball (of radius alpha) around the surface of the shape. Allow the ball to penetrate any facet of the surface, but do NOT allow it to pass through a node. So each simplex in the triangulation has 4 vertices, three of which can comprise an exterior facet. If the alpha ball can touch the 4TH vertex of a simplex on the surface, thus penetrating that simplex, then you delete the indicated simplex, removing it from the tessellation.

The idea is to remove large surface simplexes that form when we have a concave domain, because a basic delaunay triangulation ALWAYS describes a convex region.

Essentially, we allow the alpha ball to erode the tessellation. But what happens when you have a large hole in the cloud? A BAD IDEA. You will not be able to use an alpha shape with a small value of alpha, because you will see exactly what happens here.

xy = rand(1000,2);
k = sqrt(sum((xy - [1, .5]).^2,2)) < 0.4;
xy(k,:) = [];
plot(xy(:,1),xy(:,2),'.')

So a roughly square region, not too finely sampled. So there are some problems areound the edges, but mainly a hole of radius 0.4 intruding into one side. First, the delaunay triangulation will be equivalent to an alpha shape with an infinite radius alpha ball.

S0 = alphaShape(xy,inf)
area(S0)
ans =
       0.9783
       
plot(S0)

So a roughly square region. But we expect to have a region with area a little less than:

1 - pi*0.4^2/2
ans =
      0.74867

First, let alphaShape decide the radius itself.

S = alphaShape(xy)
S = 
  alphaShape with properties:
             Points: [733×2 double]
              Alpha: 0.035248
      HoleThreshold: 0
    RegionThreshold: 0
>> area(S)
ans =
      0.50374

It went too far. Note that there are subtly different ways to perform an alpha shape. My algorithm did not allow the ball to delete internal holes, alhthough I recall I offered that as an option.) The algorithm in alphaShape does allow that to happen. But I could have changed that behavior by setting the HoleThreshold.

Regardless, you can see the result used too small a value for alpha. Try it using a more reasonable alpha radius.

S = alphaShape(xy,.4)
area(S)
ans =
      0.75336
plot(S)

Here, you can see that I probably should have gone just a bit smaller for alpha, perhaps 0.3. Or, I could have used the HoleThreshold.

S = alphaShape(xy,0.1,'HoleThreshold',1)
S = 
  alphaShape with properties:
             Points: [733×2 double]
              Alpha: 0.1
      HoleThreshold: 1
    RegionThreshold: 0
    
plot(S)

To get a better approximation to the true area expected, I needed to use a reasonable alpha, but more importantly, I needed to more carefully sample the boundaries of the region.

So finally, let me do a good job of this.

xy = rand(100000,2);
k = sqrt(sum((xy - [1, .5]).^2,2)) < 0.4;
xy(k,:) = [];
S = alphaShape(xy,0.3,'HoleThreshold',1);
area(S)
ans =
      0.74726
      
plot(S)

As you can see, it quite reasonably approximates the area, which was analytically computed as 0.74867...

MATLAB: When using alpha shape and the alpha radius is not specified. How does MATLAB obtain it? I want to understand it.

I don't have access to the code, so I can only surmise (since I've written alpha shape codes before in MATLAB, but not the one in question) ...

From the help, we see:

SHP = alphaShape(X,Y,A) creates a 2D alpha shape of the points (X,Y), 
   with alpha radius A. X and Y are column vectors that specify the 
   coordinates. A is a positive scalar that defines the alpha radius. If 
   A is not specified, the default is the smallest alpha that produces an
   alpha shape enclosing all of the points.

The basic idea behind an alpha shape is that we start with a delaunay triangulation. We can think of that as an alpha shape with an alpha ball of infinite radius

xy = randn(100,2);
>> S = alphaShape(xy,inf)
S = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: Inf
      HoleThreshold: 0
    RegionThreshold: 0

The alpha ball erodes into the triangulation. If a ball of radius alpha can penetrate a triangle and touch an interior point, then the triangle is effectively deleted from the list. (It has been many years since I wrote alpha shape codes of my own, but this is the essential idea.)

We can see however, that if we reduce alpha from infinity, that gradually we will eliminate those triangles around the edges of the domain. That is the goal of an alpha shape. But how far do we go? The default, as we are told, is the radius that does not cause any points to be deleted from the shape.

S = alphaShape(xy)
S = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: 0.59772
      HoleThreshold: 0
    RegionThreshold: 0
    
plot(S)
hold on
plot(xy(:,1),xy(:,2),'ro')

So the default value was chosen to be 0.59772. Smaller than that? The alpha ball can erode some points effectively out of existence, by deleting all triangles tht contined those points s vertices.

S2 = alphaShape(xy,0.4)
S2 = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: 0.4
      HoleThreshold: 0
    RegionThreshold: 0
plot(S2)
hold on
plot(xy(:,1),xy(:,2),'ro')

However, again, I don't have the code. And it looks like the choice of a default for alpha is a moderately crude one. Since if I choose a slightly smaller value for alpha than the default, the shape still contains all the original points as vertices.

S3 = alphaShape(xy,0.58)
S3 = 
  alphaShape with properties:
             Points: [100×2 double]
              Alpha: 0.58
      HoleThreshold: 0
    RegionThreshold: 0

But if you plot it, visually, this seems to be identical to that from the default value. That means the phrase we read:

the default is the smallest alpha that produces an
   alpha shape enclosing all of the points

was only an approximation.

So clearly the scheme used to compute a default alpha is a moderately crude one, ad hoc, but not an optimization of any sort. I can imagine some schemes I might choose. What I need to recognize is the author of the code wanted to choose a scheme that was efficient. That suggests a sort of greedy algorithm. Once an alpha is identified that does reasonably well by some criteria, the search stops.

The question might be, what would I have done? I might start with the list of edges.

tri = DT.ConnectivityList;
edgelist = [tri(:,[1 2]);tri(:,[1 3]);tri(:,[2 3])];
edgelist = unique(sort(edgelist,2),'rows');
edgelens = sqrt(sum(xy(edgelist(:,1),:) - xy(edgelist(:,2),:)).^2,2);

Now, get the median edgelength.

median(edgelens)
ans =
      0.36336

So this gives me some measure of what a number that is surely going to be too small for alpha. I might try creating an alpha shape with that as alpha. If no points get deleted, then I would be done. Otherwise, I might just increase that value by some fraction, testing to look for a value where no points get deleted. But clearly the search that is done is an ad hoc, approximate one. And sadly, I cannot know better than that. (Since I am not an employee.)

Best Answer

Related Solutions

MATLAB: Alpha shape vs convex hull volume value

MATLAB: When using alpha shape and the alpha radius is not specified. How does MATLAB obtain it? I want to understand it.

Related Question