Changeset 1580 for XIOS/trunk

Timestamp:

09/19/18 14:28:26 (6 years ago)

Author:

ymipsl

Message:

Computation of the area of a polygon on the sphere :
The total area was obtain by triangulation of vertex with the barycenter of the polygons. for some pathological polygones, baricenter was outside or triangle side was cutting polygons side, leading to numerical imprecision.

Now the polygon is triangulate with by external lib (earcut), and correct area is computed in any case.

YM

Location:

XIOS/trunk/extern/remap/src

Files:

• earcut.hpp (added)
• intersection_ym.cpp (modified) (7 diffs)

XIOS/trunk/extern/remap/src/intersection_ym.cpp

@@ -8 +8 @@
 #include <stdlib.h>
 #include <limits>
+#include <array>
+#include <cstdint> 
+#include "earcut.hpp"
+#include <fstream> 
+
 
 #define epsilon 1e-3  // epsilon distance ratio over side lenght for approximate small circle by great circle
@@ -16 +21 @@
 using namespace std;
 using namespace ClipperLib ;
+        
 
 double intersect_ym(Elt *a, Elt *b)
 {
+
+  using N = uint32_t;
+  using Point = array<double, 2>;
+  vector<Point> vect_points;
+  vector< vector<Point> > polyline;
 
 // transform small circle into piece of great circle if necessary
@@ -45 +56 @@
   double cos_alpha;
 
+  /// vector<p2t::Point*> polyline;
   for(int n=0; n<na;n++)
   {
@@ -51 +63 @@
     a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
     a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = a_gno[n].x;
+    vect_points[n][1] = a_gno[n].y;
+
+  }
+
+  polyline.push_back(vect_points);
+  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
+  
+  double area_a_gno=0 ;
+  for(int i=0;i<indices_a_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
+      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_a_gno.clear();
+
+ 
 
   for(int n=0; n<nb;n++)
@@ -59 +94 @@
     b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ;
     b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1
-  }
-
+
+    vect_points.push_back( array<double, 2>() );
+    vect_points[n][0] = b_gno[n].x;
+    vect_points[n][1] = b_gno[n].y;
+  }
+
+
+  polyline.push_back(vect_points);
+  vector<N> indices_b_gno = mapbox::earcut<N>(polyline);
+
+  double area_b_gno=0 ;
+  for(int i=0;i<indices_b_gno.size()/3;++i)
+    {
+      Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ;
+      Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ;
+      Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ;
+      area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+  vect_points.clear();
+  polyline.clear();
+  indices_b_gno.clear();
 
 
@@ -109 +164 @@
   clip.AddPaths(dst, ptClip, true);
   clip.Execute(ctIntersection, intersection);
-
+ 
   double area=0 ;
 
   for(int ni=0;ni<intersection.size(); ni++)
   {
-    // go back into real coordinate on the sphere
     Coord* intersectPolygon=new Coord[intersection[ni].size()] ;
     for(int n=0; n < intersection[ni].size(); n++)
     {
-      double x=intersection[ni][n].X/xscale+xoffset ;
-      double y=intersection[ni][n].Y/yscale+yoffset ;
-
-      intersectPolygon[n]=Ox*x+Oy*y+Oz ;
-      intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ;
-    }
-
-// remove redondants vertex
-    int nv=0 ;
+      intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ;
+      intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ;
+    }
+    
+
+    int nv=0;
+
     for(int n=0; n < intersection[ni].size(); n++)
     {
-      if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex)
+       double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ;
+       double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ;
+     
+       if (dx*dx+dy*dy>fusion_vertex*fusion_vertex)
+       {
+          intersectPolygon[nv]=intersectPolygon[n] ;
+
+          vect_points.push_back( array<double, 2>() );
+          vect_points[nv][0] = intersectPolygon[n].x;
+          vect_points[nv][1] = intersectPolygon[n].y;
+
+          nv++ ;
+       }
+      
+
+    }
+
+    polyline.push_back(vect_points);
+    vect_points.clear();
+
+    if (nv>2)
+    {
+ 
+      vector<N> indices = mapbox::earcut<N>(polyline);
+
+      double area2=0 ;
+      for(int i=0;i<indices.size()/3;++i)
       {
-        intersectPolygon[nv]=intersectPolygon[n] ;
-        nv++ ;
+          Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ;
+          Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ;
+          Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ;
+          area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
       }
-    }
-
-
-    if (nv>2)
-    {
+
+      polyline.clear();
+
+      for(int n=0; n < nv; n++)
+      {
+        intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz;
+      }
+
+
 //     assign intersection to source and destination polygons
        Polyg *is = new Polyg;
        is->x = exact_barycentre(intersectPolygon,nv);
-       is->area = polygonarea(intersectPolygon,nv) ;
+//       is->area = polygonarea(intersectPolygon,nv) ;
+       is->area = area2 ;
+
 //        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ;
        if (is->area==0.) delete is ;
@@ -161 +247 @@
   delete[] b_gno ;
   return area ;
+
 }
+
+
 
 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates)