Create intersection between surface and oriented planes

Hi, i am new in Matlab. I would like to know the easiest way to obtain intersection profiles (sort of cross sections) between a complex surface obtained by a point cloud and a series of vertical planes oriented at defined angles compared to the grid and spacing. The purpose will be to perform a best fit analysis of each intersection profile, compare and analyse the results. Thanks and regards

Fabio

Answer 1

[Mike Garrity]

How is your surface defined? A patch object? Are the faces triangles, or perhaps quads?

This blog post covers a related case. You'll need to change things a bit to fit what you're doing.

The basic idea is that you use the plane equation to compute a signed distance from the plane at each vertex in your mesh. Depending on the representation of your planes, that's probably something like:

for i=1:nvert
  d(i) = dot(vert(i,:), plane(1:3)) - plane(4);
end

Now you need to look at each face. I'll consider the triangular face case, since it's the simplest. If the 3 vertices of the face all have positive or all have negative distances, then you can ignore the face. If the one vertex has the opposite sign of the other two, then you need to look at the two edges that connect that vertex to the other two. Let's say that we know that vertex i2 has opposite sign from vertex i1 (i1 & i2 are integer indices into the vertex array). Then that edge crosses the plane at the vertex defined by:

t = (0 - d(i2)) / (d(i1) - d(i2));  % between 0 & 1
pt = vert(i2,:) + t*(vert(i1,:) - vert(i2,:));

Since you're doing this for two edges on that triangular face, you'll get two vertices. Connect those with a line. Once you've done this for all of your faces, you'll have a set of lines which trace out the intersection.

If your faces are quads, then it's basically the same, but you'll have some special cases and degenerate to worry about.

Depending on your usage, you might want to connect the lines to form a continuous curve. That can be a bit tricky. One way is to do a spatial search to connect them all. That gets really expensive.

A more efficient way is to an old MATLAB graphics programming trick. Create a sparse array. Every time you create a line segment, save it in an array and put the index into that array into the sparse matrix, indexed by the i1 and i2 you used to create each of the end points. Now, when you create another line from an edge between those two vertices, you can use the sparse array to find the first one and then you can connect them. Did that make any sense?

I hope that helps. Just ask here on MATLAB Answers if you get stuck.

[Fabio Lapponi]

Hi, thanks for your answer. Lots to figure out. My surface was created with meshgrid and griddata from a x,y,z from a topographic surface. Soit is not expressed by an implicit function. My intention is to figure out a way to create a series of intersection with vertical planes spaced regularly and at a certain angle (azimuth) n order to obtain a series of 2d plits to perform curve fitting and other calculations. Maybe the easiest way wold be to rotate the grid to my angle of interest? Thanks and regards

[Mike Garrity]

Oh, that sounds good. You have a special case where there are shortcuts available.

Let's see if I can explain this quickly. I should probably turn into another post for the graphics blog because it's a bit involved.

You've got a 2D array of X,Y, & Z coordinates. You can substitute these into your plane equation and get a 2D array of signed distances from the plane. Now, there's a function called contour which knows how to generate isocurves from 2D arrays. We can ask it to create the curves for where your 2D signed distance array equals 0 like this:

contour(x,y,d,'LevelList',0)

The problem is that it's just going to draw this, and it's not very useful because it would be the projection of the intersection onto the XY plane. But the object that contour returns has a property called ContourMatrix. We can use this to generate the curves you want. This property has kind of a strange format. It has two rows. The value in c(1,1) is the level (in this case 0), and the value in c(2,1) is an integer which is the number of points in the curve. Let's call that n. In the next n columns, the first row is X and the second row is Y. After that, you'll have another column which is a value and n again, if you have more than one intersection curve.

The resulting code would look something like this:

% Generate X,Y, & Z values and draw them as a surface.
grey = [.875 .875 .875];
[x,y,z] = peaks;
surf(x,y,z,'FaceColor',grey)

% Create a random plane.
rng default
plane = [randn(1,3), 0];

% Use X,Y,Z and the plane to compute the signed distance.
d = x*plane(1) + y*plane(2) + z*plane(3) - plane(4);

% Use contour to get the ContourMatrix.
hold on
[~,hcontour] = contour(x,y,d,'LevelList',0);
c = hcontour.ContourMatrix;
delete(hcontour);

% Loop through the ContourMatrix
i = 1;
while i<size(c,2)
  % Get the X & Y for the next curve
  n = c(2,i);
  x2 = c(1,i+(1:n));
  y2 = c(2,i+(1:n));
  % Use interp2 to compute the matching Z values
  z2 = interp2(x,y,z, x2,y2);
  % Draw that line
  line(x2,y2,z2,'Color','red','LineWidth',2);
  % Advance to next curve
  i = i+n+1;
end

I hope that helps.

[Mike Garrity]

I should point out that I could also have computed z2 there by substituting x2 and y2 into the plane equation and solving for z. That would probably be much more robust than using interp2, except for the case where the plane is almost normal to the XY plane, which is the case I think you're talking about.

[Vignesh Suresh]

Hi Mike, I am kind of stuck with the same problem. I need to extract the intersection profile of the plane and surface. How to extract the intersection zone as a separate 2D graph.