- Timestamp:
- 01/23/19 10:31:44 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp
r1335 r1642 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 1e-3 // epsilon distance ratio over side lenght for approximate small circle by great circle … … 14 19 namespace sphereRemap { 15 20 16 extern CRemapGrid srcGrid;17 #pragma omp threadprivate(srcGrid)18 19 extern CRemapGrid tgtGrid;20 #pragma omp threadprivate(tgtGrid)21 22 23 21 using namespace std; 24 22 using namespace ClipperLib ; 23 25 24 26 25 double intersect_ym(Elt *a, Elt *b) 27 26 { 27 28 using N = uint32_t; 29 using Point = array<double, 2>; 30 vector<Point> vect_points; 31 vector< vector<Point> > polyline; 28 32 29 33 // transform small circle into piece of great circle if necessary … … 52 56 double cos_alpha; 53 57 58 /// vector<p2t::Point*> polyline; 54 59 for(int n=0; n<na;n++) 55 60 { … … 58 63 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 59 64 a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 60 } 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 61 89 62 90 for(int n=0; n<nb;n++) … … 66 94 b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 67 95 b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1 68 } 69 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(); 70 118 71 119 … … 116 164 clip.AddPaths(dst, ptClip, true); 117 165 clip.Execute(ctIntersection, intersection); 118 166 119 167 double area=0 ; 120 168 121 169 for(int ni=0;ni<intersection.size(); ni++) 122 170 { 123 // go back into real coordinate on the sphere124 171 Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 125 172 for(int n=0; n < intersection[ni].size(); n++) 126 173 { 127 double x=intersection[ni][n].X/xscale+xoffset ; 128 double y=intersection[ni][n].Y/yscale+yoffset ; 129 130 intersectPolygon[n]=Ox*x+Oy*y+Oz ; 131 intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 132 } 133 134 // remove redondants vertex 135 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 136 181 for(int n=0; n < intersection[ni].size(); n++) 137 182 { 138 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) 139 210 { 140 intersectPolygon[nv]=intersectPolygon[n] ; 141 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))) ; 142 215 } 143 } 144 145 146 if (nv>2) 147 { 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 intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 223 } 224 225 148 226 // assign intersection to source and destination polygons 149 227 Polyg *is = new Polyg; 150 228 is->x = exact_barycentre(intersectPolygon,nv); 151 is->area = polygonarea(intersectPolygon,nv) ; 229 // is->area = polygonarea(intersectPolygon,nv) ; 230 is->area = area2 ; 231 152 232 // if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 153 233 if (is->area==0.) delete is ; … … 168 248 delete[] b_gno ; 169 249 return area ; 250 170 251 } 252 253 171 254 172 255 void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates)
Note: See TracChangeset
for help on using the changeset viewer.