求空間直線的交點

1、因為一個小任務(wù)需要用到求空間直線交點的 MATLAB函數(shù)和求空間中兩個 平面的相交線的函數(shù)，但是在網(wǎng)上找了一下沒有，只好自己寫了幾個函數(shù)，自 己覺得還能用，在這里給大家分享一下。1.MATLAB求空間中的兩個平面的相交線functionflag,n,p = Plane2Line(n1,p1,n2,p2)% calulate the line shared by two intersecting plane% input:% n1 normal vector of plane one% p1 any point on plane one% n2 normal vector of plane tw

2、o% p2 any point on plane two% output:% flag whether the two planes are intersecting ( 1 or 0)% n the direction vector of the expected line% p any point in the expected line% author:Lai Zhenzhou from Harbin Institute of Technology% email:% date:2014.1.16%if(isvector(n1) & isvector(p1) & isvector(n2)

3、&isvector(p2)error(Plane2Line:1 / 12the parameter is not vector);endif(length(n1)=3)|(length(p1)=3)|(length(n2)=3)|(length(p2)=3) error(Plane2Line:the parameter is not 3d vector);endA = n1(1) n1(2) n1(3) ;n2(1) n2(2) n2(3) ;if(rank(A)2)flag = 0;elseflag = 1;endif(flag=1)% calculate the normal vector

4、n = cross(n1,n2);2 / 12c1 = n1(1) n1(2) n1(3) -dot(n1,p1);c2 = n2(1) n2(2) n2(3) -dot(n2,p2);% calculate the simplest Row echelon matrix temp1 = rref(A);temp2 = 1 2 3;index(1) = find(temp1(1,:),1,first);% the first nonzero elementindex(2) = find(temp1(2,:),1,first);% the first nonzero elementD = A(:

5、 ,index(1) ) A(:,index(2) );Y = dot(n1,p1);dot(n2,p2);3 / 12X = inv(D)*Y; fori=1:3 if(i=index(1) & i=index(2) ) index(3) = i;endend p(index(1) ) = X(1) ; p(index(2) ) = X(2) ; p(index(3) ) = 0; else n=;p=; end4 / 12for test% flag n p = Plane2Line(1 2 3,1 0 1,2 3 4,0 -1 0)% flag n p = Plane2Line(0 0

6、1,0.5 0.5 0.5,1 0 0,0 0 1)2.MATLAB求空間中兩條直線的交點functionflag,p = Line2Point(n1,p1,n2,p2)% determine the relation between two straight lines and% calulate the intersection point if they are intersecting% input:% n1 direction vector of line one% p1 any point in line one% n2 direction vector of line two%

7、p2 any point in line two% output:% flag the relation of the two line% flag = 0 the two line are on different plane% flag=1 thetwo line are on the same plane and they are parallel% flag = 2 the two line are on the same plane and they areintersecting% pthe point shared by the two intersecting line% au

8、thor:Lai Zhenzhou from Harbin Institute of Technology% email:% date:2014.1.17% reference:http:5 / 12%if(isvector(n1) & isvector(p1) & isvector(n2) & isvector(p2)error(Line2Point:the parameter is not vector);endif(length(n1)=3)|(length(p1)=3)|(length(n2)=3)|(length(p2)=3)error(Line2Point:the paramete

9、r is not 3d vector);endA = p2(1) -p1(1) p2(2) -p1(2) p2(3) -p1(3) ;n1(1) n1(2) n1(3) ;n2(1) n26 / 12n2(2)(3) ; if(det(A)=0) flag = 0; else if(rank(A(2:3,: )2) flag = 1; else flag = 2; end end if(flag = 2) B = rref(A(2:3,: );index(1) = find(B(1,: ),1,first); index(2) = find(B(2,: ),1,first);7 / 12for

10、i=1:3if(i=index(1) & i=index(2) )index(3) = i;endendY(1,1) = -p1(index(1) ) + p2(index(1) );Y(2,1) = -p1(index(2) ) + p2(index(2) );D = n1(index(1) -n2(index(1);n1(index(2) ) -n2(index(2);t = inv(D)*Y;p(1) = p1(1) + n1(1)*t8 / 12(1);p(2) = p1(2) + n1(2)*t(1);p(3) = p1(3) + n1(3) *t(1);elsep = ;enden

11、d%for test% flag p = Line2Point(1 0 0,0 0 0,0 1 0,1 1 0)% flag p = Line2Point(1 1 0,0 0 0,0 1 0,1 1 0)% flag p = Line2Point(1 1 -1,0 0 1,0 1 0,1 1 0)% flag p = Line2Point(1 1 - 1,0 0 1,1 1 -1,1 1 0)% flag p = Line2Point(1 1 1,0 0 0,2 0 0,2 1 1)% flagp=Line2Point(0-42050 -21025,-100,00145,00 0)3.MATL

12、AB求空間中直線和某一線段的交點functionflag,p = LineSegment2Point(n1,p1,p21,p22)% determine the relation between a straight line and% a line segment and calulate the intersecting point9 / 12% if they are intersecting% input:% n1 direction vector of the straight line% p1 any point in the straight line% p21 any poin

13、t in the line segment% p22 any another different point in the line segment% output:% flag the relation of the two line% flag=0 thestraight line and line segment are not intersecting% flag = 1 the straight line and line segment are intersecting% p the point shared by the two intersecting line% author

14、:Lai Zhenzhou from Harbin Institute of Technology% email:% date:2014.1.17%if(isvector(n1) & isvector(p1) & isvector(p21) & isvector(p22)error(LineSegment2Point:the parameter is not vector);endif(length(n1)=3)|(length(p1)=3)|(length(p21)=3)|(length(p22)=3)error(LineSegment2Point:the parameter is not

15、3d vector);endif(p21(1)=p2210 / 12(1)&(p21(2)=p22(2)&(p21(3) =p22(3) )error(LineSegment2Point:the two appointed points of line segment are the same);end%calculate the direction vector of line segmentn2 = p22(1)-p21(1),p22(2)-p21(2),p22(3) -p21(3) ;flag0,p0 = Line2Point(n1,p1,n2,p21);if(flag0 = 2)if(norm(p0-p21)=norm(p22-p21) & norm(p0-p22)=norm(p22-p21)flag = 1; p = p0;else11 / 12flag = 0;p=;endelseflag