Changeset 22758 for issm/trunk/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

Timestamp:

05/10/18 10:24:27 (7 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 22757

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/modules (modified) (1 prop)
• src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (modified) (6 diffs)

issm/trunk/src/c

@@ -20 +20 @@
 kriging
 issm_slr
+issm_ocean
+lnb_param.mod
+lovenb_sub.mod
+model.mod
+util.mod

@@ -20 +20 @@
 kriging
 issm_slr
+issm_ocean
+lnb_param.mod
+lovenb_sub.mod
+model.mod
+util.mod

issm/trunk/src/c/modules

@@ -1 +1 @@
 .deps
+InputControlUpdatex
+InputScalex
+TriaSearchx

@@ -1 +1 @@
 .deps
+InputControlUpdatex
+InputScalex
+TriaSearchx

issm/trunk/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

@@ -9 +9 @@
 #endif
 
-/*Include files: {{{*/
+/*Include files*/
 #include "./ExpToLevelSetx.h"
 double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0);
-void ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
-/*}}}*/
+void   ContourToLevelSet(double* distance,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i10);
 
 void* ExpToLevelSetxt(void* vpthread_handle){
@@ -29 +28 @@
 
         /*recover parameters :*/
-        Contours* contours  = gate->contours;
-        int       nods      = gate->nods;
-        double   *distance    = gate->distance;
-        double   *x         = gate->x;
-        double   *y         = gate->y;
+        Contours* contours = gate->contours;
+        int       nods     = gate->nods;
+        double   *distance = gate->distance;
+        double   *x        = gate->x;
+        double   *y        = gate->y;
 
         /*distribute indices across threads :*/
@@ -39 +38 @@
 
         /*Loop through all contours: */
-        for (i=0;i<contours->Size();i++){
+        for(i=0;i<contours->Size();i++){
                 Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
                 ContourToLevelSet(distance,contour->x,contour->y,contour->nods,x,y,i0,i1);
@@ -48 +47 @@
 
 void ContourToLevelSet(double* dist,double* contourx, double* contoury, int contournods, double* x, double* y, int i0, int i1){/*{{{*/
-        int i,j;
+
         double x0,y0;
         double x1,y1;
@@ -54 +53 @@
         double mind;
         
-        for(i=i0;i<i1;i++){
+        for(int i=i0;i<i1;i++){
+
+      /*Get current point*/
                 x0=x[i]; y0=y[i];
 
                 /*Figure out distance from (x0,y0) to contour: */
                 mind=1e+50;
-                for (j=0;j<contournods-1;j++){
-                        x1=contourx[j]; y1=contoury[j];
+                for(int j=0;j<contournods-1;j++){
+         /*Get distance from current segment*/
+                        x1=contourx[j];   y1=contoury[j];
                         x2=contourx[j+1]; y2=contoury[j+1];
                         mind=min(mind,minimum_distance(x1,y1,x2,y2,x0,y0));
@@ -66 +68 @@
                 dist[i]=min(dist[i],mind);
         }
+}
+double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
+        /* Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0)
+         * We use the following notations:
+         * segment: v=(x1,y1), w=(x2,y2)
+         * point:   p=(x0,y0) 
+         */
 
-}
-double ddistance(double x1,double y1,double x2,double y2){
-        return sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-}
-double ddot(double x1, double y1, double x2, double y2){
-        return x1*x2+y1*y2;
-}
+   /*Get segment length square (avoid sqrt) |w-v|^2*/
+        double l2 = pow(x2-x1,2)+pow(y2-y1,2); 
 
-bool isPointLeftOfRay(double x, double y, double raySx, double raySy, double rayEx, double rayEy) {
-          return (y-raySy)*(rayEx-raySx)
-                    >      (x-raySx)*(rayEy-raySy); 
-}
+   /*segment is single point: v == w*/
+        if(l2 == 0.) return sqrt(pow(x1-x0,2)+pow(y1-y0,2));
 
-double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
+        /*Consider the line extending the segment, parameterized as v + t (w - v).
+        We find projection of point p onto the line.  It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+        double t = ((x0-x1)*(x2-x1) + (y0-y1)*(y2-y1)) / l2;
+        if(t < 0.){
+      /*Beyond the 'v' end of the segment*/
+      return sqrt(pow(x1-x0,2)+pow(y1-y0,2));
+   }
+        else if(t > 1.){
+      /*Beyond the 'w' end of the segment*/
+      return sqrt(pow(x2-x0,2)+pow(y2-y0,2));
+   }
         
-        // Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)
-        double projectionx; 
-        double projectiony; 
-        double l2;
-        double t;
-
-        l2 = pow(x2-x1,2)+pow(y2-y1,2); // i.e. |w-v|^2 -  avoid a sqrt
-
-        if (l2 == 0.0) return ddistance(x0,y0, x1,y1); // v == w case
-        // Consider the line extending the segment, parameterized as v + t (w - v).
-        //         // We find projection of point p onto the line. 
-        //           // It falls where t = [(p-v) . (w-v)] / |w-v|^2
-        t = ddot(x0-x1,y0-y1, x2-x1, y2-y1) / l2;
-        if (t < 0.0) return ddistance(x0,y0, x1, y1);       // Beyond the 'v' end of the segment
-        else if (t > 1.0) return ddistance(x0,y0, x2,y2);  // Beyond the 'w' end of the segment
-        
-        projectionx= x1 + t* (x2-x1);  // Projection falls on the segment
-        projectiony= y1 + t* (y2-y1);
-        return ddistance(x0, y0, projectionx, projectiony);
+   /*Projection falls on segment*/
+        double projx= x1 + t* (x2-x1);
+        double projy= y1 + t* (y2-y1);
+        return sqrt(pow(projx-x0,2)+pow(projy-y0,2));
 }
 /*}}}*/

@@ -20 +20 @@
 kriging
 issm_slr
+issm_ocean
+lnb_param.mod
+lovenb_sub.mod
+model.mod
+util.mod

@@ -1 +1 @@
 .deps
+InputControlUpdatex
+InputScalex
+TriaSearchx