How can I find a fit to my data points and using triangulation to approximate the planar area?

I got these data points, I want to approx the planar area of the shape.

I am thinking of finding a fit then approxmate from the function but I am not sure.

Answer 1

[John D'Errico]

I don't have your data, while they say a picture is worth a thousand words, it is still just words, not numbers. Anyway, I think you are asking to find the cancave region where the points live, but not including those triangular chaped voids.

For this purpose, an alpha shape is often the best approach. The idea is to visualize a circle or ball, of radius r. Now imagine that ball rolling around the perimeter of your points, in a way where the ball cannot cross any data point. Allow the alpha ball to erode a triangulation to the extent that it can. So any triangles around the perimeter with long edges will be deleted. Again, since I lack any data, I'll make some up.

xy = rand(5000,2);
k = (sqrt(sum((xy - [1 .5]).^2,2)) < 0.3) | ...
    (sqrt(sum((xy - [0 .5]).^2,2)) < 0.3) | ...
    (sqrt(sum((xy - [0.5 1]).^2,2)) < 0.3) | ...
    (sqrt(sum((xy - [0.5 0]).^2,2)) < 0.3);
xy(k,:) = [];
plot(xy(:,1),xy(:,2),'.')
grid on

Admittedly, not a very good approximation to what you had. But it is early in the morning, and I don't care about the difference.

S = alphaShape(xy(:,1),xy(:,2),.20)
S = 
  alphaShape with properties:

             Points: [2163×2 double]
              Alpha: 0.2
      HoleThreshold: 0
    RegionThreshold: 0
plot(S)

I choise 0.2 for the alpha radius arbitrarily, assuming I did not know the diameter of the circle I had actually used. Anyway, in real life, you will never have exactly circular cutouts as I did here. Regardless, it is not a poor approximation, but only that. See that if I chose the alpha radius to be too large, it will start to leave large triangles aound the perimeter. Too small an alpha radius, and we will get too much erosion of the real domain we are interested in. Luckily with alpha shapes, there tends to be a fuzzy region where you can be a little wrong and it won't matter much. A problem can arise where you have sharp internal corners, but thats life.

What is the area contained though?

area(S)
ans =
      0.42664

Lets see, since I know the actual shape of the region, I can easily compute the area as:

1 - 4*(pi*0.3^2/2)
ans =
      0.43451

Not bad at all.

By way of comparison, consider what a Delaunay triangulation does on the same data.

D = delaunayTriangulation(xy);
triplot(D)

As you can see, it is those triangles with the long edges around the pererimeter that we need to erode away, and the alpha shape did that spectacularly well. Remember that a Delaunay triangulation will always create a convex domain.

I would point out that with fewer points, the underlying delaunay triangulation that is eroded to create the alpha shape will be a coarse one, and therefore you may have some issues if you make alpha too small. But since you do have moderately sharp internal corners there, you want a small alpha ball. Hey, the real world is like that, and I've only once ever met a client in my career as a consultant who offered more data than I really wanted to do the best job possible. Remember that your eyes/brain form a very good pattern recognition utility. They do easily what computers often need real work to do well.

[fengsen huang]

Thank You so much and sorry about posting it twice.

About not enough data I attached the x and y data so you can understand what the shape look like when plotted. However, is it still an alpha shape?

[fengsen huang]

I used alphashape but the graph has holes in it?

shp = alphaShape(JJ,KK);

pc = criticalAlpha(shp,'all points');

shp.Alpha = pc;

plot(shp)

[John D'Errico]

Sigh. Did you read what I said? If you choose too small of an alpha, then you get crap. You need to use a larger alpha, or get more finely resolved data.

So first, you need to use asufficiently large alpha radius. Next, you can use the alphaShape properties to kill off the internal voids. READ THE HELP!

So I tried an alpha radius of 3. Then since pi*3^2 is roughly 30, a HoleThreshold of 30 should do a decent job.

S = alphaShape(JJ,KK,3,'HoleThreshold',30),plot(S)
S = 
  alphaShape with properties:

             Points: [1144×2 double]
              Alpha: 3
      HoleThreshold: 30
    RegionThreshold: 0
    

area(S)
ans =
          3568.08620000001

If you feel you needed a better approximation to thearea, then you need to sample more points to resolve that region more accurately.

[fengsen huang]

i think i just need more points plotted to get an accurate answers.

Now im just considering a single quadrant and I added a lot more points to that quadrant because basically the grapgh is a fold fold symmetry graph. I used the absolute code to convert it all to the positive quadrant. The question is now i got the area of it but how can i find out the volume of it. I got the depth of each points so basically i got the curve plotting on the same graph, one on top and one below, same x axis but y axis changes because of the thickness. How can I and use a surf plot to create a volume graph and thus find out the volume?

[Bruno Luong]

You can compute the volume by summing the volume of each triangle+depth #i

V_i = 1/3 surface(T_i) * (H_i1 + H_i2 + Hi_3)

Where (H_i1,H_i2,Hi_3) are the depths at the 3 vertexes of the triangle V_i.

The sum them for i that returns by alphashape or my method.

V = sum_i (V_i)

The area of the triangles T where the vertexes are three 2x1 vectors A,B,C is

1/2 * abs(det([B-A,C-B]))

[fengsen huang]

Hi,

I got a matrix with 3532 x 6 cell

where collumn 2 is the x axis

collumn 3 is the y axis

collumn 4 shows if there are anything underneath

collumn 5 shows the top surface depth below ground

collumn 6 shows the depth of the object

so collumn 5 and 6 shows the thickness of the object

so i dont know if your method works

[Bruno Luong]

My method describe the mathematics formula. Just plug your data appropriately.

You hav x,y (columns 2 & 3), which triangles is on the area (T in my code) and the connectivity vertexes; the height of your object (columns 5 & 6). So everything are there.

[fengsen huang]

sorry I don't quite understand.

I baiscally have these datas attached then below is my code and a plot at the end. Once you run the plot you have some idea of what the 3D shape look like. I baiscally want to make a trisurf plot or something to show its 3d solid shape thus calcualting the volume.

X1 = cell2mat(C(:,2));

Y1 = cell2mat(C(:,3));

d1 = cell2mat(C(:,5));

d2 = cell2mat(C(:,6));

t = d1 - d2;

T = nonzeros(t);

m = (t>0);

XX = X1(m); JJ = double(XX);

YY = Y1(m); KK = double(YY);

D1=nonzeros(d1);

D2=nonzeros(d2);

jj=abs(JJ); kk=abs(KK);

plot3(jj,kk,D1,'kx',jj,kk,D2,'bx');

[Bruno Luong]

I simply suggest you the method to compute a volume. You do not have to plot anything to compute the volume.

What s the volume of the box of sides a,,b,c. It's V=a*b*c. Do you need to plot the box to know the volume? No.

If you want use something else then feel free to do it, and you just ignore what I suggest to you.

[fengsen huang]

With the variables I have now I dont know how to find the volume, I am not trying to ignore your code, I am tryingto keep up to it but dont understand. You said something about triangles??

[Bruno Luong]

Answer you bellow, on comment of my code

Answer 2

[KSSV]

Read about boundary and polyarea. Let (x,y) be your data.

idx = boundry(x,y) ;
bx = x(idx)  ; by = y(idx) ;
A = polyarea(bx,by) ;

If the above is not working....have a look on delaunaytriangulation. It will triangulate the area and also gives you area.

[fengsen huang]

Hi,

I try to use delaunay triangulation but it also connected the boundary point

As youc an see the X shape, thats what I wanted only not the square too.

thank you

[Bruno Luong]

Arguably one can see where is the area outside with a more or less clear X shape.

But if you ask 10 peoples to draw the border, you probably get 10 different answers in the inside part.

The computer might be get confused like those 10 persons.

[KSSV]

alphashape is the solution....

[Bruno Luong]

Honestly, such kind of problem is not very well defined and I'm sure I want to pretend anything is THE solution.

Beside that I have never been a big fan of alphashape myself.

[John D'Errico]

I've used them heavily over many years in both 2 and 3 dimensions. In fact, I wrote working versions of them around 20 years ago in MATLAB. They work quite well.

[Bruno Luong]

Hi John, I don't deny it works well. But

I have problem with the curvature criteria.

To me a outer surface of an object can have big curvaturre, like a corner of the box. Fixing a threhold on that seems to me then just "try it and see until you are happy", very operator dependent. The second reason is also on curvature but in 3D, anypoint have 2 separate principal curvatures, and alphashape just smooths them out by a sphere.

But it's just me.

Answer 3

[Bjorn Gustavsson]

One thing you could try to squeeze around some difficult parts of your problem is to run through a coordinate transformation - for example to cylindrical coordinates (r, theta), after plotting your points in that coordinate system perhaps a mathematically simpler paramterization of your boundary can be found that way.

HTH

[John D'Errico]

This might work if done carefully. It would probably be a periodic spline approximation to the envelope in r, were I to do it. Of course then, the area computation becomes an easy thing. A virtue of the approach is you can also do smoothing on that envelope if you wanted. The biggest problem I would foresee is the envelope function, thus r(theta) may not be a single valued function of polar angle theta. That appears in fact to be a potential problem for the 4 leafed clover region we have seen.

[Bjorn Gustavsson]

If one are a bit lucky it might work even with a QD first stab. I got something reasonably similar with

x_perim = cos(theta).*sin(n*theta);

y_perim = sin(theta).*sin(n*theta);

...and technically that is not a preimeter of a simple closed region, but it would be pretty

straightforward to fix that and adjust things for the regions where the distance to the centre

is small. If this is how the original problem was designed it is not all that very interesting - the general problem is both more difficult and interesting .

Answer 4

[Bruno Luong]

Quick and dirty fix, just discard outside delaunay triangles with long sides

load('x-data.mat');
load('y-data.mat');
xy = [JJ KK];

DT = delaunayTriangulation(xy);
K = convexHull(DT);
% adjustable threshold
dmax = 6;

T = DT.ConnectivityList;
T = sortrows(sort(T,2));
while true
    xyb = xy(K,:);
    d = diff(xyb,1,1);
    d = sqrt(sum(d.^2,2));
    Tbad = find(d > dmax);
    ibad = Tbad+[0 1];
    Kbad = reshape(K(ibad),size(ibad));
    K12 = sort(Kbad,2);
    [b12,i12] = ismember(T(:,[1 2]),K12,'rows');
    [b23,i23] = ismember(T(:,[2 3]),K12,'rows');
    [b13,i13] = ismember(T(:,[1 3]),K12,'rows');
    i12 = i12(b12);
    i23 = i23(b23);
    i13 = i13(b13);
    k12 = T(b12,3);
    k23 = T(b23,1);
    k13 = T(b13,2);
    b = b12 | b23 | b13;
    K3 = accumarray([i12;i23;i13],[k12;k23;k13]);
    if isempty(K3)
        break
    end
    T(b,:) = [];
    Tbad = Tbad(:)';
    for i=length(K3):-1:1
        K3i = K3(i);
        if K3i == 0
            continue
        end
        j = Tbad(i);
        K = [K(1:j); K3i; K(j+1:end)];
    end
end

close all

trisurf(T,DT.Points(:,1),DT.Points(:,2),zeros(size(DT.Points(:,2))));
view(2)
axis equal

[Bruno Luong]

Compute the volume:

load('all data.mat')

X1 = cell2mat(C(:,2));
Y1 = cell2mat(C(:,3));
d1 = cell2mat(C(:,5));
d2 = cell2mat(C(:,6));
t = d1 - d2;
m = (t>0);
XX = X1(m); JJ = double(XX);
YY = Y1(m); KK = double(YY);
D1=nonzeros(d1); 
D2=nonzeros(d2);
jj=abs(JJ); kk=abs(KK);

xy = [jj,kk];
DT = delaunayTriangulation(xy);
K = convexHull(DT);
% adjustable threshold
dmax = 6;
T = DT.ConnectivityList;
T = sortrows(sort(T,2));
while true
    xyb = xy(K,:);
    d = diff(xyb,1,1);
    d = sqrt(sum(d.^2,2));
    Tbad = find(d > dmax);
    ibad = Tbad+[0 1];
    Kbad = reshape(K(ibad),size(ibad));
    K12 = sort(Kbad,2);
    [b12,i12] = ismember(T(:,[1 2]),K12,'rows');
    [b23,i23] = ismember(T(:,[2 3]),K12,'rows');
    [b13,i13] = ismember(T(:,[1 3]),K12,'rows');
    i12 = i12(b12);
    i23 = i23(b23);
    i13 = i13(b13);
    k12 = T(b12,3);
    k23 = T(b23,1);
    k13 = T(b13,2);
    b = b12 | b23 | b13;
    K3 = accumarray([i12;i23;i13],[k12;k23;k13]);
    if isempty(K3)
        break
    end
    T(b,:) = [];
    Tbad = Tbad(:)';
    for i=length(K3):-1:1
        K3i = K3(i);
        if K3i == 0
            continue
        end
        j = Tbad(i);
        K = [K(1:j); K3i; K(j+1:end)];
    end
end

% close all
% hold on
% trisurf(T,DT.Points(:,1),DT.Points(:,2),D1);
% trisurf(T,DT.Points(:,1),DT.Points(:,2),D2);
% view(3)
% axis equal

H = abs(D2-D1);
ABC = xy(T,:);
ABC = reshape(ABC,[size(T),2]);
ABC = permute(ABC,[3 2 1]);
ABC = ABC(:,2:3,:)-ABC(:,1,:);
% Triangle areas
S = 0.5*abs(ABC(1,1,:).*ABC(2,2,:)-ABC(1,2,:).*ABC(2,1,:));
HT = 1/3*sum(H(T),2);
S = reshape(S,size(HT));
V = sum(S.*HT);

fprintf('Volume = %g\n', V);