Subtract delaunay surfaces from each other

Question

I have a surface by the code below and another surface which is created by the exact same code. I want to see the height differences in another figure. How am I able to do that? Already operated with the Minus-operator but this won't work.

Furthermore the matrices have NOT the same size!

Appreciate your help!

x1 = Cx1;
y1 = Cy1;
z1 = Cz1;
tri1 = delaunay(x1,y1);


fig1 = figure%('units','normalized','outerposition',[0 0 1 1]);
trisurf(tri1,x2,y2,z2)
xlabel('x [mm] ','FontSize',30)
ylabel('y [mm] ','FontSize',30)
zlabel('z [mm] ','FontSize',30)

You first need to define how you would like the difference to be output. For instance, do you want to map one surface to the other one and then do the difference? Or do you want a polygonal surface that maintains full resolution of both? Or, something else? — TTT, Dec 07 '16 at 17:23
Thank you for your answer. Oh okay. What would you suggest? I would love to learn both ways if possible. So we can kill two birds with one stone in this thread :) — MatlabNewb, Dec 07 '16 at 17:27
Creating a new polygonal surface is no simple task. But no matter how you do it, you need to define functionally how values vary across your triangular elements. I believe the `trisurf` function assumes linear variation. Then, you can calculate values at any given point on the triangle, and thereby interpolate between the two triangulations. — TTT, Dec 07 '16 at 17:46
Can you show an example please? That's why I was asking :) I would like to see the difference in height in a seperate plot. That would be awesome! — MatlabNewb, Dec 07 '16 at 18:00

score 0 · Answer 1 · edited May 23 '17 at 12:01

The simplest way to solve this problem is to interpolate from one mesh onto the other one. Such an approach works well when one is more highly resolved than the other, or when you're not as concerned with results at individual nodes, but rather the overall pattern across elements. If that's not the case, then you have a very complicated problem because you need to create a polygonal surface that fully captures all nodes and edges of both triangulations. Consider the following pair of triangular patterns:

A surface that captured all the variations would need to have all the vertices and edges that make up both of them, which is not a purely triangular surface. So, let us instead assume the easier case. To map results from one triangulation to the other, you simply need to formulate functions that define how the values vary along the triangles, which are more broadly called basis functions. It is often assumed that values betweeen the nodes (i.e. vertices) of the triangles vary linearly along the surfaces of the triangles. You can do it differently if you want, it just requires defining new basis functions. If we go for linear functions, then the equations in 2D are pretty simple. Let's say you make an array trimap that has which triangle each of the vertices of the other triangulation is inside of. This can be accomplished using the info here. Then, we set the coordinates of the vertices of the current triangle to (x1,y1), (x2,y2), and (x3,y3), and then do the math:

for cnt1=1,npoints
  x1=x(tri1(trimap(cnt1),1));
  x2=x(tri1(trimap(cnt1),2));
  x3=x(tri1(trimap(cnt1),3));
  y1=y(tri1(trimap(cnt1),1));
  y2=y(tri1(trimap(cnt1),2));
  y3=y(tri1(trimap(cnt1),3));
  delta=x2*y3+x1*y2+x3*y1-x2*y1-x1*y3-x3*y2;
  delta1=(x2*y3-x3*y2+xstat(cnt1)*(y2-y3)+ystat(cnt1)*(x3-x2));
  delta2=(x3*y1-x1*y3+xstat(cnt1)*(y3-y1)+ystat(cnt1)*(x1-x3));
  delta3=(x1*y2-x2*y1+xstat(cnt1)*(y1-y2)+ystat(cnt1)*(x2-x1));
  weights(cnt1,1)=delta1/delta;
  weights(cnt1,2)=delta2/delta;
  weights(cnt1,3)=delta3/delta;
  z1=z(tri1(trimap(cnt1),1));
  z2=z(tri1(trimap(cnt1),2));
  z3=z(tri1(trimap(cnt1),3));
  valinterp(cnt1)=sum(weights(cnt1,:).*[z1,z2,z3]);
end

valinterp is the interpolated value for each point. Here and here are some nice slides explaining the mathematics behind all this. Note that I've not tested any of this code. Note also that you will need to do something to assign to values outside of the triangulation. Perhaps a null value, or an inverse-distance weighted value.

Subtract delaunay surfaces from each other

1 Answers1