Changeset 22456
- Timestamp:
- 02/23/18 11:52:37 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp
r21931 r22456 9 9 #endif 10 10 11 /*Include files : {{{*/11 /*Include files*/ 12 12 #include "./ExpToLevelSetx.h" 13 13 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0); 14 void ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10); 15 /*}}}*/ 14 void ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10); 16 15 17 16 void* ExpToLevelSetxt(void* vpthread_handle){ … … 29 28 30 29 /*recover parameters :*/ 31 Contours* contours 32 int nods 33 double *distance 34 double *x 35 double *y 30 Contours* contours = gate->contours; 31 int nods = gate->nods; 32 double *distance = gate->distance; 33 double *x = gate->x; 34 double *y = gate->y; 36 35 37 36 /*distribute indices across threads :*/ … … 39 38 40 39 /*Loop through all contours: */ 41 for 40 for(i=0;i<contours->Size();i++){ 42 41 Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i); 43 42 ContourToLevelSet(distance,contour->x,contour->y,contour->nods,x,y,i0,i1); … … 48 47 49 48 void ContourToLevelSet(double* dist,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i1){/*{{{*/ 50 int i,j;51 49 double x0,y0; 52 50 double x1,y1; … … 54 52 double mind; 55 53 56 for(i=i0;i<i1;i++){ 54 for(int i=i0;i<i1;i++){ 55 56 /*Get current point*/ 57 57 x0=x[i]; y0=y[i]; 58 58 59 59 /*Figure out distance from (x0,y0) to contour: */ 60 60 mind=1e+50; 61 for (j=0;j<contournods-1;j++){ 62 x1=contourx[j]; y1=contoury[j]; 61 for(int j=0;j<contournods-1;j++){ 62 /*Get distance from current segment*/ 63 x1=contourx[j]; y1=contoury[j]; 63 64 x2=contourx[j+1]; y2=contoury[j+1]; 64 65 mind=min(mind,minimum_distance(x1,y1,x2,y2,x0,y0)); … … 66 67 dist[i]=min(dist[i],mind); 67 68 } 68 69 69 } 70 70 double ddistance(double x1,double y1,double x2,double y2){ return sqrt(pow(x2-x1,2)+pow(y2-y1,2)); } 71 71 double ddot(double x1, double y1, double x2, double y2){ return x1*x2+y1*y2; } 72 72 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){ 73 /* Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)*/ 74 75 /*Get segment length square (avoid sqrt) |w-v|^2*/ 76 double l2 = pow(x2-x1,2)+pow(y2-y1,2); 77 78 /*segment is single point: v == w*/ 79 if(l2 == 0.) return ddistance(x0,y0,x1,y1); 80 81 /*Consider the line extending the segment, parameterized as v + t (w - v). 82 We find projection of point p onto the line. 83 It falls where t = [(p-v) . (w-v)] / |w-v|^2*/ 84 double t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2; 85 if(t < 0.){ 86 /*Beyond the 'v' end of the segment*/ 87 return ddistance(x0,y0, x1, y1); 88 } 89 else if(t > 1.){ 90 /*Beyond the 'w' end of the segment*/ 91 return ddistance(x0,y0, x2,y2); 92 } 73 93 74 // Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0) 75 double projectionx; 76 double projectiony; 77 double l2; 78 double t; 79 80 l2 = pow(x2-x1,2)+pow(y2-y1,2); // i.e. |w-v|^2 - avoid a sqrt 81 82 if (l2 == 0.0) return ddistance(x0,y0, x1,y1); // v == w case 83 // Consider the line extending the segment, parameterized as v + t (w - v). 84 // // We find projection of point p onto the line. 85 // // It falls where t = [(p-v) . (w-v)] / |w-v|^2 86 t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2; 87 if (t < 0.0) return ddistance(x0,y0, x1, y1); // Beyond the 'v' end of the segment 88 else if (t > 1.0) return ddistance(x0,y0, x2,y2); // Beyond the 'w' end of the segment 89 90 projectionx= x1 + t* (x2-x1); // Projection falls on the segment 91 projectiony= y1 + t* (y2-y1); 94 /*Projection falls on segment!*/ 95 double projectionx= x1 + t* (x2-x1); 96 double projectiony= y1 + t* (y2-y1); 92 97 return ddistance(x0, y0, projectionx, projectiony); 93 98 }
Note:
See TracChangeset
for help on using the changeset viewer.