Intersection of a surface generated by scattered points and a line

44 views (last 30 days)
Hello all,
Actually i am new to MatLab and i am trying to figure out if there is a way of finding the intersection of a 3d surface (coordinates of the surface are read from excel and does not follow a specific formula, therefore i cant calculate the equation of the surface) and a line? I am using below code to generate my surface.
Y1= [xlsread(filename,'AI6:AI63')];
Z1= [xlsread(filename,'AG6:AG63')];
X1= [xlsread(filename,'AH6:AH63')];
dt = DelaunayTri(X1,Y1,Z1);
[tri Xb]= freeBoundary(dt);
trisurf(tri,Xb(:,1),Xb(:,2),Xb(:,3), 'FaceColor', 'cyan', 'faceAlpha', 0.8);
Any help will be much apperacitated.

Accepted Answer

Roger Stafford
Roger Stafford on 7 May 2013
Edited: Roger Stafford on 7 May 2013
The intersection of a triangle and a line in three dimensional space can be found as follows. Using three cartesian coordinates, let vectors P1, P2, and P3 be the triangle's three vertices, and let vectors Q1 and Q2 be any two points along the given line. Then do this in matlab:
N = cross(P2-P1,P3-P1); % Normal to the plane of the triangle
P0 = Q1 + dot(P1-Q1,N)/dot(Q2-Q1,N)*(Q2-Q1); % The point of intersection
P0 will be the intersection of the line and the infinite plane of the triangle. To see if P0 lies within the triangle's perimeter, the following must hold true:
dot(cross(P2-P1,P0-P1),N)>=0 & ...
dot(cross(P3-P2,P0-P2),N)>=0 & ...
To solve your problem as posed you can test the given line against each of the triangles in your triangulation in this way to see which, if any, triangles it intersects. This will require extracting each of the triangles' set of three vertices from the triangulation object output by DelaunayTri.

More Answers (3)

Bladeyus on 8 May 2013
Thanks to you Roger i could be able to make some progress. The final stage of my code is below.
X1= [xlsread(filename,'AO6:AO13')];
Y1= [xlsread(filename,'AP6:AP13')];
Z1= [xlsread(filename,'AQ6:AQ13')];
dt = DelaunayTri(X1,Y1,Z1);
[tri Xb]= freeBoundary(dt);
numtri = size(tri,1);
Q1=[0 0 0]
Q2=[0 0 50]
for j=1:numtri
P1=[X1(tri(j,1),1) Y1(tri(j,1),1) Z1(tri(j,1),1)]
P2=[X1(tri(j,2),1) Y1(tri(j,2),1) Z1(tri(j,2),1)]
P3=[X1(tri(j,3),1) Y1(tri(j,3),1) Z1(tri(j,3),1)]
N = cross(P2-P1,P3-P1); % Normal to the plane of the triangle
P0 = Q1 + dot(P1-Q1,N)/dot(Q2-Q1,N)*(Q2-Q1); % The point of intersection
dot(cross(P2-P1,P0-P1),N)>=0 & ...
dot(cross(P3-P2,P0-P2),N)>=0 & ...
and my X1, Y1,Z1 values are
20 20 100
-20 20 100
-8 -20 100
8 -20 100
20 20 -100
-20 20 -100
-8 -20 -100
8 -20 -100
The last problem is even if the line does not intersect with the surfaces (does not go outside the surface) it returns answer as 1. And finds an intersecion point P0 = [0 0 -100]. I have tried other values of Q1 and Q1 but unfortunately still the same. Do you have any idea why?
  1 Comment
Roger Stafford
Roger Stafford on 9 May 2013
I have no way of determining which are the triangles from the values of X1, Y1, Z1 you show, so I can't check your work. It looks entirely possible for a triangle among three of the last four of your vertices to actually intersect your line. You should check your triangulation to see if it matches your idea of the surface. Also be aware that the line Q1 and Q2 define is considered by this algorithm to extend to infinity in both directions.
I don't understand the code you show - the test for whether the line intersects the j-th triangle doesn't appear to be stored anywhere. Is that the way you actually did your computation? You need to detect whether that logic proposition is true at each j-th pass with an 'if' statement. Whenever it is true you should store the corresponding P0 somewhere. At the end of your run if any P0's have been stored in this way, these are your intersection points.
Before I answered your original request I tested the expressions in it. To increase your confidence in the method I have included the code for that test which you can run for yourself. It randomly selects the three vertices P1, P2, and P3 for a single fixed triangle and then repeatedly generates random Q1 and Q2 points for various lines. If the line happens to pass through the triangle interior, the intersection point P0 is stored in the X,Y,Z arrays. When 8192 intersection points have been stored, it uses 'plot3' to display the triangle with a magenta outline and vertices in red, green, and blue and with the 8192 intersection points in yellow. As you can see, the "yellow" points all lie within the triangle. Of course there were a great many other intersections with the triangle's plane which lay outside the triangle but these are not shown since they didn't pass the test.
P1 = randn(3,1); P2 = randn(3,1); P3 = randn(3,1); % The triangle
N = cross(P2-P1,P3-P1); % The normal to plane of the triangle
n = 0; ne = 8192;
X = zeros(ne,1); Y = zeros(ne,1); Z = zeros(ne,1);
while n<ne
Q1 = randn(3,1); Q2 = randn(3,1); % Generate a new line each time
P0 = Q1 + dot(P1-Q1,N)/dot(Q2-Q1,N)*(Q2-Q1); % Get intersection
if dot(cross(P2-P1,P0-P1),N)>=0 & ... % Is P0 is inside triangle?
dot(cross(P3-P2,P0-P2),N)>=0 & ...
n = n+1;
X(n) = P0(1); Y(n) = P0(2); Z(n) = P0(3); % If so, store P0

Sign in to comment.

Bladeyus on 7 May 2013
Thank you very much for the response. I will try it as soon as possible and share the result here.

Bladeyus on 9 May 2013
Thanks Roger,
I finally figured it out. Sorry for disturbing so much as i am new to MatLab.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!