1 #include "GeomCoords.h"
2
3 namespace Upp {
4
5 #define LLOG(x) // RLOG(x)
6
DegreeStep(double min_step)7 double DegreeStep(double min_step)
8 {
9 static const double step[] =
10 {
11 // 90,
12 // 45,
13 // 30,
14 15,
15 10,
16 5,
17 2,
18 1,
19 30 / 60.0,
20 20 / 60.0,
21 15 / 60.0,
22 10 / 60.0,
23 5 / 60.0,
24 2 / 60.0,
25 1 / 60.0,
26 30 / 3600.0,
27 20 / 3600.0,
28 15 / 3600.0,
29 10 / 3600.0,
30 5 / 3600.0,
31 2 / 3600.0,
32 1 / 3600.0,
33 };
34 const double *p = step;
35 while(p < step + __countof(step) - 1 && p[1] >= min_step)
36 p++;
37 return *p;
38 }
39
DegreeMask(double start_angle,double end_angle)40 int DegreeMask(double start_angle, double end_angle)
41 {
42 if(end_angle < start_angle)
43 Swap(start_angle, end_angle);
44 double len = end_angle - start_angle;
45 if(len >= 360)
46 return 0xFF;
47 if(len <= 0)
48 return 0;
49 int sx = fceil(start_angle / 90);
50 int nx = 1 << ((ffloor(end_angle / 90) - sx + 1) & 7);
51 int em = (nx - 1) << (sx & 3);
52 int qm = (2 * nx - 1) << ((sx - 1) & 3);
53 return ((em | (em >> 4)) & 0x0F) | ((qm | (qm << 4)) & 0xF0);
54 }
55
DegreeToExtent(const Rectf & lonlat)56 Rectf DegreeToExtent(const Rectf& lonlat)
57 {
58 Rectf out;
59 int mask = DegreeMask(lonlat.left, lonlat.right);
60 double lrad = lonlat.left * DEGRAD, rrad = lonlat.right * DEGRAD;
61 double a, b;
62 a = sin(lrad);
63 b = sin(rrad);
64 if(a > b) Swap(a, b);
65 out.left = (mask & AM_E3 ? -lonlat.bottom : a * (a >= 0 ? lonlat.top : lonlat.bottom));
66 out.right = (mask & AM_E1 ? +lonlat.bottom : b * (b >= 0 ? lonlat.bottom : lonlat.top));
67 a = -cos(lrad);
68 b = -cos(rrad);
69 if(a > b) Swap(a, b);
70 out.top = (mask & AM_E0 ? -lonlat.bottom : a * (a >= 0 ? lonlat.top : lonlat.bottom));
71 out.bottom = (mask & AM_E2 ? +lonlat.bottom : b * (b >= 0 ? lonlat.bottom : lonlat.top));
72 return out;
73 }
74
ExtentToDegree(const Rectf & xy)75 Rectf ExtentToDegree(const Rectf& xy)
76 {
77 double mineps, maxeps, minrho, maxrho;
78 if(xy.left >= 0)
79 {
80 mineps = atan2(xy.left, -xy.top);
81 maxeps = atan2(xy.left, -xy.bottom);
82 minrho = hypot(xy.left, xy.top > 0 ? xy.top : xy.bottom < 0 ? xy.bottom : 0);
83 maxrho = hypot(xy.right, max(fabs(xy.top), fabs(xy.bottom)));
84 }
85 else if(xy.right <= 0)
86 {
87 mineps = atan2(xy.right, -xy.bottom);
88 maxeps = atan2(xy.right, -xy.top);
89 minrho = hypot(xy.right, xy.top > 0 ? xy.top : xy.bottom < 0 ? xy.bottom : 0);
90 maxrho = hypot(xy.left, max(fabs(xy.top), fabs(xy.bottom)));
91 }
92 else if(xy.top >= 0)
93 {
94 mineps = atan2(xy.right, -xy.top);
95 maxeps = atan2(xy.left, -xy.top);
96 minrho = hypot(xy.top, xy.left > 0 ? xy.left : xy.right < 0 ? xy.right : 0);
97 maxrho = hypot(xy.bottom, max(fabs(xy.left), fabs(xy.right)));
98 }
99 else if(xy.bottom <= 0)
100 {
101 mineps = atan2(xy.left, -xy.bottom);
102 maxeps = atan2(xy.right, -xy.bottom);
103 minrho = hypot(xy.bottom, xy.left > 0 ? xy.left : xy.right < 0 ? xy.right : 0);
104 maxrho = hypot(xy.top, max(fabs(xy.left), fabs(xy.right)));
105 }
106 else
107 {
108 mineps = -M_PI;
109 maxeps = +M_PI;
110 minrho = 0;
111 maxrho = hypot(max(-xy.left, xy.right), max(-xy.top, xy.bottom));
112 }
113 return Rectf(mineps / DEGRAD, minrho, maxeps / DEGRAD, maxrho);
114 }
115
GisLengthDecimals(double pixel_len)116 int GisLengthDecimals(double pixel_len)
117 {
118 return minmax<int>(1 - ilog10(pixel_len), 0, 8);
119 }
120
DegreeDecimals(double pixel_angle)121 int DegreeDecimals(double pixel_angle)
122 {
123 pixel_angle = fabs(pixel_angle);
124 if(pixel_angle >= 1)
125 return -2;
126 if(pixel_angle >= 1 / 60.0)
127 return -1;
128 return minmax<int>(-ilog10(pixel_angle / (1 / 3600.0)), 0, 3);
129 }
130
FormatDegree(double d,int decimals,bool spaces)131 String FormatDegree(double d, int decimals, bool spaces)
132 {
133 if(IsNull(d))
134 return Null;
135 d = modulo(d + 180, 360) - 180;
136 char sign = (d < 0 ? '-' : '+');
137 if(d < 0) d = -d;
138 int deg = ffloor(d);
139 String cd = ToCharset(CHARSET_DEFAULT, "%c%d°", CHARSET_UTF8);
140 if(decimals <= -2)
141 return NFormat(cd, sign, deg);
142 d = (d - deg) * 60;
143 int min = ffloor(d);
144 if(decimals <= -1)
145 return NFormat(cd + (spaces ? " %02d\'" : "%02d\'"), sign, deg, min);
146 d = (d - min) * 60;
147 String sec = FormatDoubleFix(d, decimals);
148 if(!IsDigit(sec[1]))
149 sec.Insert(0, '0');
150 return NFormat(cd + (spaces ? " %02d\' %s\"" : "%02d\'%s\""), sign, deg, min, sec);
151 }
152
ScanDegree(const char * p)153 Value ScanDegree(const char *p)
154 {
155 int deg = ScanInt(p, &p);
156 int min = 0;
157 double sec = 0;
158 if(IsNull(deg))
159 return Null;
160 if(deg < -360 || deg > 360)
161 return ErrorValue(NFormat("Neplatný poèet úhlových stupòù: %d", deg));
162 while(*p && !IsDigit(*p))
163 p++;
164 if(*p)
165 {
166 min = ScanInt(p, &p);
167 if(min < 0 || min >= 60)
168 return ErrorValue(NFormat("Neplatný poèet úhlových minut: %d", min));
169 while(*p && !IsDigit(*p))
170 p++;
171 if(*p)
172 {
173 sec = ScanDouble(p);
174 if(IsNull(sec) || sec < 0 || sec > 60)
175 return ErrorValue(NFormat("Neplatný poèet úhlových sekund: %d", sec));
176 }
177 }
178 return ((sec / 60.0 + min) / 60.0 + tabs(deg)) * (deg >= 0 ? 1 : -1);
179 }
180
ConvertDegree(int decimals,bool not_null,double min,double max)181 ConvertDegree::ConvertDegree(int decimals, bool not_null, double min, double max)
182 : decimals(decimals), not_null(not_null), min(min), max(max)
183 {}
184
~ConvertDegree()185 ConvertDegree::~ConvertDegree() {}
186
Format(const Value & v) const187 Value ConvertDegree::Format(const Value& v) const
188 {
189 if(IsNull(v) || !IsNumber(v))
190 return v;
191 return FormatDegree(v, decimals);
192 }
193
Scan(const Value & v) const194 Value ConvertDegree::Scan(const Value& v) const
195 {
196 if(IsNull(v) || !IsString(v))
197 return v;
198 double d = ScanDegree((String)v);
199 if(IsNull(d))
200 {
201 if(not_null)
202 return ErrorValue("Hodnota nesmí být prázdná.");
203 return Null;
204 }
205 if(!IsNull(min) && d < min)
206 return ErrorValue(NFormat("Zadaný úhel je menší než povolená dolní mez, %s.", FormatDegree(min, 0)));
207 if(!IsNull(max) && d > max)
208 return ErrorValue(NFormat("Zadaný úhel je vìtší než povolená horní mez, %s.", FormatDegree(max, 0)));
209 return d;
210 }
211
Filter(int c) const212 int ConvertDegree::Filter(int c) const
213 {
214 return IsDigit(c) || c == ':' || c == '\'' || c == '\"' || c == ',' || c == '.' || c == '/'
215 || c == '+' || c == '-' ? c : 0;
216 }
217
Dump(double xmin,double xmax,int steps) const218 void GisFunction::Dump(double xmin, double xmax, int steps) const
219 {
220 double dy = 0;
221 for(int i = 0; i <= steps; i++)
222 {
223 double x = xmin + (xmax - xmin) * i / steps;
224 double y = Get(x);
225 RLOG(NFormat("%10nl\t%10nl\t%10nl", x, y, y - dy));
226 dy = y;
227 }
228 }
229
GisInverse(double xmin_,double xmax_,const GisFunction & fn_,int sections_,double epsilon_)230 GisInverse::GisInverse(double xmin_, double xmax_, const GisFunction& fn_, int sections_, double epsilon_)
231 : fn(fn_)
232 {
233 xstep = (xmax_ - xmin_) / sections_;
234 rawxmin = (xmin = xmin_) - xstep;
235 rawxmax = (xmax = xmax_) + xstep;
236 epsilon = epsilon_;
237 ymap.SetCount(sections_ + 3);
238 ymap[0] = fn(xmin - xstep);
239 ymap.Top() = fn(xmax + xstep);
240 rawymin = min(ymap[0], ymap.Top());
241 rawymax = max(ymap[0], ymap.Top());
242 ymin = ymax = ymap[1] = fn(xmin);
243 for(int i = 1; i <= sections_; i++)
244 {
245 double y = fn(xmin + i * xstep);
246 if(y < ymin) ymin = y;
247 if(y > ymax) ymax = y;
248 ymap[i + 1] = y;
249 }
250 if(ymin < rawymin) rawymin = ymin;
251 if(ymax > rawymax) rawymax = ymax;
252 ystep = (rawymax - rawymin) / sections_;
253 int prev = minmax(ffloor((ymap[0] - rawymin) / ystep), 0, sections_ + 1);
254 for(int i = 0; i < sections_ + 2; i++)
255 {
256 int next = minmax(ffloor((ymap[i + 1] - rawymin) / ystep), 0, sections_ + 1);
257 if(prev <= next)
258 while(prev < next)
259 accel.Add(prev++, i);
260 else
261 while(prev > next)
262 accel.Add(prev--, i);
263 accel.Add(prev = next, i);
264 }
265 }
266
Get(double y) const267 double GisInverse::Get(double y) const
268 {
269 for(int f = accel.Find(minmax<int>((int)((y - rawymin) / ystep), 0, accel.GetCount() - 1)); f >= 0;
270 f = accel.FindNext(f))
271 {
272 int sec = accel[f];
273 if(ymap[sec] == y)
274 return rawxmin + xstep * sec;
275 else if(ymap[sec + 1] == y)
276 return rawxmin + xstep * (sec + 1);
277 else if(ymap[sec] > y && ymap[sec + 1] < y)
278 {
279 double lx = rawxmin + xstep * sec, hx = lx + xstep;
280 // double ly = ymap[sec], hy = ymap[sec + 1];
281 while(hx - lx > epsilon)
282 {
283 /* double dh = hy - ly, mx, my;
284 if(fabs(dh) > epsilon)
285 mx = lx + (y - ly) * (hx - lx) / dh;
286 else
287 */
288 double mx = (lx + hx) / 2;
289 double my = fn(mx);
290 if(my > y)
291 {
292 lx = mx;
293 // ly = my;
294 }
295 else
296 {
297 hx = mx;
298 // hy = my;
299 }
300 }
301 return (lx + hx) / 2;
302 }
303 else if(ymap[sec] < y && ymap[sec + 1] > y)
304 {
305 double lx = rawxmin + xstep * sec, hx = lx + xstep;
306 double ly = ymap[sec], hy = ymap[sec + 1];
307 while(hx - lx > epsilon)
308 {
309 double dh = hy - ly;
310 double mx = (fabs(dh) > epsilon ? lx + (y - ly) * (hx - lx) / dh : (lx + hx) / 2);
311 if(mx - lx <= epsilon)
312 return mx;
313 if(hx - mx <= epsilon)
314 return mx;
315 double my = fn(mx);
316 if(my < y)
317 {
318 lx = mx;
319 ly = my;
320 }
321 else
322 {
323 hx = mx;
324 hy = my;
325 }
326 }
327 return (lx + hx) / 2;
328 }
329 }
330 return 0;
331 }
332
GisInverseDelta(double xmin,double xmax,const GisFunction & fn,int sections,double epsilon,int samples)333 String GisInverseDelta(double xmin, double xmax, const GisFunction& fn, int sections, double epsilon, int samples)
334 {
335 String out;
336 GisInverse inverse(xmin, xmax, fn, sections, epsilon);
337 int show_samples = min(samples, 1000);
338 double show_step = (xmax - xmin) / show_samples;
339 for(int i = 0; i < show_samples; i++)
340 {
341 double x = xmin + i * show_step;
342 double y = fn(x);
343 double ix = inverse(y);
344 out << NFormat("%15>10!nf %15>10!nf %15>10!nf %15>10!nf\n", x, y, ix, ix - x);
345 }
346 double step = (xmax - xmin) / samples;
347 double max_error = 0;
348 for(int i = 0; i < samples; i++)
349 {
350 double x = xmin + i * step;
351 double y = fn(x);
352 double ix = inverse(y);
353 max_error = max(max_error, fabs(ix - x));
354 }
355 return NFormat("max delta = %10n\n%s", max_error, out);
356 }
357
GisInverseTiming(double xmin,double xmax,const GisFunction & fn,int sections,double epsilon)358 String GisInverseTiming(double xmin, double xmax, const GisFunction& fn, int sections, double epsilon)
359 {
360 String out;
361 GisInverse inverse(xmin, xmax, fn, sections, epsilon);
362 Buffer<double> yval(1000);
363 for(int i = 0; i < 1000; i++)
364 yval[i] = inverse.GetYMin() + (inverse.GetYMax() - inverse.GetYMin()) * i / 999.0;
365 int start = msecs(), duration;
366 int count = 0;
367 do
368 {
369 count++;
370 // double x;
371 for(int i = 0; i < 1000; i++)
372 /*x = */inverse(yval[i]);
373 }
374 while((duration = msecs(start)) < 500);
375 double nsecs = duration * 1000.0 / double(count);
376 return NFormat("Function inverse: %4v nsecs", nsecs);
377 }
378
QuadraticLeastSquare(const double * vx,const double * vy,int count,double coef_out[3])379 void QuadraticLeastSquare(const double *vx, const double *vy, int count, double coef_out[3])
380 {
381 double left[3][3], right[3];
382 ZeroArray(left);
383 ZeroArray(right);
384 double xpow[3] = { 1 };
385 for(int s = 0; s < count; s++)
386 {
387 xpow[1] = vx[s];
388 xpow[2] = sqr(xpow[1]);
389 double sy = vy[s];
390 for(int y = 0; y < 3; y++)
391 {
392 double xy = xpow[y];
393 const double *xp = xpow;
394 double *dest = left[y];
395 for(int x = 0; x < 3; x++)
396 *dest++ += *xp++ * xy;
397 right[y] += sy * xy;
398 }
399 }
400 double D = Determinant(
401 left[0][0], left[0][1], left[0][2],
402 left[1][0], left[1][1], left[1][2],
403 left[2][0], left[2][1], left[2][2]);
404 double D1 = Determinant(
405 right[0], left[0][1], left[0][2],
406 right[1], left[1][1], left[1][2],
407 right[2], left[2][1], left[2][2]);
408 double DX = Determinant(
409 left[0][0], right[0], left[0][2],
410 left[1][0], right[1], left[1][2],
411 left[2][0], right[2], left[2][2]);
412 double DXX = Determinant(
413 left[0][0], left[0][1], right[0],
414 left[1][0], left[1][1], right[1],
415 left[2][0], left[2][1], right[2]);
416 coef_out[0] = D1 / D;
417 coef_out[1] = DX / D;
418 coef_out[2] = DXX / D;
419 }
420
Create(double xmin_,double xmax_,const GisFunction & fn,int buckets,int sections,int samples)421 void GisInterpolator::Create(double xmin_, double xmax_, const GisFunction& fn, int buckets, int sections, int samples)
422 {
423 samples &= ~1;
424 ASSERT(sections >= buckets);
425 bucket_index.SetCount(sections);
426 abc.SetCount(3 * sections);
427 xmin = xmin_;
428 xmax = xmax_;
429 step = (xmax_ - xmin_) / buckets;
430 limit = buckets - 1;
431 Vector<int> bucket_sections;
432 bucket_sections.SetCount(buckets, 1);
433 int section_count = buckets;
434 Vector<double> bucket_error;
435 bucket_error.SetCount(buckets);
436 Buffer<double> xsmpl(samples + 1), ysmpl(samples + 1);
437 for(int b = 0; b < buckets; b++)
438 {
439 double error = 0;
440 int nsec = bucket_sections[b];
441 double buck_begin = xmin + b * step;
442 double sec_step = step / nsec;
443 for(int s = 0; s < nsec; s++)
444 {
445 double sec_begin = buck_begin + s * sec_step;
446 double sample_step = sec_step / samples;
447 for(int m = 0; m <= samples; m++)
448 {
449 xsmpl[m] = m / double(samples);
450 ysmpl[m] = fn(sec_begin + sample_step * m);
451 }
452 double abc[3];
453 abc[0] = ysmpl[0];
454 abc[1] = 4 * ysmpl[samples >> 1] - ysmpl[samples] - 3 * ysmpl[0];
455 abc[2] = ysmpl[samples] - ysmpl[0] - abc[1];
456 // QuadraticLeastSquare(xsmpl, ysmpl, samples + 1, abc);
457 for(int m = 0; m <= samples; m++)
458 error = max(error, fabs(abc[0] + xsmpl[m] * (abc[1] + xsmpl[m] * abc[2]) - ysmpl[m]));
459 }
460 bucket_error[b] = error;
461 LLOG(NFormat("bucket[%d] (%4v .. %4v) error = %4ne", b, buck_begin, buck_begin + step, error));
462 }
463 Vector<int> order = GetSortOrder(bucket_error);
464 while(section_count + 1 < sections)
465 {
466 int worst = order.Top();
467 int add_sections = min(bucket_sections[worst], sections - section_count);
468 bucket_sections[worst] += add_sections;
469 section_count += add_sections;
470 double error = 0;
471 int nsec = bucket_sections[worst];
472 double buck_begin = xmin + worst * step;
473 double sec_step = step / nsec;
474 for(int s = 0; s < nsec; s++)
475 {
476 double sec_begin = buck_begin + s * sec_step;
477 double sample_step = sec_step / samples;
478 for(int m = 0; m <= samples; m++)
479 {
480 xsmpl[m] = m / double(samples);
481 ysmpl[m] = fn(sec_begin + sample_step * m);
482 }
483 double abc[3];
484 abc[0] = ysmpl[0];
485 abc[1] = 4 * ysmpl[samples >> 1] - ysmpl[samples] - 3 * ysmpl[0];
486 abc[2] = ysmpl[samples] - ysmpl[0] - abc[1];
487 // QuadraticLeastSquare(xsmpl, ysmpl, samples + 1, abc);
488 for(int m = 0; m <= samples; m++)
489 {
490 double abcval = abc[0] + xsmpl[m] * (abc[1] + xsmpl[m] * abc[2]);
491 double new_error = fabs(abcval - ysmpl[m]);
492 if(new_error > error)
493 {
494 // LLOG(NFormat("error at %10nf: abc = %10nf, y = %10nf, error = %10nf",
495 // sec_begin + sample_step * xsmpl[m], abcval, ysmpl[m], new_error));
496 error = new_error;
497 }
498 }
499 }
500 LLOG("total = " << section_count << ": bucket[" << worst << "] expand to " << nsec << ", error: "
501 << FormatDoubleExp(bucket_error[worst], 4) << " -> " << FormatDoubleExp(error, 4));
502 bucket_error[worst] = error;
503 for(int b = order.GetCount() - 2; b >= 0 && bucket_error[order[b]] > bucket_error[order[b + 1]]; b--)
504 Swap(order[b], order[b + 1]);
505 }
506
507 int bucket = 0;
508 // double y0, ym, y1;
509 for(int b = 0; b < buckets; b++)
510 {
511 int nsec = bucket_sections[b];
512 LLOG("# bucket sections[" << b << "] = " << bucket_sections[b]);
513 double buck_begin = xmin + b * step;
514 double sec_step = step / nsec;
515 for(int s = 0; s < nsec; s++)
516 {
517 double sec_begin = buck_begin + s * sec_step;
518 double y0 = fn(sec_begin);
519 double ym = fn(sec_begin + sec_step / 2);
520 double y1 = fn(sec_begin + sec_step);
521 // double sample_step = sec_step / samples;
522 // for(int m = 0; m <= samples; m++)
523 // {
524 // xsmpl[m] = m / double(samples);
525 // fn(sec_begin + sample_step * m, ysmpl[m]);
526 // }
527 double a0 = y0;
528 double a1 = 4 * ym - y1 - 3 * y0;
529 double a2 = y1 - y0 - a1;
530 abc[3 * (bucket + s) + 0] = a0 - a1 * s + a2 * s * s;
531 abc[3 * (bucket + s) + 1] = a1 * nsec - 2 * a2 * s * nsec;
532 abc[3 * (bucket + s) + 2] = a2 * nsec * nsec;
533 // QuadraticLeastSquare(xsmpl, ysmpl, samples + 1, abc.GetIter(3 * (bucket_index[b] + s)));
534 }
535
536 bucket_index[b] = bucket | (nsec << 16);
537 bucket += nsec;
538 }
539
540 /*
541 for(int i = 0; i < split; i++)
542 {
543 double left[3][3], right[3];
544 ZeroArray(left);
545 ZeroArray(right);
546 double xpow[3] = { 1 };
547 for(int s = 0; s <= samples; s++)
548 {
549 xpow[1] = s / (double)samples;
550 xpow[2] = sqr(xpow[1]);
551 double sy;
552 fn(xmin + (i + xpow[1]) * step, sy);
553 for(int y = 0; y < 3; y++)
554 {
555 double xy = xpow[y];
556 const double *xp = xpow;
557 double *dest = left[y];
558 for(int x = 0; x < 3; x++)
559 *dest++ += *xp++ * xy;
560 right[y] += sy * xy;
561 }
562 double D = Determinant(
563 left[0][0], left[0][1], left[0][2],
564 left[1][0], left[1][1], left[1][2],
565 left[2][0], left[2][1], left[2][2]);
566 double D1 = Determinant(
567 right[0], left[0][1], left[0][2],
568 right[1], left[1][1], left[1][2],
569 right[2], left[2][1], left[2][2]);
570 double DX = Determinant(
571 left[0][0], right[0], left[0][2],
572 left[1][0], right[1], left[1][2],
573 left[2][0], right[2], left[2][2]);
574 double DXX = Determinant(
575 left[0][0], left[0][1], right[0],
576 left[1][0], left[1][1], right[1],
577 left[2][0], left[2][1], right[2]);
578 abc[3 * i + 0] = D1 / D;
579 abc[3 * i + 1] = DX / D;
580 abc[3 * i + 2] = DXX / D;
581 }
582 }
583 */
584 /*
585 Vector<double> ys;
586 split &= ~1;
587 ys.SetCount(split + 1);
588 double step = (xmax_ - xmin_) / split;
589 for(int i = 0; i <= split; i++)
590 fn(xmin_ + i * step, ys[i]);
591 Create(xmin_, xmax_, ys, epsilon);
592 */
593 }
594
CreateInverse(double xmin,double xmax,const GisFunction & fn,int buckets,int sections,int samples)595 void GisInterpolator::CreateInverse(double xmin, double xmax, const GisFunction& fn, int buckets, int sections, int samples)
596 {
597 GisInverse inverse(xmin, xmax, fn, buckets, 1e-10);
598 Create(inverse.GetYMin(), inverse.GetYMax(), inverse, buckets, sections, samples);
599 }
600
CreateDump(double xmin_,double xmax_,const GisFunction & fn,int buckets,int sections,int samples,int check)601 String GisInterpolator::CreateDump(double xmin_, double xmax_, const GisFunction& fn, int buckets, int sections, int samples, int check)
602 {
603 Create(xmin_, xmax_, fn, buckets, sections, samples);
604 return Dump(fn, check);
605 }
606
CreateInverseDump(double xmin_,double xmax_,const GisFunction & fn,int buckets,int sections,int samples,int check)607 String GisInterpolator::CreateInverseDump(double xmin_, double xmax_, const GisFunction& fn, int buckets, int sections, int samples, int check)
608 {
609 CreateInverse(xmin_, xmax_, fn, buckets, sections, samples);
610 return Dump(fn, check);
611 }
612
Dump(const GisFunction & fn,int check)613 String GisInterpolator::Dump(const GisFunction& fn, int check)
614 {
615 /*
616 String out = NFormat("Interpolator: index(%d), abc(%d), xmin = %4v, xmax = %4v\n"
617 " X Y I D\n",
618 index.GetCount(), abc.GetCount(), xmin, xmax);
619 //*/
620 // String out = NFormat("Interpolator: y(%d), xmin %4v, xmax %4v\n"
621 // " X Y I D\n", y.GetCount(), xmin, xmax);
622 String out = NFormat("Interpolator: abc(%d), xmin %4v, xmax %4v\n"
623 " X Y I D\n", abc.GetCount(), xmin, xmax);
624 for(int t = 0; t < check; t++)
625 {
626 double x = xmin + t * (xmax - xmin) / (check - 1);
627 double f = fn(x);
628 double i = Get(x);
629 out << NFormat("%15>10!nf %15>10!nf %15>10!nf %15>10!nf\n", x, f, i, f - i);
630 }
631 return out;
632 }
633
Get(double x) const634 double GisInterpolator::Get(double x) const
635 {
636 int i = (int)(x = (x - xmin) / step);
637 if(i < 0) i = 0;
638 else if(i > limit) i = limit;
639 x -= i;
640 unsigned buck_sec = bucket_index[i];
641 int nsec = (buck_sec >> 16);
642 int bucket = (word)buck_sec;
643 // int b = bucket_index[i], nsec = bucket_index[i + 1] - b;
644 // int b = 0, nsec = 1;
645 if(nsec > 1)
646 {
647 int s = (int)(x * nsec);
648 if(s < 0)
649 s = 0;
650 else if(s >= nsec)
651 s = nsec - 1;
652 // x -= s;
653 bucket += s;
654 }
655 const double *a = &abc[3 * bucket];
656 return a[0] + x * (a[1] + x * a[2]);
657
658 /*
659 x = (x - xmin) / step;
660 int i = x < 0 ? 0 : x >= limit ? limit : (int)x;
661 const double *a = abc.GetIter(index[i]);
662 x = (x - begin[i]) / len[i];
663 return a[0] + x * (a[1] + x * a[2]);
664 //*/
665
666 /*
667 x = (x - xmin) * divisor;
668 int ifac = (x < 0 ? 0 : x >= limit ? limit : (int)x);
669 x -= ifac;
670 const double *yfac = &y[2 * ifac];
671 return yfac[0] + x * (yfac[1] + x * (yfac[2] - yfac[1] - yfac[0]));
672 //*/
673 }
674
GisInterpolatorDelta(double xmin,double xmax,const GisFunction & fn,int buckets,int sections,int samples,int check)675 String GisInterpolatorDelta(double xmin, double xmax, const GisFunction& fn, int buckets, int sections, int samples, int check)
676 {
677 GisInterpolator interpolator;
678 String dump = interpolator.CreateDump(xmin, xmax, fn, buckets, sections, samples, min(check, 1000));
679 double dmax = 0;
680 double step_check = (xmax - xmin) / check;
681 for(int ix = 0; ix <= check; ix++)
682 {
683 double fx = xmin + ix * step_check;
684 double fy = fn(fx);
685 double fiy = interpolator(fx);
686 double d = fabs(fiy - fy);
687 if(d > dmax)
688 dmax = d;
689 }
690 return NFormat("d-max = %4v\n\n%s", dmax, dump);
691 }
692
GisInterpolatorTiming(double xmin,double xmax,const GisFunction & fn,int buckets,int sections,int samples,int check)693 String GisInterpolatorTiming(double xmin, double xmax, const GisFunction& fn, int buckets, int sections, int samples, int check)
694 {
695 GisInterpolator interpolator;
696 String dump = interpolator.CreateDump(xmin, xmax, fn, buckets, sections, samples, check);
697 double dmax = 0;
698 double step_check = (xmax - xmin) / check;
699 // double step_64K = (xmax - xmin) / 65536;
700 Buffer<double> check_table(1000);
701 for(int c = 0; c < 1000; c++)
702 check_table[c] = xmin + c / 999.0;
703 for(int ix = 0; ix <= check; ix++)
704 {
705 double fx = xmin + ix * step_check;
706 double fy = fn(fx);
707 double fiy = interpolator(fx);
708 dmax = max(dmax, fabs(fiy - fy));
709 }
710 int start, duration_e, duration_f, duration_i;
711 int count_e = 0, count_f = 0, count_i = 0;
712 start = msecs();
713 while((duration_e = msecs(start)) == 0);
714 start += duration_e;
715 Callback2<double, double&> empty;
716 double y;
717 do
718 {
719 count_e++;
720 for(int t = 0; t < 1000; t++)
721 empty(check_table[t], y);
722 }
723 while((duration_e = msecs(start)) < 500);
724 start += duration_e;
725 do
726 {
727 count_f++;
728 for(int t = 0; t < 1000; t++)
729 y = fn(check_table[t]);
730 }
731 while((duration_f = msecs(start)) < 500);
732 start += duration_f;
733 do
734 {
735 count_i++;
736 for(int t = 0; t < 1000; t++)
737 interpolator(check_table[t]);
738 }
739 while((duration_i = msecs(start)) < 500);
740 double e_nsecs = duration_e * 1000 / count_e;
741 double f_nsecs = duration_f * 1000 / count_f - e_nsecs;
742 double i_nsecs = duration_i * 1000 / count_i - e_nsecs;
743
744 return NFormat("d-max = %4v, f = %4v nsecs, i = %4v nsecs\n\n%s", dmax, f_nsecs, i_nsecs, dump);
745 }
746
747 /*
748 void Interpolator::Create(Callback2<double, double&> calc, double xmin, double xmax, int min_depth, int max_depth, double epsilon)
749 {
750 calc(extent.left = xmin, extent.top);
751 calc(extent.right = xmax, extent.bottom);
752 Pointf mid;
753 calc(mid.x = (extent.left + extent.right) * 0.5, mid.y);
754 scale = 1 << max_depth;
755 divisor = extent.Width() / scale;
756 tree = CreateTree(calc, extent, mid, 1, min_depth, max_depth, epsilon);
757 }
758 */
759
760 /*
761 One<Interpolator::Tree> Interpolator::CreateTree(Callback2<double, double&> calc, const Rectf& extent, const Pointf& mid,
762 int depth, int min_depth, int max_depth, double epsilon)
763 {
764 One<Tree> out = new Tree;
765 out->mid_y = mid.y;
766 if(++depth <= max_depth)
767 {
768 Pointf lmid, rmid;
769 calc(lmid.x = (extent.left + mid.x) / 2, lmid.y);
770 calc(rmid.x = (mid.x + extent.right) / 2, rmid.y);
771 double a1_2 = (extent.bottom - extent.top) / 4;
772 double a2_4_a0 = (extent.top + extent.bottom + mid.y * 6) / 8;
773 if(depth <= min_depth || fabs(a2_4_a0 - a1_2 - lmid.y) > epsilon)
774 out->left = CreateTree(calc, Rectf(extent.left, extent.top, mid.x, mid.y), lmid, depth, min_depth, max_depth, epsilon);
775 if(depth <= min_depth || fabs(a2_4_a0 + a1_2 - rmid.y) > epsilon)
776 out->right = CreateTree(calc, Rectf(mid.x, mid.y, extent.right, extent.bottom), rmid, depth, min_depth, max_depth, epsilon);
777 }
778 return out;
779 }
780 */
781
782 /*
783 double Interpolator::operator [] (double x) const
784 {
785 x = (x - extent.left) / divisor;
786 int ifac = ffloor(minmax<double>(x, 0, scale - 1));
787 x -= ifac;
788 int bit = scale >> 1;
789 const Tree *node = ~tree;
790 double ymin = extent.top, ymax = extent.bottom;
791 for(;;)
792 if(ifac & bit)
793 {
794 if(node->right)
795 {
796 ymin = node->mid_y;
797 node = ~node->right;
798 bit >>= 1;
799 }
800 else
801 break;
802 }
803 else
804 {
805 if(node->left)
806 {
807 ymax = node->mid_y;
808 node = ~node->left;
809 bit >>= 1;
810 }
811 else
812 break;
813 }
814 x = (x + ((ifac & (2 * bit - 1)) - bit)) / bit;
815 double a1 = (ymax - ymin) / 2;
816 double a2 = (ymin + ymax) / 2 - node->mid_y;
817 return node->mid_y + x * (a1 + x * a2);
818 }
819 */
820
821 /*
822 double Interpolator::Linear(double x) const
823 {
824 x = (x - extent.left) / divisor;
825 int ifac = ffloor(minmax<double>(x, 0, scale - 1));
826 x -= ifac;
827 int bit = scale;
828 double ymin = extent.top, ymax = extent.bottom;
829 for(const Tree *node = ~tree; node;)
830 if(ifac & (bit >>= 1))
831 {
832 ymin = node->mid_y;
833 node = ~node->right;
834 }
835 else
836 {
837 ymax = node->mid_y;
838 node = ~node->left;
839 }
840 x = (x + ((ifac & (bit - 1)))) / bit;
841 return ymin + (ymax - ymin) * x;
842 }
843 */
844
GisOrientation(Pointf p)845 GisOrientation::GisOrientation(Pointf p)
846 {
847 static const double EPS = 1 / 3600.0;
848 pole = p;
849 delta_phi = M_PI / 2 - pole.y * DEGRAD;
850 // int lquad = ffloor((pole.y + EPS) / (M_PI / 2));
851 // double reduced = pole.y - lquad * (M_PI / 2);
852 pole.x = modulo(pole.x + 180, 360) - 180;
853 // pole.y -= lquad * (M_PI / 2);
854 // int gquad = lquad;
855 identity = false;
856 if(fabs(delta_phi) <= EPS)
857 {
858 if(fabs(pole.x) <= EPS)
859 {
860 identity = true;
861 localproc = globalproc = &GisOrientation::Identity;
862 localextent = globalextent = &GisOrientation::IdentityExtent;
863 }
864 else
865 {
866 localproc = &GisOrientation::LocalDelta;
867 globalproc = &GisOrientation::GlobalDelta;
868 localextent = &GisOrientation::LocalDeltaExtent;
869 globalextent = &GisOrientation::GlobalDeltaExtent;
870 }
871 }
872 else
873 {
874 // lquad++;
875 suk = sin(pole.y * DEGRAD);
876 sukneg = (suk < 0);
877 cuk = cos(pole.y * DEGRAD);
878 cukneg = (cuk < 0);
879 localproc = &GisOrientation::LocalAny;
880 globalproc = &GisOrientation::GlobalAny;
881 localextent = &GisOrientation::LocalAnyExtent;
882 globalextent = &GisOrientation::GlobalAnyExtent;
883 }
884 /*
885 switch(lquad & 3)
886 {
887 case 1: localproc = localquad; break;
888 case 2: localproc = &WorldTransform::Local1; break;
889 case 3: localproc = &WorldTransform::Local2; break;
890 case 0: localproc = &WorldTransform::Local3; break;
891 }
892 switch(gquad & 3)
893 {
894 case 1: globalproc = globalquad; break;
895 case 0: globalproc = &WorldTransform::Global1; break;
896 case 3: globalproc = &WorldTransform::Global2; break;
897 case 2: globalproc = &WorldTransform::Global3; break;
898 }
899 */
900 }
901
LocalAny(double lon,double lat) const902 Pointf GisOrientation::LocalAny(double lon, double lat) const
903 {
904 double dv = (lon - pole.x) * DEGRAD;
905 lat *= DEGRAD;
906 double su = sin(lat), cu = cos(lat), sv = sin(dv), cv = cos(dv);
907 double cuv = cu * cv;
908 double s = suk * su + cuk * cuv;
909 double d = dv;
910 if(s <= -1.0)
911 s = -M_PI / 2;
912 else if(s >= +1.0)
913 s = +M_PI / 2;
914 else
915 {
916 double cs = sqrt(1 - s * s);
917 s = asin(s);
918 d = SafeArcCos((suk * cuv - cuk * su) / cs, sv < 0);
919 }
920 return Pointf(d / DEGRAD, s / DEGRAD);
921 }
922
GlobalAny(double lon,double lat) const923 Pointf GisOrientation::GlobalAny(double lon, double lat) const
924 {
925 lon *= DEGRAD;
926 lat *= DEGRAD;
927 double su = sin(lat), cu = cos(lat), sv = sin(lon), cv = cos(lon);
928 double cuv = cu * cv;
929 double s = suk * su - cuk * cuv;
930 double d = lon;
931 if(s <= -1.0)
932 s = -M_PI / 2;
933 else if(s >= +1.0)
934 s = +M_PI / 2;
935 else
936 {
937 double cs = sqrt(1 - s * s);
938 s = asin(s);
939 d = SafeArcCos((suk * cuv + cuk * su) / cs, sv < 0);
940 }
941 return Pointf(d / DEGRAD + pole.x, s / DEGRAD);
942 }
943
LocalDelta(double lon,double lat) const944 Pointf GisOrientation::LocalDelta(double lon, double lat) const
945 {
946 return Pointf(lon - pole.x, lat);
947 }
948
GlobalDelta(double lon,double lat) const949 Pointf GisOrientation::GlobalDelta(double lon, double lat) const
950 {
951 return Pointf(lon + pole.x, lat);
952 }
953
Identity(double lon,double lat) const954 Pointf GisOrientation::Identity(double lon, double lat) const
955 {
956 return Pointf(lon, lat);
957 }
958
CalcRatio(double x,double y)959 static inline double CalcRatio(double x, double y)
960 {
961 double den = 1 - y * y;
962 return den > 0 ? x / sqrt(den) : double(Null);
963 }
964
LocalAnyExtent(const Rectf & lonlat) const965 Rectf GisOrientation::LocalAnyExtent(const Rectf& lonlat) const
966 {
967 Rectf out = Null;
968 out |= Local(lonlat.TopLeft());
969 out |= Local(lonlat.TopRight());
970 out |= Local(lonlat.BottomLeft());
971 out |= Local(lonlat.BottomRight());
972 return out;
973 /*
974 // if(lonlat.Width() >= 2 * M_PI)
975 // return Rectf(-M_PI, -M_PI / 2, +M_PI, +M_PI / 2);
976 double dv1 = (lonlat.left - pole.x) * DEGRAD, sv1 = sin(dv1), cv1 = cos(dv1);
977 double dv2 = (lonlat.right - pole.x) * DEGRAD, sv2 = sin(dv2), cv2 = cos(dv2);
978 double trad = lonlat.top * DEGRAD, brad = lonlat.bottom * DEGRAD;
979 int xmask = DegreeMask(dv1, dv2);
980 int yfmask = DegreeMask(trad + delta_phi, brad + delta_phi);
981 int ybmask = DegreeMask(trad - delta_phi, brad - delta_phi);
982 if(xmask & AM_E0 && yfmask & AM_E1 && yfmask & AM_E3
983 || xmask & AM_E2 && ybmask & AM_E1 && ybmask & AM_E3)
984 return Rectf(-M_PI, -M_PI / 2, +M_PI, +M_PI / 2);
985 double su1 = sin(trad), su2 = sin(brad);
986 double cu1 = cos(trad), cu2 = cos(brad);
987 double ccv1 = cuk * cv1, ccv2 = cuk * cv2;
988 double cvmin = ccv1, cvmax = ccv2;
989 if(cvmin > cvmax)
990 Swap(cvmin, cvmax);
991 if(xmask & AM_E2) (cukneg ? cvmax : cvmin) = -cuk;
992 if(xmask & AM_E0) (cukneg ? cvmin : cvmax) = +cuk;
993 double suks1 = suk * su1, suks2 = suk * su2;
994 double smin = min(suks1 + cu1 * cvmin, suks2 + cu2 * cvmin);
995 double smax = max(suks1 + cu1 * cvmax, suks2 + cu2 * cvmax);
996 if(xmask & (AM_E0 | AM_E2))
997 {
998 if(xmask & AM_E0 ? yfmask & AM_E1 : ybmask & AM_E3)
999 return Rectf(-M_PI, SafeArcSin(smin), +M_PI, +M_PI / 2);
1000 if(xmask & AM_E0 ? yfmask & AM_E3 : ybmask & AM_E1)
1001 return Rectf(-M_PI, -M_PI / 2, +M_PI, SafeArcSin(smin));
1002 }
1003 // if(ymask & AM_E1)
1004 // smax = (sukneg ? min(su1, su2) : max(su1, su2)) / suk;
1005 // if(ymask & AM_E3)
1006 // smin = (sukneg ? max(su1, su2) : min(su1, su2)) / suk;
1007 double cuks1 = cuk * su1, cuks2 = cuk * su2;
1008 double scv1 = cv1 * suk, scv2 = cv2 * suk;
1009 double lt = CalcRatio(scv1 * cu1 - cuks1, ccv1 * cu1 + suks1);
1010 double lb = CalcRatio(scv1 * cu2 - cuks2, ccv1 * cu2 + suks2);
1011 double rt = CalcRatio(scv2 * cu1 - cuks1, ccv2 * cu1 + suks1);
1012 double rb = CalcRatio(scv2 * cu2 - cuks2, ccv2 * cu2 + suks2);
1013 if(IsNull(lt)) lt = lb; else if(IsNull(lb)) lb = lt;
1014 if(IsNull(rt)) rt = rb; else if(IsNull(rb)) rb = rt;
1015 double cmin = -M_PI, cmax = +M_PI;
1016 if(yfmask & AM_E1)
1017 {
1018 cmin = (lt >= rt ? SafeArcCos(lt, sv1 < 0) : SafeArcCos(rt, sv2 < 0));
1019 cmax = (lb <= rb ? SafeArcCos(lb, sv1 < 0) : SafeArcCos(rb, sv2 < 0));
1020 if(sv1 < 0)
1021 Swap(cmin, cmax);
1022 }
1023 else if(yfmask & AM_E3)
1024 {
1025 cmin = (lb >= rb ? SafeArcCos(lb, sv1 < 0) : SafeArcCos(rb, sv2 < 0));
1026 cmax = (lt <= rt ? SafeArcCos(lt, sv1 < 0) : SafeArcCos(rt, sv2 < 0));
1027 if(sv1 < 0)
1028 Swap(cmin, cmax);
1029 }
1030 else if(yfmask & (AM_Q0 | AM_Q3))
1031 { // front octants
1032 cmin = SafeArcCos(sv1 >= 0 ? max(lt, lb) : min(lt, lb), sv1 < 0);
1033 cmax = SafeArcCos(sv2 >= 0 ? min(rt, rb) : max(rt, rb), sv2 < 0);
1034 }
1035 else
1036 {
1037 cmin = SafeArcCos(sv1 >= 0 ? max(rt, rb) : min(rt, rb), sv2 < 0);
1038 cmax = SafeArcCos(sv2 >= 0 ? min(lt, lb) : max(lt, lb), sv1 < 0);
1039 }
1040 // return Rectf(cmin, SafeArcSin(smin), cmax >= cmin ? cmax : cmax + 2 * M_PI, SafeArcSin(smax));
1041 cmin /= DEGRAD;
1042 cmax /= DEGRAD;
1043 return Rectf(cmin, -90, cmax >= cmin ? cmax : cmax + 360, 90);
1044 */
1045 }
1046
GlobalAnyExtent(const Rectf & lonlat) const1047 Rectf GisOrientation::GlobalAnyExtent(const Rectf& lonlat) const
1048 {
1049 Rectf out = Null;
1050 out |= Global(lonlat.TopLeft());
1051 out |= Global(lonlat.TopRight());
1052 out |= Global(lonlat.BottomLeft());
1053 out |= Global(lonlat.BottomRight());
1054 return out;
1055 }
1056
LocalDeltaExtent(const Rectf & lonlat) const1057 Rectf GisOrientation::LocalDeltaExtent(const Rectf& lonlat) const
1058 {
1059 return lonlat.OffsetedHorz(-pole.x);
1060 }
1061
GlobalDeltaExtent(const Rectf & lonlat) const1062 Rectf GisOrientation::GlobalDeltaExtent(const Rectf& lonlat) const
1063 {
1064 return lonlat.OffsetedHorz(pole.x);
1065 }
1066
IdentityExtent(const Rectf & lonlat) const1067 Rectf GisOrientation::IdentityExtent(const Rectf& lonlat) const
1068 {
1069 return lonlat;
1070 }
1071
1072 /*
1073 Pointf WorldTransform::Local1(const Pointf& pt) const
1074 {
1075 return Pointf((pt.x >= 0
1076 Pointf out = (this->*localquad)(pt);
1077 return Pointf((out.x >= 0 ? M_PI / 2 + out.y : -M_PI / 2 - out.y), M_PI / 2 - fabs(out.x));
1078 }
1079
1080 Pointf WorldTransform::Local2(const Pointf& pt) const
1081 {
1082 Pointf out = (this->*localquad)(pt);
1083 return Pointf(M_PI - out.x, -out.y);
1084 }
1085
1086 Pointf WorldTransform::Local3(const Pointf& pt) const
1087 {
1088 Pointf out = (this->*localquad)(pt);
1089 return Pointf((out.x >= 0 ? M_PI / 2 - out.y : -M_PI / 2 + out.y), fabs(out.x) - M_PI / 2);
1090 }
1091 */
1092
1093 /*
1094 Pointf WorldTransform::Global1(const Pointf& pt) const
1095 {
1096 Pointf out = (this->*globalquad)(pt);
1097 return Pointf((out.x >= 0 ? M_PI / 2 - out.y : -M_PI / 2 + out.y), fabs(out.x) - M_PI / 2);
1098 }
1099
1100 Pointf WorldTransform::Global2(const Pointf& pt) const
1101 {
1102 Pointf out = (this->*globalquad)(pt);
1103 return Pointf(M_PI - out.x, -out.y);
1104 }
1105
1106 Pointf WorldTransform::Global3(const Pointf& pt) const
1107 {
1108 Pointf out = (this->*globalquad)(pt);
1109 return Pointf((out.x >= 0 ? M_PI / 2 + out.y : -M_PI / 2 - out.y), M_PI / 2 - fabs(out.x));
1110 }
1111 */
1112
Calculate(const GisTransform & transform,const Rectf & src)1113 void Gis2DPolynome::Calculate(const GisTransform& transform, const Rectf& src)
1114 {
1115 int xinter = 10, yinter = 10;
1116 LinearSolver xsolv(COEF_COUNT), ysolv(COEF_COUNT);
1117 double bases[COEF_COUNT];
1118 for(int ix = 0; ix <= xinter; ix++)
1119 for(int iy = 0; iy <= yinter; iy++) {
1120 double x = ix / (double)xinter, y = iy / (double)yinter;
1121 double x2 = x * x, y2 = y * y;
1122 Pointf dest = transform.Target(src.TopLeft() + src.Size() * Sizef(x, y));
1123 bases[COEF_1] = 1;
1124 bases[COEF_X] = x;
1125 bases[COEF_Y] = y;
1126 bases[COEF_X2] = x2;
1127 bases[COEF_XY] = x * y;
1128 bases[COEF_Y2] = y2;
1129 bases[COEF_X3] = x2 * x;
1130 bases[COEF_X2Y] = x2 * y;
1131 bases[COEF_XY2] = x * y2;
1132 bases[COEF_Y3] = y2 * y;
1133 xsolv.AddLSI(bases, dest.x);
1134 ysolv.AddLSI(bases, dest.y);
1135 }
1136 Vector<double> xcoef = xsolv.Solve();
1137 Vector<double> ycoef = ysolv.Solve();
1138 for(int i = 0; i < COEF_COUNT; i++)
1139 coef[i] = Sizef(xcoef[i], ycoef[i]);
1140 }
1141
Transform(double x,double y) const1142 Pointf Gis2DPolynome::Transform(double x, double y) const
1143 {
1144 double x2 = x * x, y2 = y * y;
1145 return coef[COEF_1]
1146 + coef[COEF_X] * x
1147 + coef[COEF_Y] * y
1148 + coef[COEF_XY] * (x * y)
1149 + coef[COEF_X2] * x2
1150 + coef[COEF_Y2] * y2
1151 + coef[COEF_X3] * (x2 * x)
1152 + coef[COEF_X2Y] * (x2 * y)
1153 + coef[COEF_XY2] * (x * y2)
1154 + coef[COEF_Y3] * (y2 * y)
1155 ;
1156 }
1157
Gis2DGrid(const Sizef & block_size_,const Rect & block_limit_)1158 Gis2DGrid::Gis2DGrid(const Sizef& block_size_, const Rect& block_limit_)
1159 : block_size(block_size_)
1160 , block_limit(block_limit_)
1161 {
1162 block_span = Rectf(0, 0, 0, 0);
1163 }
1164
GetBlockIndex(const Pointf & point) const1165 Point Gis2DGrid::GetBlockIndex(const Pointf& point) const
1166 {
1167 return Point(ffloor(point.x / block_size.cx), ffloor(point.y / block_size.cy));
1168 }
1169
GetBlockSpan(const Rectf & rc) const1170 Rect Gis2DGrid::GetBlockSpan(const Rectf& rc) const
1171 {
1172 return Rect(ffloor(rc.left / block_size.cx), ffloor(rc.top / block_size.cy),
1173 ffloor(rc.right / block_size.cx) + 1, ffloor(rc.bottom / block_size.cy) + 1);
1174 }
1175
Transform(const Pointf & pt) const1176 Pointf Gis2DGrid::Transform(const Pointf& pt) const
1177 {
1178 Point block = GetBlockIndex(pt);
1179 if(const Gis2DPolynome *poly = GetBlock(block))
1180 return poly->Transform(pt.x / block_size.cx - block.x, pt.y / block_size.cy - block.y);
1181 return Null;
1182 }
1183
GetBlock(int x,int y) const1184 const Gis2DPolynome *Gis2DGrid::GetBlock(int x, int y) const
1185 {
1186 return (x >= block_span.left && x < block_span.right && y >= block_span.top && y < block_span.bottom
1187 ? &block_rows[y - block_span.top][x - block_span.left]
1188 : NULL);
1189 }
1190
SizeOf() const1191 int Gis2DGrid::SizeOf() const
1192 {
1193 return block_span.Width() * block_span.Height() * (sizeof(Gis2DPolynome) + 32) + sizeof(*this);
1194 }
1195
Grow(const GisTransform & transform,const Rectf & extent)1196 void Gis2DGrid::Grow(const GisTransform& transform, const Rectf& extent)
1197 {
1198 Rect target_span = GetBlockSpan(extent) & block_limit;
1199 if(block_span.Contains(target_span))
1200 return;
1201 if(block_span.IsEmpty())
1202 block_span = Rect(target_span.left, target_span.top, target_span.left, target_span.top);
1203 target_span |= block_span;
1204 int add_left = block_span.left - target_span.left;
1205 int add_right = target_span.right - block_span.right;
1206 ASSERT(add_left >= 0 && add_right >= 0);
1207 Rectf blk_extent(block_size.cx * block_span.left, block_size.cy * block_span.top,
1208 block_size.cx * block_span.right, block_size.cy * block_span.bottom);
1209 if(add_left || add_right) {
1210 Rectf row_extent = blk_extent;
1211 row_extent.bottom = row_extent.top + block_size.cy;
1212 for(int i = 0; i < block_rows.GetCount(); i++) {
1213 BiArray<Gis2DPolynome>& row = block_rows[i];
1214 if(add_left) {
1215 Rectf cell = row_extent;
1216 for(int n = 0; n < add_left; n++) {
1217 cell.right = cell.left;
1218 cell.left -= block_size.cx;
1219 row.AddHead().Calculate(transform, cell);
1220 }
1221 }
1222 if(add_right) {
1223 Rectf cell = row_extent;
1224 for(int n = 0; n < add_right; n++) {
1225 cell.left = cell.right;
1226 cell.right += block_size.cx;
1227 row.AddTail().Calculate(transform, cell);
1228 }
1229 }
1230 row_extent.OffsetVert(block_size.cy);
1231 }
1232 block_span.Inflate(add_left, 0, add_right, 0);
1233 blk_extent.Inflate(-block_size.cx * add_left, 0, block_size.cx * add_right, 0);
1234 }
1235 int add_top = block_span.top - target_span.top;
1236 if(add_top) {
1237 Rectf cell = blk_extent;
1238 for(int i = 0; i < add_top; i++) {
1239 cell.bottom = cell.top;
1240 cell.top -= block_size.cy;
1241 BiArray<Gis2DPolynome>& top = block_rows.AddHead();
1242 cell.right = blk_extent.left;
1243 for(int j = block_span.left; j < block_span.right; j++) {
1244 cell.left = cell.right;
1245 cell.right += block_size.cx;
1246 top.AddTail().Calculate(transform, cell);
1247 }
1248 }
1249 block_span.top -= add_top;
1250 blk_extent.top -= add_top * block_size.cy;
1251 }
1252 int add_bottom = target_span.bottom - block_span.bottom;
1253 if(add_bottom) {
1254 Rectf cell = blk_extent;
1255 for(int i = 0; i < add_bottom; i++) {
1256 cell.top = cell.bottom;
1257 cell.bottom += block_size.cy;
1258 BiArray<Gis2DPolynome>& bottom = block_rows.AddTail();
1259 cell.right = blk_extent.left;
1260 for(int j = block_span.left; j < block_span.right; j++) {
1261 cell.left = cell.right;
1262 cell.right += block_size.cx;
1263 bottom.AddTail().Calculate(transform, cell);
1264 }
1265 }
1266 block_span.bottom += add_bottom;
1267 blk_extent.bottom += add_bottom * block_size.cy;
1268 }
1269 }
1270
CreateLinearSplit(Point s1,Point s2,Point t1,Point t2,const SegmentTreeInfo & info,int depth)1271 static One<LinearSegmentTree::Node> CreateLinearSplit(Point s1, Point s2, Point t1, Point t2, const SegmentTreeInfo& info, int depth)
1272 {
1273 double m = info.src_trg.SourceDeviation(Pointf(s1) * info.img_src, Pointf(s2) * info.img_src);
1274 if(m <= info.max_deviation)
1275 return NULL;
1276 One<LinearSegmentTree::Node> split = new LinearSegmentTree::Node;
1277 split->source = (s1 + s2) >> 1;
1278 split->target = split->source * info;
1279 if(++depth <= info.max_depth)
1280 {
1281 split->below = CreateLinearSplit(s1, split->source, t1, split->target, info, depth);
1282 split->above = CreateLinearSplit(split->source, s2, split->target, t2, info, depth);
1283 }
1284 return split;
1285 }
1286
CreateLinearTree(Point s1,Point s2,const SegmentTreeInfo & info)1287 LinearSegmentTree CreateLinearTree(Point s1, Point s2, const SegmentTreeInfo& info)
1288 {
1289 LinearSegmentTree out;
1290 out.source1 = s1;
1291 out.source2 = s2;
1292 out.target1 = s1 * info;
1293 out.target2 = s2 * info;
1294 out.split = CreateLinearSplit(out.source1, out.source2, out.target1, out.target2, info, 0);
1295 return out;
1296 }
1297
CreatePlanarSplit(PlanarSegmentTree::Node & node,const LinearSegmentTree::Node * left,const LinearSegmentTree::Node * top,const LinearSegmentTree::Node * right,const LinearSegmentTree::Node * bottom,const Rect & src,const SegmentTreeInfo & info,int depth,Point trg_topleft,Point trg_topright,Point trg_bottomleft,Point trg_bottomright)1298 static void CreatePlanarSplit(PlanarSegmentTree::Node& node,
1299 const LinearSegmentTree::Node *left, const LinearSegmentTree::Node *top,
1300 const LinearSegmentTree::Node *right, const LinearSegmentTree::Node *bottom,
1301 const Rect& src, const SegmentTreeInfo& info, int depth,
1302 Point trg_topleft, Point trg_topright, Point trg_bottomleft, Point trg_bottomright)
1303 {
1304 double m = info.src_trg.SourceExtentDeviation(Rectf(src) * info.img_src);
1305 node.source = src;
1306 node.trg_topleft = trg_topleft;
1307 node.trg_topright = trg_topright;
1308 node.trg_bottomleft = trg_bottomleft;
1309 node.trg_bottomright = trg_bottomright;
1310 if(m > info.max_deviation && ++depth <= info.max_depth) {
1311 node.split = new PlanarSegmentTree::Split;
1312 Point mid = src.CenterPoint();
1313 Point mtrg = mid * info;
1314 Point lmid = (left ? left->target : (trg_topleft + trg_bottomleft) >> 1); //src.CenterLeft() * info);
1315 Point tmid = (top ? top->target : (trg_topleft + trg_topright) >> 1); // src.TopCenter() * info);
1316 Point rmid = (right ? right->target : (trg_topright + trg_bottomright) >> 1); //src.CenterRight() * info);
1317 Point bmid = (bottom ? bottom->target : (trg_bottomleft + trg_bottomright) >> 1); // src.BottomCenter() * info);
1318 CreatePlanarSplit(node.split->topleft, left ? ~left->below : NULL, top ? ~top->below : NULL, NULL, NULL,
1319 Rect(src.left, src.top, mid.x, mid.y), info, depth, trg_topleft, tmid, lmid, mtrg);
1320 CreatePlanarSplit(node.split->topright, NULL, top ? ~top->above : NULL, right ? ~right->below : NULL, NULL,
1321 Rect(mid.x, src.top, src.right, mid.y), info, depth, tmid, trg_topright, mtrg, rmid);
1322 CreatePlanarSplit(node.split->bottomleft, left ? ~left->above : NULL, NULL, NULL, bottom ? ~bottom->below : NULL,
1323 Rect(src.left, mid.y, mid.x, src.bottom), info, depth, lmid, mtrg, trg_bottomleft, bmid);
1324 CreatePlanarSplit(node.split->bottomright, NULL, NULL, right ? ~right->above : NULL, bottom ? ~bottom->above : NULL,
1325 Rect(mid.x, mid.y, src.right, src.bottom), info, depth, mtrg, rmid, bmid, trg_bottomright);
1326 }
1327 }
1328
CreatePlanarTree(const LinearSegmentTree & left,const LinearSegmentTree & top,const LinearSegmentTree & right,const LinearSegmentTree & bottom,const SegmentTreeInfo & info)1329 PlanarSegmentTree CreatePlanarTree(const LinearSegmentTree& left, const LinearSegmentTree& top,
1330 const LinearSegmentTree& right, const LinearSegmentTree& bottom, const SegmentTreeInfo& info)
1331 {
1332 PlanarSegmentTree out;
1333 CreatePlanarSplit(out.root, ~left.split, ~top.split, ~right.split, ~bottom.split,
1334 Rect(left.source1, right.source2), info, 0,
1335 left.target1, right.target1, left.target2, right.target2);
1336 return out;
1337 }
1338
GisCoordsGaussLatitude()1339 GisCoordsGaussLatitude::GisCoordsGaussLatitude()
1340 {
1341 }
1342
Get(double phi) const1343 double SphericalLatitudeFunction::Get(double phi) const
1344 {
1345 // RTIMING("SphericalLatitudeFunction::Get");
1346 phi *= DEGRAD;
1347 double esx = e * sin(phi);
1348 double eps = pow((1 - esx) / (1 + esx), e * alpha / 2) / k;
1349 double dpi = M_PI / 4 - phi / 2;
1350 if(dpi <= 0.001)
1351 {
1352 // RLOG("first dpi = " << FormatDouble(x, 5));
1353 // RLOG("saturation: " << dpi);
1354 return 90 - 2 / DEGRAD * (pow(fabs(dpi), alpha) / (dpi >= 0 ? eps : -eps));
1355 }
1356 else
1357 {
1358 double rho = phi / 2 + M_PI / 4;
1359 return 2 / DEGRAD * atan(pow(fabs(tan(rho)), alpha) * (rho >= 0 ? eps : -eps)) - 90;
1360 }
1361 }
1362
Create(double a,double e2,double base_parallel)1363 void GisCoordsGaussLatitude::Create(double a, double e2, double base_parallel)
1364 {
1365 double e = sqrt(e2);
1366 double phi0 = base_parallel * DEGRAD;
1367 double alpha = sqrt(1 + (e2 * sqr(sqr(cos(phi0)))) / (1 - e2));
1368 double sinphi = sin(phi0);
1369 double U0 = asin(sinphi / alpha);
1370 double k = exp(alpha * (log(tan(phi0 / 2 + M_PI / 4)) + e / 2 * log((1 - e * sinphi) / (1 + e * sinphi))))
1371 / tan(U0 / 2 + M_PI / 4);
1372 // k = pow(tan(base_parallel / 2 + M_PI / 4), alpha)
1373 // * pow((1 - e * sinphi) / (1 + e * sinphi), alpha * e / 2)
1374 // / tan(U0 / 2 + M_PI / 4);
1375 radius = a * sqrt(1 - e2) / (1 - e2 * sqr(sinphi));
1376 gauss_projected.Clear();
1377 gauss_latitude.Clear();
1378
1379 SphericalLatitudeFunction gslf(alpha, k, radius, e, U0);
1380 //gslf.Dump(-1.58, +1.58, 1000);
1381 //gslf.Dump(-1.58, -1.56, 1000);
1382 //gslf.Dump(+1.56, +1.58, 1000);
1383 gauss_projected.Create(base_parallel - 30, base_parallel + 30, gslf, 300, 5000, 4);
1384 gauss_latitude.CreateInverse(base_parallel - 30, base_parallel + 30, gslf, 300, 5000, 4);
1385 }
1386
1387 }
1388