Changeset 1580
 Timestamp:
 09/19/18 14:28:26 (3 years ago)
 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 8 8 #include <stdlib.h> 9 9 #include <limits> 10 #include <array> 11 #include <cstdint> 12 #include "earcut.hpp" 13 #include <fstream> 14 10 15 11 16 #define epsilon 1e3 // epsilon distance ratio over side lenght for approximate small circle by great circle … … 16 21 using namespace std; 17 22 using namespace ClipperLib ; 23 18 24 19 25 double intersect_ym(Elt *a, Elt *b) 20 26 { 27 28 using N = uint32_t; 29 using Point = array<double, 2>; 30 vector<Point> vect_points; 31 vector< vector<Point> > polyline; 21 32 22 33 // transform small circle into piece of great circle if necessary … … 45 56 double cos_alpha; 46 57 58 /// vector<p2t::Point*> polyline; 47 59 for(int n=0; n<na;n++) 48 60 { … … 51 63 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 52 64 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 54 89 55 90 for(int n=0; n<nb;n++) … … 59 94 b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 60 95 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(); 63 118 64 119 … … 109 164 clip.AddPaths(dst, ptClip, true); 110 165 clip.Execute(ctIntersection, intersection); 111 166 112 167 double area=0 ; 113 168 114 169 for(int ni=0;ni<intersection.size(); ni++) 115 170 { 116 // go back into real coordinate on the sphere117 171 Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 118 172 for(int n=0; n < intersection[ni].size(); n++) 119 173 { 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 129 181 for(int n=0; n < intersection[ni].size(); n++) 130 182 { 131 if (norm(intersectPolygon[n]intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 183 double dx=intersectPolygon[n].xintersectPolygon[(n+1)%intersection[ni].size()].x ; 184 double dy=intersectPolygon[n].yintersectPolygon[(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) 132 210 { 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))) ; 135 215 } 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 141 225 // assign intersection to source and destination polygons 142 226 Polyg *is = new Polyg; 143 227 is>x = exact_barycentre(intersectPolygon,nv); 144 is>area = polygonarea(intersectPolygon,nv) ; 228 // is>area = polygonarea(intersectPolygon,nv) ; 229 is>area = area2 ; 230 145 231 // if (is>area < 1e12) cout<<"Small intersection : "<<is>area<<endl ; 146 232 if (is>area==0.) delete is ; … … 161 247 delete[] b_gno ; 162 248 return area ; 249 163 250 } 251 252 164 253 165 254 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates)
Note: See TracChangeset
for help on using the changeset viewer.