source: issm/trunk/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp@ 6088

Last change on this file since 6088 was 6088, checked in by Eric.Larour, 14 years ago

Updated all routines to take into account ISSM_DIR AND ISSM_TIER environment variables
File size: 2.6 KB
Line 
1/*! \file SegmentIntersect.cpp
2*/
3
4#include "./MeshProfileIntersectionx.h"
5
6int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){
7
8 /*See ISSM_TIER/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
9
10 /*output: */
11 double alpha=-1;
12 double beta=-1;
13
14 double xA,xB,xC,xD,yA,yB,yC,yD;
15 double O2A[2],O2B[2],O1C[2],O1D[2];
16 double n1[2],n2[2];
17 double test1, test2, test3, test4;
18 double det;
19 double O2O1[2];
20 double pO1A,pO1B,pO1C,pO1D;
21
22 xA=x1[0]; yA=y1[0];
23 xB=x1[1]; yB=y1[1];
24 xC=x2[0]; yC=y2[0];
25 xD=x2[1]; yD=y2[1];
26
27 O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
28 O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
29 O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
30 O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
31
32
33 n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA
34 n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB
35
36 test1=n2[0]*O2A[0]+n2[1]*O2A[1];
37 test2=n2[0]*O2B[0]+n2[1]*O2B[1];
38
39 if (test1*test2>0){
40 return SeparateEnum;
41 }
42
43 test3=n1[0]*O1C[0]+n1[1]*O1C[1];
44 test4=n1[0]*O1D[0]+n1[1]*O1D[1];
45
46 if (test3*test4>0){
47 return SeparateEnum;
48 }
49
50 /*If colinear: */
51 det=n1[0]*n2[1]-n2[0]*y1[1];
52
53 if(test1*test2==0 && test3*test4==0 && det==0){
54
55 //projection on the axis O1O2
56 O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
57 O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
58
59 pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
60 pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
61 pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
62 pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
63
64 //test if one point is included in the other segment (->intersects=true)
65 if ((pO1C-pO1A)*(pO1D-pO1A)<0){
66 alpha=0; beta=0;
67 *palpha=alpha;*pbeta=beta;
68 return ColinearEnum;
69 }
70 if ((pO1C-pO1B)*(pO1D-pO1B)<0){
71 alpha=0; beta=0;
72 *palpha=alpha;*pbeta=beta;
73 return ColinearEnum;
74 }
75 if ((pO1A-pO1C)*(pO1B-pO1C)<0){
76 alpha=0; beta=0;
77 *palpha=alpha;*pbeta=beta;
78 return ColinearEnum;
79 }
80 if ((pO1A-pO1D)*(pO1B-pO1D)<0){
81 alpha=0; beta=0;
82 *palpha=alpha;*pbeta=beta;
83 return ColinearEnum;
84 }
85
86 //test if the 2 segments have the same middle (->intersects=true)
87 if (O2O1==0){
88 alpha=0; beta=0;
89 *palpha=alpha;*pbeta=beta;
90 return ColinearEnum;
91 }
92
93 //if we are here, both segments are colinear, but do not interset:
94 alpha=-1; beta=-1;
95 *palpha=alpha;*pbeta=beta;
96 return SeparateEnum;
97 }
98
99 /*if we are here, both segments intersect. Determine where in the segment coordinate
100 * system: */
101 beta=-1;
102 alpha=-(xA*yB-xC*yB+yC*xB-yC*xA+xC*yA-yA*xB)/(-xD*yB+xD*yA+xC*yB-xC*yA-yD*xA+yD*xB+yC*xA-yC*xB); //from intersect.m in formal calculus
103
104 *palpha=alpha;*pbeta=beta;
105 return IntersectEnum;
106}
Note: See TracBrowser for help on using the repository browser.