Changeset 1580


Ignore:
Timestamp:
09/19/18 14:28:26 (3 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:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/extern/remap/src/intersection_ym.cpp

    r849 r1580  
    88#include <stdlib.h> 
    99#include <limits> 
     10#include <array> 
     11#include <cstdint>  
     12#include "earcut.hpp" 
     13#include <fstream>  
     14 
    1015 
    1116#define epsilon 1e-3  // epsilon distance ratio over side lenght for approximate small circle by great circle 
     
    1621using namespace std; 
    1722using namespace ClipperLib ; 
     23         
    1824 
    1925double intersect_ym(Elt *a, Elt *b) 
    2026{ 
     27 
     28  using N = uint32_t; 
     29  using Point = array<double, 2>; 
     30  vector<Point> vect_points; 
     31  vector< vector<Point> > polyline; 
    2132 
    2233// transform small circle into piece of great circle if necessary 
     
    4556  double cos_alpha; 
    4657 
     58  /// vector<p2t::Point*> polyline; 
    4759  for(int n=0; n<na;n++) 
    4860  { 
     
    5163    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 
    5264    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 
    53   } 
     65 
     66    vect_points.push_back( array<double, 2>() ); 
     67    vect_points[n][0] = a_gno[n].x; 
     68    vect_points[n][1] = a_gno[n].y; 
     69 
     70  } 
     71 
     72  polyline.push_back(vect_points); 
     73  vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 
     74   
     75  double area_a_gno=0 ; 
     76  for(int i=0;i<indices_a_gno.size()/3;++i) 
     77    { 
     78      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 
     79      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 
     80      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 
     81      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
     82    } 
     83 
     84  vect_points.clear(); 
     85  polyline.clear(); 
     86  indices_a_gno.clear(); 
     87 
     88  
    5489 
    5590  for(int n=0; n<nb;n++) 
     
    5994    b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 
    6095    b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1 
    61   } 
    62  
     96 
     97    vect_points.push_back( array<double, 2>() ); 
     98    vect_points[n][0] = b_gno[n].x; 
     99    vect_points[n][1] = b_gno[n].y; 
     100  } 
     101 
     102 
     103  polyline.push_back(vect_points); 
     104  vector<N> indices_b_gno = mapbox::earcut<N>(polyline); 
     105 
     106  double area_b_gno=0 ; 
     107  for(int i=0;i<indices_b_gno.size()/3;++i) 
     108    { 
     109      Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ; 
     110      Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ; 
     111      Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ; 
     112      area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
     113    } 
     114 
     115  vect_points.clear(); 
     116  polyline.clear(); 
     117  indices_b_gno.clear(); 
    63118 
    64119 
     
    109164  clip.AddPaths(dst, ptClip, true); 
    110165  clip.Execute(ctIntersection, intersection); 
    111  
     166  
    112167  double area=0 ; 
    113168 
    114169  for(int ni=0;ni<intersection.size(); ni++) 
    115170  { 
    116     // go back into real coordinate on the sphere 
    117171    Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 
    118172    for(int n=0; n < intersection[ni].size(); n++) 
    119173    { 
    120       double x=intersection[ni][n].X/xscale+xoffset ; 
    121       double y=intersection[ni][n].Y/yscale+yoffset ; 
    122  
    123       intersectPolygon[n]=Ox*x+Oy*y+Oz ; 
    124       intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 
    125     } 
    126  
    127 // remove redondants vertex 
    128     int nv=0 ; 
     174      intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ; 
     175      intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ; 
     176    } 
     177     
     178 
     179    int nv=0; 
     180 
    129181    for(int n=0; n < intersection[ni].size(); n++) 
    130182    { 
    131       if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 
     183       double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ; 
     184       double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ; 
     185      
     186       if (dx*dx+dy*dy>fusion_vertex*fusion_vertex) 
     187       { 
     188          intersectPolygon[nv]=intersectPolygon[n] ; 
     189 
     190          vect_points.push_back( array<double, 2>() ); 
     191          vect_points[nv][0] = intersectPolygon[n].x; 
     192          vect_points[nv][1] = intersectPolygon[n].y; 
     193 
     194          nv++ ; 
     195       } 
     196       
     197 
     198    } 
     199 
     200    polyline.push_back(vect_points); 
     201    vect_points.clear(); 
     202 
     203    if (nv>2) 
     204    { 
     205  
     206      vector<N> indices = mapbox::earcut<N>(polyline); 
     207 
     208      double area2=0 ; 
     209      for(int i=0;i<indices.size()/3;++i) 
    132210      { 
    133         intersectPolygon[nv]=intersectPolygon[n] ; 
    134         nv++ ; 
     211          Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ; 
     212          Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ; 
     213          Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ; 
     214          area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
    135215      } 
    136     } 
    137  
    138  
    139     if (nv>2) 
    140     { 
     216 
     217      polyline.clear(); 
     218 
     219      for(int n=0; n < nv; n++) 
     220      { 
     221        intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz; 
     222      } 
     223 
     224 
    141225//     assign intersection to source and destination polygons 
    142226       Polyg *is = new Polyg; 
    143227       is->x = exact_barycentre(intersectPolygon,nv); 
    144        is->area = polygonarea(intersectPolygon,nv) ; 
     228//       is->area = polygonarea(intersectPolygon,nv) ; 
     229       is->area = area2 ; 
     230 
    145231//        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 
    146232       if (is->area==0.) delete is ; 
     
    161247  delete[] b_gno ; 
    162248  return area ; 
     249 
    163250} 
     251 
     252 
    164253 
    165254void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates) 
Note: See TracChangeset for help on using the changeset viewer.