Intersection of a surface generated by scattered points and a line

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); 

figure 
trisurf(tri,Xb(:,1),Xb(:,2),Xb(:,3), 'FaceColor', 'cyan', 'faceAlpha', 0.8);

Any help will be much apperacitated.

Answer 1

[Roger Stafford]

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 & ...
 dot(cross(P1-P3,P0-P3),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.

Answer 2

[Bladeyus]

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 & ...
   dot(cross(P1-P3,P0-P3),N)>=0
end
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?

[Roger Stafford]

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.

 clear
 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 & ...
     dot(cross(P1-P3,P0-P3),N)>=0
   n = n+1;
   X(n) = P0(1); Y(n) = P0(2); Z(n) = P0(3); % If so, store P0
  end
 end
 plot3([P1(1),P2(1),P3(1),P1(1)],[P1(2),P2(2),P3(2),P1(2)],...
       [P1(3),P2(3),P3(3),P1(3)],'m-',...
       P1(1),P1(2),P1(3),'ro',P2(1),P2(2),P2(3),'go',...
       P3(1),P3(2),P3(3),'bo',...
       X,Y,Z,'y.')

Answer 3

[Bladeyus]

Thank you very much for the response. I will try it as soon as possible and share the result here.

Answer 4

[Bladeyus]

Thanks Roger,

I finally figured it out. Sorry for disturbing so much as i am new to MatLab.