1 /* -*-c++-*- */
2 /* osgEarth - Geospatial SDK for OpenSceneGraph
3 * Copyright 2019 Pelican Mapping
4 * http://osgearth.org
5 *
6 * osgEarth is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>
18 */
19 #include <osgEarthUtil/FlatteningLayer>
20 #include <osgEarth/Registry>
21 #include <osgEarth/HeightFieldUtils>
22 #include <osgEarthFeatures/FeatureCursor>
23
24 using namespace osgEarth;
25 using namespace osgEarth::Util;
26 using namespace osgEarth::Features;
27 using namespace osgEarth::Symbology;
28
29 #define LC "[FlatteningLayer] "
30
31 #define OE_TEST OE_DEBUG
32
33 namespace
34 {
35 // linear interpolation between a and b
mix(double a,double b,double t)36 double inline mix(double a, double b, double t)
37 {
38 return a + (b-a)*t;
39 }
40
41 // smoothstep (cos approx) interpolation between a and b
smoothstep(double a,double b,double t)42 double inline smoothstep(double a, double b, double t)
43 {
44 // smoothstep (approximates cosine):
45 t = t*t*(3.0-2.0*t);
46 return a + (b-a)*t;
47 }
48
smootherstep(double a,double b,double t)49 double inline smootherstep(double a, double b, double t)
50 {
51 t = t*t*t*(t*(t*6.0 - 15.0)+10.0);
52 return a + (b-a)*t;
53 }
54
55 // clamp "a" to [lo..hi].
clamp(double a,double lo,double hi)56 double inline clamp(double a, double lo, double hi)
57 {
58 return osg::maximum(osg::minimum(a, hi), lo);
59 }
60
61 typedef osg::Vec3d POINT;
62 typedef osg::Vec3d VECTOR;
63
64 // iterator that permits the use of looping indexes.
65 template<typename T>
66 struct CircleIterator {
CircleIterator__anon71f210340111::CircleIterator67 CircleIterator(const std::vector<T>& v) : _v(v) { }
size__anon71f210340111::CircleIterator68 int size() const { return _v.size(); }
operator []__anon71f210340111::CircleIterator69 const T& operator[](int i) const { return _v[i % _v.size()]; }
70 const std::vector<T>& _v;
71 };
72
73 // is P inside the CCW triangle ABC?
triangleContains(const POINT & A,const POINT & B,const POINT & C,const POINT & P)74 bool triangleContains(const POINT& A, const POINT& B, const POINT& C, const POINT& P)
75 {
76 VECTOR AB = B-A, BC = C-B, CA = A-C;
77 if (((P - A) ^ AB).z() > 0.0) return false;
78 if (((P - B) ^ BC).z() > 0.0) return false;
79 if (((P - C) ^ CA).z() > 0.0) return false;
80 return true;
81 }
82
83 // Find a point internal to the polygon.
84 // This will not always work with polygons that contain holes,
85 // so we need to come up with a different algorithm if this becomes a problem.
86 // Maybe try a random point generator and profile it.
getInternalPoint(const Symbology::Polygon * p)87 osg::Vec3d inline getInternalPoint(const Symbology::Polygon* p)
88 {
89 // Simple test: if the centroid is in the polygon, use it.
90 osg::Vec3d centroid = p->getBounds().center();
91 if (p->contains2D(centroid.x(), centroid.y()))
92 return centroid;
93
94 // Concave/holey polygon, so try the hard way.
95 // Ref: http://apodeline.free.fr/FAQ/CGAFAQ/CGAFAQ-3.html
96
97 CircleIterator<POINT> vi(p->asVector());
98
99 for (int i = 0; i < vi.size(); ++i)
100 {
101 const POINT& V = vi[i];
102 const POINT& A = vi[i-1];
103 const POINT& B = vi[i+1];
104
105 if (((V - A) ^ (B - V)).z() > 0.0) // Convex vertex? (assume CCW winding)
106 {
107 double minDistQV2 = DBL_MAX; // shortest distance from test point to candidate point
108 int indexBestQ; // index of best candidate point
109
110 // loop over all other verts (besides A, V, B):
111 for (int j = i + 2; j < i + 2 + vi.size() - 3; ++j)
112 {
113 const POINT& Q = vi[j];
114
115 if (triangleContains(A, V, B, Q))
116 {
117 double distQV2 = (Q - V).length2();
118 if (distQV2 < minDistQV2)
119 {
120 minDistQV2 = distQV2;
121 indexBestQ = j;
122 }
123 }
124 }
125
126 POINT result;
127
128 // If no inside point was found, return the midpoint of AB.
129 if (minDistQV2 == DBL_MAX)
130 {
131 result = (A+B)*0.5;
132 }
133
134 // Otherwise, use the midpoint of QV.
135 else
136 {
137 const POINT& Q = vi[indexBestQ];
138 result = (Q+V)*0.5;
139 }
140
141 // make sure the resulting point doesn't fall within any of the
142 // polygon's holes.
143 if (p->contains2D(result.x(), result.y()))
144 {
145 return result;
146 }
147 }
148 }
149
150 // Will only happen is holes prevent us from finding an internal point.
151 OE_WARN << LC << "getInternalPoint failed miserably\n";
152 return p->getBounds().center();
153 }
154
getDistanceSquaredToClosestEdge(const osg::Vec3d & P,const Symbology::Polygon * poly)155 double getDistanceSquaredToClosestEdge(const osg::Vec3d& P, const Symbology::Polygon* poly)
156 {
157 double Dmin = DBL_MAX;
158 ConstSegmentIterator segIter(poly, true);
159 while (segIter.hasMore())
160 {
161 const Segment segment = segIter.next();
162 const POINT& A = segment.first;
163 const POINT& B = segment.second;
164 const VECTOR AP = P-A, AB = B-A;
165 double t = clamp((AP*AB)/AB.length2(), 0.0, 1.0);
166 VECTOR PROJ = A + AB*t;
167 double D = (P - PROJ).length2();
168 if (D < Dmin) Dmin = D;
169 }
170 return Dmin;
171 }
172
173
174 struct Widths {
Widths__anon71f210340111::Widths175 Widths(const Widths& rhs) {
176 bufferWidth = rhs.bufferWidth;
177 lineWidth = rhs.lineWidth;
178 }
179
Widths__anon71f210340111::Widths180 Widths(double bufferWidth, double lineWidth) {
181 this->bufferWidth = bufferWidth;
182 this->lineWidth = lineWidth;
183 }
184
185 double bufferWidth;
186 double lineWidth;
187 };
188
189 typedef std::vector<Widths> WidthsList;
190
191 // Creates a heightfield that flattens an area intersecting the input polygon geometry.
192 // The height of the area is found by sampling a point internal to the polygon.
193 // bufferWidth = width of transition from flat area to natural terrain.
integratePolygons(const TileKey & key,osg::HeightField * hf,const MultiGeometry * geom,const SpatialReference * geomSRS,WidthsList & widths,ElevationEnvelope * envelope,bool fillAllPixels,ProgressCallback * progress)194 bool integratePolygons(const TileKey& key, osg::HeightField* hf, const MultiGeometry* geom, const SpatialReference* geomSRS,
195 WidthsList& widths, ElevationEnvelope* envelope,
196 bool fillAllPixels, ProgressCallback* progress)
197 {
198 bool wroteChanges = false;
199
200 const GeoExtent& ex = key.getExtent();
201
202 double col_interval = ex.width() / (double)(hf->getNumColumns()-1);
203 double row_interval = ex.height() / (double)(hf->getNumRows()-1);
204
205 POINT Pex, P, internalP;
206
207 bool needsTransform = ex.getSRS() != geomSRS;
208
209 for (unsigned col = 0; col < hf->getNumColumns(); ++col)
210 {
211 Pex.x() = ex.xMin() + (double)col * col_interval;
212
213 for (unsigned row = 0; row < hf->getNumRows(); ++row)
214 {
215 // check for cancelation periodically
216 //if (progress && progress->isCanceled())
217 // return false;
218
219 Pex.y() = ex.yMin() + (double)row * row_interval;
220
221 if (needsTransform)
222 ex.getSRS()->transform(Pex, geomSRS, P);
223 else
224 P = Pex;
225
226 bool done = false;
227 double minD2 = DBL_MAX;//bufferWidth * bufferWidth; // minimum distance(squared) to closest polygon edge
228 double bufferWidth = 0.0;
229
230 const Symbology::Polygon* bestPoly = 0L;
231
232 for (unsigned int geomIndex = 0; geomIndex < geom->getNumComponents(); geomIndex++)
233 {
234 Geometry* component = geom->getComponents()[geomIndex].get();
235 Widths width = widths[geomIndex];
236 ConstGeometryIterator giter(component, false);
237 while (giter.hasMore() && !done)
238 {
239 const Symbology::Polygon* polygon = dynamic_cast<const Symbology::Polygon*>(giter.next());
240 if (polygon)
241 {
242 // Does the point P fall within the polygon?
243 if (polygon->contains2D(P.x(), P.y()))
244 {
245 // yes, flatten it to the polygon's centroid elevation;
246 // and we're dont with this point.
247 done = true;
248 bestPoly = polygon;
249 minD2 = -1.0;
250 bufferWidth = width.bufferWidth;
251 }
252
253 // If not in the polygon, how far to the closest edge?
254 else
255 {
256 double D2 = getDistanceSquaredToClosestEdge(P, polygon);
257 if (D2 < minD2)
258 {
259 minD2 = D2;
260 bestPoly = polygon;
261 bufferWidth = width.bufferWidth;
262 }
263 }
264 }
265 }
266 }
267
268 if (bestPoly && minD2 != 0.0)
269 {
270 float h;
271 POINT internalP = getInternalPoint(bestPoly);
272 float elevInternal = envelope->getElevation(internalP.x(), internalP.y());
273
274 if (minD2 < 0.0)
275 {
276 h = elevInternal;
277 }
278 else
279 {
280 float elevNatural = envelope->getElevation(P.x(), P.y());
281 double blend = clamp(sqrt(minD2)/bufferWidth, 0.0, 1.0); // [0..1] 0=internal, 1=natural
282 h = smootherstep(elevInternal, elevNatural, blend);
283 }
284
285 hf->setHeight(col, row, h);
286 wroteChanges = true;
287 }
288
289 else if (fillAllPixels)
290 {
291 float h = envelope->getElevation(P.x(), P.y());
292 hf->setHeight(col, row, h);
293 // do not set wroteChanges
294 }
295 }
296 }
297
298 return wroteChanges;
299 }
300
301
302
303 struct Sample {
304 double D2; // distance to segment squared
305 osg::Vec3d A; // endpoint of segment
306 osg::Vec3d B; // other endpoint of segment;
307 double T; // segment parameter of closest point
308
309 // used later:
310 float elevPROJ; // elevation at point on segment
311 float elev; // flattened elevation
312 double D; // distance
313
314 double innerRadius;
315 double outerRadius;
316
operator <__anon71f210340111::Sample317 bool operator < (struct Sample& rhs) const { return D2 < rhs.D2; }
318 };
319
320 typedef std::vector<Sample> Samples;
321
EQ2(const osg::Vec3d & a,const osg::Vec3d & b)322 bool EQ2(const osg::Vec3d& a, const osg::Vec3d& b) {
323 return osg::equivalent(a.x(), b.x()) && osg::equivalent(a.y(), b.y());
324 }
325
isSampleValid(const Sample * b1,Samples & samples)326 bool isSampleValid(const Sample* b1, Samples& samples)
327 {
328 for (unsigned i = 0; i < samples.size(); ++i)
329 {
330 Sample* b2 = &samples[i];
331 if (b1 == b2) continue;
332 if (b1->T == 0.0 && (EQ2(b1->A, b2->A) || EQ2(b1->A, b2->B))) return false;
333 if (b1->T == 1.0 && (EQ2(b1->B, b2->A) || EQ2(b1->B, b2->B))) return false;
334 }
335 return true;
336 }
337
338 // Inverse Direct Weighting (IDW) interpolation
interpolateSamplesIDW(Samples & samples)339 float interpolateSamplesIDW(Samples& samples)
340 {
341 // higher powerParam corresponds to increased proximity favoring
342 const double powerParam = 2.5;
343
344 if (samples.size() == 0) return 0.0f;
345 if (samples.size() == 1) return samples[0].elev;
346
347 double numer = 0.0;
348 double denom = 0.0;
349 for (unsigned i = 0; i < samples.size(); ++i)
350 {
351 if (osg::equivalent(samples[i].D, 0.0)) {
352 numer = samples[i].elev, denom = 1.0;
353 break;
354 }
355 else {
356 double w = pow(1.0/samples[i].D, powerParam);
357 numer += w * samples[i].elev;
358 denom += w;
359 }
360 }
361
362 return numer / denom;
363 }
364
365 // Linear interpolation (simple mean)
interpolateSamplesLinear(Samples & samples)366 float interpolateSamplesLinear(Samples& samples)
367 {
368 double numer = 0.0;
369 for (unsigned i = 0; i<samples.size(); ++i)
370 numer += samples[i].elev;
371 return samples.size() > 0 ? (numer / (double)(samples.size())) : FLT_MAX;
372 }
373
374 /**
375 * Create a heightfield that flattens the terrain around linear geometry.
376 * lineWidth = width of completely flat area
377 * bufferWidth = width of transition from flat area to natural terrain
378 *
379 * Note: this algorithm only samples elevation data from the source (elevation pool).
380 * As it progresses, however, it is creating new modified elevation data -- but later
381 * points will continue to derive their source data from the original data. This means
382 * there will be some discontinuities in the final data, especially along the edges of
383 * the flattening buffer.
384 *
385 * There is not perfect solution for this, but one improvement would be to copy the
386 * source elevation into the heightfield as a starting point, and then sample that
387 * modifiable heightfield as we go along.
388 */
integrateLines(const TileKey & key,osg::HeightField * hf,const MultiGeometry * geom,const SpatialReference * geomSRS,WidthsList & widths,ElevationEnvelope * envelope,bool fillAllPixels,ProgressCallback * progress)389 bool integrateLines(const TileKey& key, osg::HeightField* hf, const MultiGeometry* geom, const SpatialReference* geomSRS,
390 WidthsList& widths, ElevationEnvelope* envelope,
391 bool fillAllPixels, ProgressCallback* progress)
392 {
393 bool wroteChanges = false;
394
395 const GeoExtent& ex = key.getExtent();
396
397 double col_interval = ex.width() / (double)(hf->getNumColumns()-1);
398 double row_interval = ex.height() / (double)(hf->getNumRows()-1);
399
400 osg::Vec3d Pex, P, PROJ;
401
402 bool needsTransform = ex.getSRS() != geomSRS;
403
404 // Loop over the new heightfield.
405 for (unsigned col = 0; col < hf->getNumColumns(); ++col)
406 {
407 Pex.x() = ex.xMin() + (double)col * col_interval;
408
409 for (unsigned row = 0; row < hf->getNumRows(); ++row)
410 {
411 // check for cancelation periodically
412 //if (progress && progress->isCanceled())
413 // return false;
414
415 Pex.y() = ex.yMin() + (double)row * row_interval;
416
417 // Move the point into the working SRS if necessary
418 if (needsTransform)
419 ex.getSRS()->transform(Pex, geomSRS, P);
420 else
421 P = Pex;
422
423 // For each point, we need to find the closest line segments to that point
424 // because the elevation values on these line segments will be the flattening
425 // value. There may be more than one line segment that falls within the search
426 // radius; we will collect up to MaxSamples of these for each heightfield point.
427 static const unsigned Maxsamples = 4;
428 Samples samples;
429
430 for (unsigned int geomIndex = 0; geomIndex < geom->getNumComponents(); geomIndex++)
431 {
432 Widths w = widths[geomIndex];
433 double innerRadius = w.lineWidth * 0.5;
434 double outerRadius = innerRadius + w.bufferWidth;
435 double outerRadius2 = outerRadius * outerRadius;
436
437 Geometry* component = geom->getComponents()[geomIndex].get();
438 // Search for line segments.
439 ConstGeometryIterator giter(component);
440 while (giter.hasMore())
441 {
442 const Geometry* part = giter.next();
443
444 for (int i = 0; i < part->size()-1; ++i)
445 {
446 // AB is a candidate line segment:
447 const osg::Vec3d& A = (*part)[i];
448 const osg::Vec3d& B = (*part)[i+1];
449
450 osg::Vec3d AB = B - A; // current segment AB
451
452 double t; // parameter [0..1] on segment AB
453 double D2; // shortest distance from point P to segment AB, squared
454 double L2 = AB.length2(); // length (squared) of segment AB
455 osg::Vec3d AP = P - A; // vector from endpoint A to point P
456
457 if (L2 == 0.0)
458 {
459 // trivial case: zero-length segment
460 t = 0.0;
461 D2 = AP.length2();
462 }
463 else
464 {
465 // Calculate parameter "t" [0..1] which will yield the closest point on AB to P.
466 // Clamping it means the closest point won't be beyond the endpoints of the segment.
467 t = clamp((AP * AB)/L2, 0.0, 1.0);
468
469 // project our point P onto segment AB:
470 PROJ.set( A + AB*t );
471
472 // measure the distance (squared) from P to the projected point on AB:
473 D2 = (P - PROJ).length2();
474 }
475
476 // If the distance from our point to the line segment falls within
477 // the maximum flattening distance, store it.
478 if (D2 <= outerRadius2)
479 {
480 // see if P is a new sample.
481 Sample* b;
482 if (samples.size() < Maxsamples)
483 {
484 // If we haven't collected the maximum number of samples yet,
485 // just add this to the list:
486 samples.push_back(Sample());
487 b = &samples.back();
488 }
489 else
490 {
491 // If we are maxed out on samples, find the farthest one we have so far
492 // and replace it if the new point is closer:
493 unsigned max_i = 0;
494 for (unsigned i=1; i<samples.size(); ++i)
495 if (samples[i].D2 > samples[max_i].D2)
496 max_i = i;
497
498 b = &samples[max_i];
499
500 if (b->D2 < D2)
501 b = 0L;
502 }
503
504 if (b)
505 {
506 b->D2 = D2;
507 b->A = A;
508 b->B = B;
509 b->T = t;
510 b->innerRadius = innerRadius;
511 b->outerRadius = outerRadius;
512 }
513 }
514 }
515 }
516 }
517
518 // Remove unnecessary sample points that lie on the endpoint of a segment
519 // that abuts another segment in our list.
520 for (unsigned i = 0; i < samples.size();) {
521 if (!isSampleValid(&samples[i], samples)) {
522 samples[i] = samples[samples.size() - 1];
523 samples.resize(samples.size() - 1);
524 }
525 else ++i;
526 }
527
528 // Now that we are done searching for line segments close to our point,
529 // we will collect the elevations at our sample points and use them to
530 // create a new elevation value for our point.
531 if (samples.size() > 0)
532 {
533 // The original elevation at our point:
534 float elevP = envelope->getElevation(P.x(), P.y());
535
536 for (unsigned i = 0; i < samples.size(); ++i)
537 {
538 Sample& sample = samples[i];
539
540 sample.D = sqrt(sample.D2);
541
542 // Blend factor. 0 = distance is less than or equal to the inner radius;
543 // 1 = distance is greater than or equal to the outer radius.
544 double blend = clamp(
545 (sample.D - sample.innerRadius) / (sample.outerRadius - sample.innerRadius),
546 0.0, 1.0);
547
548 if (sample.T == 0.0)
549 {
550 sample.elevPROJ = envelope->getElevation(sample.A.x(), sample.A.y());
551 if (sample.elevPROJ == NO_DATA_VALUE)
552 sample.elevPROJ = elevP;
553 }
554 else if (sample.T == 1.0)
555 {
556 sample.elevPROJ = envelope->getElevation(sample.B.x(), sample.B.y());
557 if (sample.elevPROJ == NO_DATA_VALUE)
558 sample.elevPROJ = elevP;
559 }
560 else
561 {
562 float elevA = envelope->getElevation(sample.A.x(), sample.A.y());
563 if (elevA == NO_DATA_VALUE)
564 elevA = elevP;
565
566 float elevB = envelope->getElevation(sample.B.x(), sample.B.y());
567 if (elevB == NO_DATA_VALUE)
568 elevB = elevP;
569
570 // linear interpolation of height from point A to point B on the segment:
571 sample.elevPROJ = mix(elevA, elevB, sample.T);
572 }
573
574 // smoothstep interpolation of along the buffer (perpendicular to the segment)
575 // will gently integrate the new value into the existing terrain.
576 sample.elev = smootherstep(sample.elevPROJ, elevP, blend);
577 }
578
579 // Finally, combine our new elevation values and set the new value in the output.
580 float finalElev = interpolateSamplesIDW(samples);
581 if (finalElev < FLT_MAX)
582 hf->setHeight(col, row, finalElev);
583 else
584 hf->setHeight(col, row, elevP);
585
586 wroteChanges = true;
587 }
588
589 else if (fillAllPixels)
590 {
591 // No close segments were found, so just copy over the source data.
592 float h = envelope->getElevation(P.x(), P.y());
593 hf->setHeight(col, row, h);
594
595 // Note: do not set wroteChanges to true.
596 }
597 }
598 }
599
600 return wroteChanges;
601 }
602
603
integrate(const TileKey & key,osg::HeightField * hf,const MultiGeometry * geom,const SpatialReference * geomSRS,WidthsList & widths,ElevationEnvelope * envelope,bool fillAllPixels,ProgressCallback * progress)604 bool integrate(const TileKey& key, osg::HeightField* hf, const MultiGeometry* geom, const SpatialReference* geomSRS,
605 WidthsList& widths, ElevationEnvelope* envelope,
606 bool fillAllPixels, ProgressCallback* progress)
607 {
608 if (geom->isLinear())
609 return integrateLines(key, hf, geom, geomSRS, widths, envelope, fillAllPixels, progress);
610 else
611 return integratePolygons(key, hf, geom, geomSRS, widths, envelope, fillAllPixels, progress);
612 }
613 }
614
615 //........................................................................
616
617
FlatteningLayer(const FlatteningLayerOptions & options)618 FlatteningLayer::FlatteningLayer(const FlatteningLayerOptions& options) :
619 ElevationLayer(&_optionsConcrete),
620 _options(&_optionsConcrete),
621 _optionsConcrete(options)
622 {
623 // Experiment with this and see what will work.
624 _pool = new ElevationPool();
625 _pool->setTileSize(257u);
626
627 init();
628 }
629
630 void
init()631 FlatteningLayer::init()
632 {
633 setTileSourceExpected(false);
634 ElevationLayer::init();
635 }
636
637 const Status&
open()638 FlatteningLayer::open()
639 {
640 // ensure the caller named a feature source:
641 if (!options().featureSource().isSet() &&
642 !options().featureSourceLayer().isSet())
643 {
644 return setStatus(Status::Error(Status::ConfigurationError, "Missing required feature source"));
645 }
646
647 // If the feature source is inline, open it now.
648 if (options().featureSource().isSet())
649 {
650 // Create the feature source instance:
651 FeatureSource* fs = FeatureSourceFactory::create(options().featureSource().get());
652 if (!fs)
653 {
654 return setStatus(Status::Error(Status::ServiceUnavailable, "Unable to create feature source as defined"));
655 }
656
657 // Open it:
658 setStatus(fs->open(getReadOptions()));
659 if (getStatus().isError())
660 return getStatus();
661
662 setFeatureSource(fs);
663
664 if (getStatus().isError())
665 return getStatus();
666 }
667
668 const Profile* profile = getProfile();
669 if ( !profile )
670 {
671 profile = Registry::instance()->getGlobalGeodeticProfile();
672 setProfile( profile );
673 }
674
675 return ElevationLayer::open();
676 }
677
678 void
setFeatureSource(FeatureSource * fs)679 FlatteningLayer::setFeatureSource(FeatureSource* fs)
680 {
681 if (fs)
682 {
683 _featureSource = fs;
684 if (_featureSource)
685 {
686 if (!_featureSource->getFeatureProfile())
687 {
688 setStatus(Status::Error(Status::ConfigurationError, "No feature profile (is the source open?)"));
689 _featureSource = 0L;
690 return;
691 }
692 }
693 }
694 }
695
~FlatteningLayer()696 FlatteningLayer::~FlatteningLayer()
697 {
698 _featureLayerListener.clear();
699 }
700
701 void
setFeatureSourceLayer(FeatureSourceLayer * layer)702 FlatteningLayer::setFeatureSourceLayer(FeatureSourceLayer* layer)
703 {
704 if (layer)
705 {
706 if (layer->getStatus().isOK())
707 setFeatureSource(layer->getFeatureSource());
708 }
709 else
710 {
711 setFeatureSource(0L);
712 }
713 }
714
715 void
addedToMap(const Map * map)716 FlatteningLayer::addedToMap(const Map* map)
717 {
718 // Initialize the elevation pool with our map:
719 OE_INFO << LC << "Attaching elevation pool to map\n";
720 _pool->setMap( map );
721
722
723 // Listen for our feature source layer to arrive, if there is one.
724 if (options().featureSourceLayer().isSet())
725 {
726 _featureLayerListener.listen(
727 map,
728 options().featureSourceLayer().get(),
729 this, &FlatteningLayer::setFeatureSourceLayer);
730 }
731
732 // Collect all elevation layers preceding this one and use them for flattening.
733 ElevationLayerVector layers;
734 map->getLayers(layers);
735 for (ElevationLayerVector::iterator i = layers.begin(); i != layers.end(); ++i) {
736 if (i->get() == this) {
737 layers.erase(i);
738 break;
739 }
740 else {
741 OE_INFO << LC << "Using: " << i->get()->getName() << "\n";
742 }
743 }
744 if (!layers.empty())
745 {
746 _pool->setElevationLayers(layers);
747 }
748 }
749
750 void
removedFromMap(const Map * map)751 FlatteningLayer::removedFromMap(const Map* map)
752 {
753 _featureLayerListener.clear();
754 }
755
756 void
createImplementation(const TileKey & key,osg::ref_ptr<osg::HeightField> & hf,osg::ref_ptr<NormalMap> & normalMap,ProgressCallback * progress)757 FlatteningLayer::createImplementation(const TileKey& key,
758 osg::ref_ptr<osg::HeightField>& hf,
759 osg::ref_ptr<NormalMap>& normalMap,
760 ProgressCallback* progress)
761 {
762 if (getStatus().isError())
763 {
764 return;
765 }
766
767 if (!_featureSource.valid())
768 {
769 setStatus(Status(Status::ServiceUnavailable, "No feature source"));
770 return;
771 }
772
773 const FeatureProfile* featureProfile = _featureSource->getFeatureProfile();
774 if (!featureProfile)
775 {
776 setStatus(Status(Status::ConfigurationError, "Feature profile is missing"));
777 return;
778 }
779
780 const SpatialReference* featureSRS = featureProfile->getSRS();
781 if (!featureSRS)
782 {
783 setStatus(Status(Status::ConfigurationError, "Feature profile has no SRS"));
784 return;
785 }
786
787 if (_pool->getElevationLayers().empty())
788 {
789 OE_WARN << LC << "Internal error - Pool layer set is empty\n";
790 return;
791 }
792
793 OE_START_TIMER(create);
794
795 // If the feature source has a tiling profile, we are going to have to map the incoming
796 // TileKey to a set of intersecting TileKeys in the feature source's tiling profile.
797 GeoExtent queryExtent = key.getExtent().transform(featureSRS);
798
799 // Lat/Long extent:
800 GeoExtent geoExtent = queryExtent.transform(featureSRS->getGeographicSRS());
801
802 // Buffer the query extent to include the potentially flattened area.
803 /*
804 double linewidth = SpatialReference::transformUnits(
805 options().lineWidth().get(),
806 featureSRS,
807 geoExtent.getCentroid().y());
808
809 double bufferwidth = SpatialReference::transformUnits(
810 options().bufferWidth().get(),
811 featureSRS,
812 geoExtent.getCentroid().y());
813 */
814 // TODO: JBFIX. Add a "max" setting somewhere.
815 double linewidth = 10.0;
816 double bufferwidth = 10.0;
817 double queryBuffer = 0.5*linewidth + bufferwidth;
818 queryExtent.expand(queryBuffer, queryBuffer);
819
820 // We must do all the feature processing in a projected system since we're using vector math.
821 #if 0
822 osg::ref_ptr<const SpatialReference> workingSRS = queryExtent.getSRS();
823 //if (workingSRS->isGeographic())
824 {
825 osg::Vec3d refPos = queryExtent.getCentroid();
826 workingSRS = workingSRS->createTangentPlaneSRS(refPos);
827 }
828 #else
829 const SpatialReference* workingSRS = queryExtent.getSRS()->isGeographic() ? SpatialReference::get("spherical-mercator") :
830 queryExtent.getSRS();
831 #endif
832
833 bool needsTransform = !featureSRS->isHorizEquivalentTo(workingSRS);
834
835 osg::ref_ptr< StyleSheet > styleSheet = new StyleSheet();
836 styleSheet->setScript(options().getScript());
837 osg::ref_ptr< Session > session = new Session( _map.get(), styleSheet.get());
838
839 // We will collection all the feature geometries in this multigeometry:
840 MultiGeometry geoms;
841 WidthsList widths;
842
843 if (featureProfile->getProfile())
844 {
845 // Tiled source, must resolve complete set of intersecting tiles:
846 std::vector<TileKey> intersectingKeys;
847 featureProfile->getProfile()->getIntersectingTiles(queryExtent, key.getLOD(), intersectingKeys);
848
849 std::set<TileKey> featureKeys;
850 for (int i = 0; i < intersectingKeys.size(); ++i)
851 {
852 if (intersectingKeys[i].getLOD() > featureProfile->getMaxLevel())
853 featureKeys.insert(intersectingKeys[i].createAncestorKey(featureProfile->getMaxLevel()));
854 else
855 featureKeys.insert(intersectingKeys[i]);
856 }
857
858 // Query and collect all the features we need for this tile.
859 for (std::set<TileKey>::const_iterator i = featureKeys.begin(); i != featureKeys.end(); ++i)
860 {
861 Query query;
862 query.tileKey() = *i;
863
864 osg::ref_ptr<FeatureCursor> cursor = _featureSource->createFeatureCursor(query, progress);
865 while (cursor.valid() && cursor->hasMore())
866 {
867 Feature* feature = cursor->nextFeature();
868
869 double lineWidth = 0.0;
870 double bufferWidth = 0.0;
871 if (options().lineWidth().isSet())
872 {
873 NumericExpression lineWidthExpr(options().lineWidth().get());
874 lineWidth = feature->eval(lineWidthExpr, session.get());
875 }
876
877 if (options().bufferWidth().isSet())
878 {
879 NumericExpression bufferWidthExpr(options().bufferWidth().get());
880 bufferWidth = feature->eval(bufferWidthExpr, session.get());
881 }
882
883 // Transform the feature geometry to our working (projected) SRS.
884 if (needsTransform)
885 feature->transform(workingSRS);
886
887 lineWidth = SpatialReference::transformUnits(
888 Distance(lineWidth),
889 featureSRS,
890 geoExtent.getCentroid().y());
891
892 bufferWidth = SpatialReference::transformUnits(
893 Distance(bufferWidth),
894 featureSRS,
895 geoExtent.getCentroid().y());
896
897 //TODO: optimization: test the geometry bounds against the expanded tilekey bounds
898 // in order to discard geometries we don't care about
899 geoms.getComponents().push_back(feature->getGeometry());
900 widths.push_back(Widths(bufferWidth, lineWidth));
901 }
902 }
903 }
904 else
905 {
906 // Non-tiled feaure source, just query arbitrary extent:
907 // Set up the query; bounds must be in the feature SRS:
908 Query query;
909 query.bounds() = queryExtent.bounds();
910
911 // Run the query and fill the list.
912 osg::ref_ptr<FeatureCursor> cursor = _featureSource->createFeatureCursor(query, progress);
913 while (cursor.valid() && cursor->hasMore())
914 {
915 Feature* feature = cursor->nextFeature();
916
917 double lineWidth = 0.0;
918 double bufferWidth = 0.0;
919 if (options().lineWidth().isSet())
920 {
921 NumericExpression lineWidthExpr(options().lineWidth().get());
922 lineWidth = feature->eval(lineWidthExpr, session.get());
923 }
924
925 if (options().bufferWidth().isSet())
926 {
927 NumericExpression bufferWidthExpr(options().bufferWidth().get());
928 bufferWidth = feature->eval(bufferWidthExpr, session.get());
929 }
930
931 // Transform the feature geometry to our working (projected) SRS.
932 if (needsTransform)
933 feature->transform(workingSRS);
934
935 lineWidth = SpatialReference::transformUnits(
936 Distance(lineWidth),
937 featureSRS,
938 geoExtent.getCentroid().y());
939
940 bufferWidth = SpatialReference::transformUnits(
941 Distance(bufferWidth),
942 featureSRS,
943 geoExtent.getCentroid().y());
944
945 // Transform the feature geometry to our working (projected) SRS.
946 if (needsTransform)
947 feature->transform(workingSRS);
948
949 //TODO: optimization: test the geometry bounds against the expanded tilekey bounds
950 // in order to discard geometries we don't care about
951
952 geoms.getComponents().push_back(feature->getGeometry());
953 widths.push_back(Widths(bufferWidth, lineWidth));
954 }
955 }
956
957 if (!geoms.getComponents().empty())
958 {
959 if (!hf.valid())
960 {
961 // Make an empty heightfield to populate:
962 hf = HeightFieldUtils::createReferenceHeightField(
963 queryExtent,
964 257, 257, // base tile size for elevation data
965 0u, // no border
966 true); // initialize to HAE (0.0) heights
967
968 // Initialize to NO DATA.
969 hf->getFloatArray()->assign(hf->getNumColumns()*hf->getNumRows(), NO_DATA_VALUE);
970 }
971
972 // Create an elevation query envelope at the LOD we are creating
973 osg::ref_ptr<ElevationEnvelope> envelope = _pool->createEnvelope(workingSRS, key.getLOD());
974
975 bool fill = (options().fill() == true);
976
977 integrate(key, hf.get(), &geoms, workingSRS, widths, envelope.get(), fill, progress);
978 }
979 }
980