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
20 #include <osgEarth/Cube>
21
22 using namespace osgEarth;
23
24 #define LC "[Cube] "
25
26 // --------------------------------------------------------------------------
27
28 bool
latLonToFaceCoords(double lat_deg,double lon_deg,double & out_x,double & out_y,int & out_face,int faceHint)29 CubeUtils::latLonToFaceCoords(double lat_deg, double lon_deg,
30 double& out_x, double& out_y, int& out_face,
31 int faceHint )
32 {
33 // normalized latitude and longitude
34 double nlat = (lat_deg+90.0)/180.0;
35 double nlon = (lon_deg+180.0)/360.0;
36
37 // check for out-of-range:
38 if ( nlat < 0 || nlat > 1 || nlon < 0 || nlon > 1 )
39 return false;
40
41 int face_x;
42
43 if ( faceHint >= 0 )
44 {
45 out_face = faceHint;
46 if ( faceHint < 4 )
47 {
48 face_x = faceHint;
49 }
50 else
51 {
52 face_x = (int)(4 * nlon);
53 if ( face_x == 4 )
54 face_x = 3;
55 }
56 }
57 else
58 {
59 face_x = (int)(4 * nlon);
60 if ( face_x == 4 )
61 face_x = 3;
62
63 int face_y = (int)(2 * nlat + 0.5);
64 if ( face_y == 1 )
65 out_face = face_x;
66 else
67 out_face = face_y < 1 ? 5 : 4;
68
69 //GW: not sure why this was here; but I think this issue is the cause of cracks when
70 // projecting CUBE source imagery onto a WGS84 globe.
71 //
72 //if ( osg::equivalent( lat_deg, -45 ) )
73 // out_face = 5;
74 }
75
76 out_x = 4 * nlon - face_x;
77 out_y = 2 * nlat - 0.5;
78
79 if(out_face < 4) // equatorial calculations done
80 return true;
81
82 double tmp;
83 if(out_face == 4) // north polar face
84 {
85 out_y = 1.5 - out_y;
86 out_x = 2 * (out_x - 0.5) * out_y + 0.5;
87 switch(face_x)
88 {
89 case 0: // bottom
90 out_y = 0.5 - out_y;
91 break;
92 case 1: // right side, swap and reverse lat
93 tmp = out_x;
94 out_x = 0.5 + out_y;
95 out_y = tmp;
96 break;
97 case 2: // top; reverse lat and lon
98 out_x = 1 - out_x;
99 out_y = 0.5 + out_y;
100 break;
101 case 3: // left side; swap and reverse lon
102 tmp = out_x;
103 out_x = 0.5 - out_y;
104 out_y = 1 - tmp;
105 break;
106 }
107 }
108 else // south polar face
109 {
110 out_y += 0.5;
111 out_x = 2 * (out_x - 0.5) * out_y + 0.5;
112 switch(face_x)
113 {
114 case 0: // left
115 tmp = out_x;
116 out_x = 0.5 - out_y;
117 out_y = tmp;
118 break;
119 case 1: // top
120 out_y = 0.5 + out_y;
121 break;
122 case 2: // right
123 tmp = out_x;
124 out_x = 0.5 + out_y;
125 out_y = 1 - tmp;
126 break;
127 case 3: // bottom
128 out_x = 1 - out_x;
129 out_y = 0.5 - out_y;
130 break;
131 }
132 }
133 return true;
134 }
135
136 bool
faceCoordsToLatLon(double x,double y,int face,double & out_lat_deg,double & out_lon_deg)137 CubeUtils::faceCoordsToLatLon( double x, double y, int face, double& out_lat_deg, double& out_lon_deg )
138 {
139 double offset = 0.0;
140 osg::Vec2d s( x, y );
141
142 // validate coordinate range:
143 if ( x < 0 || x > 1 || y < 0 || y > 1 )
144 {
145 OE_WARN << LC << "faceCoordToLatLon: input out of range" << std::endl;
146 return false;
147 }
148
149 if ( face < 4 ) // equatorial faces
150 {
151 s.x() = (x + face) * 0.25;
152 s.y() = (y + 0.5) * 0.5;
153 }
154 else if( face == 4 ) // north polar face
155 {
156 if ( x < y ) // left or top quadrant
157 {
158 if(x + y < 1.0) // left quadrant
159 {
160 s.x() = 1.0 - y;
161 s.y() = x;
162 offset += 3;
163 }
164 else // top quadrant
165 {
166 s.y() = 1.0 - y;
167 s.x() = 1.0 - x;
168 offset += 2;
169 }
170 }
171 else if( x + y >= 1.0 ) // right quadrant
172 {
173 s.x() = y;
174 s.y() = 1.0 - x;
175 offset += 1.0;
176 }
177 s.x() -= s.y();
178 if(s.y() != 0.5)
179 s.x() *= 0.5 / (0.5 - s.y());
180
181 s.x() = (s.x() + offset) * 0.25;
182 s.y() = (s.y() + 1.5) * 0.5;
183 }
184 else if ( face == 5 ) // south polar face
185 {
186 offset = 1.0;
187 if ( x > y ) // right or bottom quadrant
188 {
189 if( x + y >= 1.0) // right quadrant
190 {
191 s.x() = 1.0 - y;
192 s.y() = x - 0.5;
193 offset += 1.0;
194 }
195 else // bottom quadrant
196 {
197 s.x() = 1.0 - x;
198 s.y() = 0.5 - y;
199 offset += 2;
200 }
201 }
202 else // left or top quadrant
203 {
204 if(x + y < 1.0) // left quadrant
205 {
206 s.x() = y;
207 s.y() = 0.5 - x;
208 offset -= 1.0;
209 }
210 else // top quadrant
211 s.y() = y - 0.5;
212 }
213 if(s.y() != 0)
214 s.x() = (s.x() - 0.5) * 0.5 / s.y() + 0.5;
215 s.x() = (s.x() + offset) * 0.25;
216 s.y() *= 0.5;
217 }
218 else
219 {
220 return false; // invalid face specification
221 }
222
223 // convert to degrees
224 out_lon_deg = s.x() * 360 - 180;
225 out_lat_deg = s.y() * 180 - 90;
226
227 return true;
228 }
229
230 bool
cubeToFace(double & in_out_x,double & in_out_y,int & out_face)231 CubeUtils::cubeToFace( double& in_out_x, double& in_out_y, int& out_face )
232 {
233 // convert from unicube space (0,0=>6,1) to face space (0,0=>1,1 + face#)
234 // too tired to compute a formula right now
235 out_face =
236 in_out_x <= 1.0 ? 0 :
237 in_out_x <= 2.0 ? 1 :
238 in_out_x <= 3.0 ? 2 :
239 in_out_x <= 4.0 ? 3 :
240 in_out_x <= 5.0 ? 4 : 5;
241
242 in_out_x = in_out_x - (double)out_face;
243 // y unchanged
244 return true;
245 }
246
247 bool
cubeToFace(double & in_out_xmin,double & in_out_ymin,double & in_out_xmax,double & in_out_ymax,int & out_face)248 CubeUtils::cubeToFace(double& in_out_xmin, double& in_out_ymin,
249 double& in_out_xmax, double& in_out_ymax,
250 int& out_face)
251 {
252 int min_face =
253 in_out_xmin < 1.0 ? 0 :
254 in_out_xmin < 2.0 ? 1 :
255 in_out_xmin < 3.0 ? 2 :
256 in_out_xmin < 4.0 ? 3 :
257 in_out_xmin < 5.0 ? 4 : 5;
258
259 int max_face =
260 in_out_xmax <= 1.0 ? 0 :
261 in_out_xmax <= 2.0 ? 1 :
262 in_out_xmax <= 3.0 ? 2 :
263 in_out_xmax <= 4.0 ? 3 :
264 in_out_xmax <= 5.0 ? 4 : 5;
265
266 if ( min_face != max_face )
267 {
268 OE_WARN << LC << "Min face <> Max face!" << std::endl;
269 return false;
270 }
271
272 out_face = min_face;
273
274 in_out_xmin -= (double)out_face;
275 in_out_xmax -= (double)out_face;
276
277 // y values are unchanged
278 return true;
279 }
280
281 bool
faceToCube(double & in_out_x,double & in_out_y,int face)282 CubeUtils::faceToCube( double& in_out_x, double& in_out_y, int face )
283 {
284 // convert from face space (0,0=>1,1 + face#) to unicube space (0,0=>6,1)
285 in_out_x = (double)face + in_out_x;
286 // y unchanged
287 return true;
288 }
289
290 // --------------------------------------------------------------------------
291
CubeFaceLocator(unsigned int face)292 CubeFaceLocator::CubeFaceLocator(unsigned int face):
293 _face(face)
294 {
295 //NOP
296 }
297
~CubeFaceLocator()298 CubeFaceLocator::~CubeFaceLocator()
299 {
300 }
301
302
303 bool
convertLocalToModel(const osg::Vec3d & local,osg::Vec3d & world) const304 CubeFaceLocator::convertLocalToModel( const osg::Vec3d& local, osg::Vec3d& world ) const
305 {
306 #if ((OPENSCENEGRAPH_MAJOR_VERSION <= 2) && (OPENSCENEGRAPH_MINOR_VERSION < 8))
307 // OSG 2.7 bug workaround: bug fix in Locator submitted by GW
308 const_cast<CubeFaceLocator*>(this)->_inverse.invert( _transform );
309 #endif
310
311 if ( _coordinateSystemType == GEOCENTRIC )
312 {
313 //Convert the NDC coordinate into face space
314 osg::Vec3d faceCoord = local * _transform;
315
316 double lat_deg, lon_deg;
317 if ( !CubeUtils::faceCoordsToLatLon( faceCoord.x(), faceCoord.y(), _face, lat_deg, lon_deg ))
318 return false;
319
320 //OE_NOTICE << "LatLon=" << latLon << std::endl;
321
322 // convert to geocentric:
323 _ellipsoidModel->convertLatLongHeightToXYZ(
324 osg::DegreesToRadians( lat_deg ),
325 osg::DegreesToRadians( lon_deg ),
326 local.z(),
327 world.x(), world.y(), world.z() );
328
329 return true;
330 }
331 return true;
332 }
333
334
335 bool
convertModelToLocal(const osg::Vec3d & world,osg::Vec3d & local) const336 CubeFaceLocator::convertModelToLocal(const osg::Vec3d& world, osg::Vec3d& local) const
337 {
338 #if ((OPENSCENEGRAPH_MAJOR_VERSION <= 2) && (OPENSCENEGRAPH_MINOR_VERSION < 8))
339 // OSG 2.7 bug workaround: bug fix in Locator submitted by GW
340 const_cast<CubeFaceLocator*>(this)->_inverse.invert( _transform );
341 #endif
342
343 switch(_coordinateSystemType)
344 {
345 case(GEOCENTRIC):
346 {
347 double longitude, latitude, height;
348
349 _ellipsoidModel->convertXYZToLatLongHeight(world.x(), world.y(), world.z(), latitude, longitude, height );
350
351 int face=-1;
352 double x, y;
353
354 double lat_deg = osg::RadiansToDegrees(latitude);
355 double lon_deg = osg::RadiansToDegrees(longitude);
356
357 bool success = CubeUtils::latLonToFaceCoords( lat_deg, lon_deg, x, y, face, _face );
358
359 if (!success)
360 {
361 OE_WARN << LC << "Couldn't convert to face coords " << std::endl;
362 return false;
363 }
364 if (face != _face)
365 {
366 OE_WARN << LC
367 << "Face should be " << _face << " but is " << face
368 << ", lat = " << lat_deg
369 << ", lon = " << lon_deg
370 << std::endl;
371 }
372
373 local = osg::Vec3d( x, y, height ) * _inverse;
374 return true;
375 }
376
377
378 case(GEOGRAPHIC):
379 case(PROJECTED):
380 // Neither of these is supported for this locator..
381 {
382 local = world * _inverse;
383 return true;
384 }
385 }
386
387 return false;
388 }
389
390 // --------------------------------------------------------------------------
391
CubeSpatialReference(void * handle)392 CubeSpatialReference::CubeSpatialReference( void* handle ) :
393 SpatialReference(handle, std::string("OSGEARTH"))
394 {
395 _key.horiz = "unified-cube";
396 _key.horizLower = "unified-cube";
397 _name = "Unified Cube";
398 }
399
~CubeSpatialReference()400 CubeSpatialReference::~CubeSpatialReference()
401 {
402 }
403
404 void
_init()405 CubeSpatialReference::_init()
406 {
407 SpatialReference::_init();
408
409 _is_user_defined = true;
410 _is_cube = true;
411 _is_contiguous = false;
412 _is_geographic = false;
413 _key.horiz = "unified-cube";
414 _key.horizLower = "unified-cube";
415 _name = "Unified Cube";
416
417 // Custom units. The big number there roughly converts [0..1] to meters
418 // on a spheroid with WGS84-ish radius. Not perfect but close enough for
419 // the purposes of this class
420 _units = Units("Cube face", "cube", Units::TYPE_LINEAR, 42949672.96/4.0);
421 }
422
423 GeoLocator*
createLocator(double xmin,double ymin,double xmax,double ymax) const424 CubeSpatialReference::createLocator(double xmin, double ymin, double xmax, double ymax) const
425 {
426 int face;
427 CubeUtils::cubeToFace( xmin, ymin, xmax, ymax, face );
428
429 GeoLocator* result = new CubeFaceLocator( face );
430
431 osg::Matrixd transform;
432 transform.set(
433 xmax-xmin, 0.0, 0.0, 0.0,
434 0.0, ymax-ymin, 0.0, 0.0,
435 0.0, 0.0, 1.0, 0.0,
436 xmin, ymin, 0.0, 1.0);
437 result->setTransform( transform );
438
439 return result;
440 }
441
442 const SpatialReference*
preTransform(std::vector<osg::Vec3d> & points) const443 CubeSpatialReference::preTransform( std::vector<osg::Vec3d>& points ) const
444 {
445 for( unsigned i=0; i<points.size(); ++i )
446 {
447 osg::Vec3d& p = points[i];
448
449 // Convert the incoming points from cube => face => lat/long.
450 int face;
451 if ( !CubeUtils::cubeToFace( p.x(), p.y(), face ) )
452 {
453 OE_WARN << LC << "Failed to convert (" << p.x() << "," << p.y() << ") into face coordinates." << std::endl;
454 return 0L;
455 }
456
457 double lat_deg, lon_deg;
458 bool success = CubeUtils::faceCoordsToLatLon( p.x(), p.y(), face, lat_deg, lon_deg );
459 if (!success)
460 {
461 OE_WARN << LC <<
462 std::fixed << std::setprecision(2)
463 << "Could not transform face coordinates ["
464 << p.x() << ", " << p.y() << ", " << face << "] to lat lon"
465 << std::endl;
466 return 0L;
467 }
468 p.x() = lon_deg;
469 p.y() = lat_deg;
470 }
471 return getGeodeticSRS();
472 }
473
474 const SpatialReference*
postTransform(std::vector<osg::Vec3d> & points) const475 CubeSpatialReference::postTransform( std::vector<osg::Vec3d>& points) const
476 {
477 for( unsigned i=0; i<points.size(); ++i )
478 {
479 osg::Vec3d& p = points[i];
480
481 //Convert the incoming points from lat/lon back to face coordinates
482 int face;
483 double out_x, out_y;
484
485 // convert from lat/long to x/y/face
486 bool success = CubeUtils::latLonToFaceCoords( p.y(), p.x(), out_x, out_y, face );
487 if (!success)
488 {
489 OE_WARN << LC
490 << std::fixed << std::setprecision(2)
491 << "Could not transform lat long ["
492 << p.y() << ", " << p.x() << "] coordinates to face"
493 << std::endl;
494 return 0L;
495 }
496
497 //TODO: what to do about boundary points?
498
499 if ( !CubeUtils::faceToCube( out_x, out_y, face ) )
500 {
501 OE_WARN << LC << "fromFace(" << out_x << "," << out_y << "," << face << ") failed" << std::endl;
502 return 0L;
503 }
504
505 p.x() = out_x;
506 p.y() = out_y;
507 }
508 return getGeodeticSRS();
509 }
510
511 #define LL 0
512 #define LR 1
513 #define UR 2
514 #define UL 3
515 #define SMALLEST( W,X,Y,Z ) osg::minimum(W, osg::minimum( X, osg::minimum( Y, Z ) ) )
516 #define LARGEST( W,X,Y,Z ) osg::maximum(W, osg::maximum( X, osg::maximum( Y, Z ) ) )
517
518 bool
transformExtentToMBR(const SpatialReference * to_srs,double & in_out_xmin,double & in_out_ymin,double & in_out_xmax,double & in_out_ymax) const519 CubeSpatialReference::transformExtentToMBR(const SpatialReference* to_srs,
520 double& in_out_xmin,
521 double& in_out_ymin,
522 double& in_out_xmax,
523 double& in_out_ymax ) const
524 {
525 // input bounds:
526 Bounds inBounds(in_out_xmin, in_out_ymin, in_out_xmax, in_out_ymax);
527
528 Bounds outBounds;
529
530 // for each CUBE face, find the intersection of the input bounds and that face.
531 for (int face = 0; face < 6; ++face)
532 {
533 Bounds faceBounds( (double)(face), 0.0, (double)(face+1), 1.0);
534
535 Bounds intersection = faceBounds.intersectionWith(inBounds);
536
537 // if they intersect (with a non-zero area; abutting doesn't count in this case)
538 // transform the intersection and include in the result.
539 if (intersection.isValid() && intersection.area2d() > 0.0)
540 {
541 double
542 xmin = intersection.xMin(), ymin = intersection.yMin(),
543 xmax = intersection.xMax(), ymax = intersection.yMax();
544
545 if (transformInFaceExtentToMBR(to_srs, face, xmin, ymin, xmax, ymax))
546 {
547 outBounds.expandBy(Bounds(xmin, ymin, xmax, ymax));
548 }
549 }
550 }
551
552 if (outBounds.valid())
553 {
554 in_out_xmin = outBounds.xMin();
555 in_out_ymin = outBounds.yMin();
556 in_out_xmax = outBounds.xMax();
557 in_out_ymax = outBounds.yMax();
558 return true;
559 }
560 else
561 {
562 return false;
563 }
564 }
565
566 bool
transformInFaceExtentToMBR(const SpatialReference * to_srs,int face,double & in_out_xmin,double & in_out_ymin,double & in_out_xmax,double & in_out_ymax) const567 CubeSpatialReference::transformInFaceExtentToMBR(const SpatialReference* to_srs,
568 int face,
569 double& in_out_xmin,
570 double& in_out_ymin,
571 double& in_out_xmax,
572 double& in_out_ymax ) const
573 {
574
575
576 // note: this method only works when the extent is isolated to one face of the cube. If you
577 // want to transform an artibrary extent, you need to break it up into separate extents for
578 // each cube face.
579 bool ok = true;
580
581 double face_xmin = in_out_xmin, face_ymin = in_out_ymin;
582 double face_xmax = in_out_xmax, face_ymax = in_out_ymax;
583
584 //int face;
585 CubeUtils::cubeToFace( face_xmin, face_ymin, face_xmax, face_ymax, face );
586
587 // for equatorial faces, the normal transformation process will suffice (since it will call into
588 // pre/postTransform).
589 if ( face < 4 )
590 {
591 ok = SpatialReference::transformExtentToMBR( to_srs, in_out_xmin, in_out_ymin, in_out_xmax, in_out_ymax );
592 }
593 else
594 {
595 // otherwise we are on one of the polar faces (4 or 5):
596
597 // four corners in face space:
598 double fx[4] = { face_xmin, face_xmax, face_xmax, face_xmin };
599 double fy[4] = { face_ymin, face_ymin, face_ymax, face_ymax };
600
601 bool crosses_pole = fx[LL] < 0.5 && fx[UR] > 0.5 && fy[LL] < 0.5 && fy[UR] > 0.5;
602
603 if ( crosses_pole ) // full x extent.
604 {
605 bool north = face == 4; // else south
606 osg::Vec3d output;
607
608 to_srs->getGeographicSRS()->transform( osg::Vec3d(-180.0, north? 45.0 : -90.0, 0), to_srs, output );
609 in_out_xmin = output.x();
610 in_out_ymin = output.y();
611
612 to_srs->getGeographicSRS()->transform( osg::Vec3d(180.0, north? 90.0 : -45.0, 0), to_srs, output );
613 in_out_xmax = output.x();
614 in_out_ymax = output.y();
615
616 //to_srs->getGeographicSRS()->transform2D( -180.0, north? 45.0 : -90.0, to_srs, in_out_xmin, in_out_ymin );
617 //to_srs->getGeographicSRS()->transform2D( 180.0, north? 90.0 : -45.0, to_srs, in_out_xmax, in_out_ymax );
618 }
619
620 else
621 {
622 double lat_deg[4];
623 double lon_deg[4];
624 double latmin, latmax, lonmin, lonmax;
625
626 for( int i=0; i<4; ++i )
627 {
628 CubeUtils::faceCoordsToLatLon( fx[i], fy[i], face, lat_deg[i], lon_deg[i] );
629 }
630
631 latmin = SMALLEST( lat_deg[0], lat_deg[1], lat_deg[2], lat_deg[3] );
632 latmax = LARGEST( lat_deg[0], lat_deg[1], lat_deg[2], lat_deg[3] );
633
634 // check to see whether the extent crosses the date line boundary. If so,
635 // make the UL corner the southwest and the LR corner the east.
636 bool crosses_date_line = fx[UL]+(1-fy[UL]) < 1.0 && (1-fx[LR])+fy[LR] < 1.0 && fx[LL]+fy[LL] < 1.0;
637 if ( crosses_date_line )
638 {
639 lonmin = lon_deg[UL];
640 lonmax = lon_deg[LR];
641 }
642 else
643 {
644 lonmin = SMALLEST( lon_deg[0], lon_deg[1], lon_deg[2], lon_deg[3] );
645 lonmax = LARGEST( lon_deg[0], lon_deg[1], lon_deg[2], lon_deg[3] );
646 }
647
648 if ( to_srs->isGeographic() )
649 {
650 in_out_xmin = lonmin;
651 in_out_xmax = lonmax;
652 in_out_ymin = latmin;
653 in_out_ymax = latmax;
654 }
655 else
656 {
657 osg::Vec3d output;
658
659 bool ok1 = transform( osg::Vec3d(lonmin, latmin, 0), to_srs, output );
660 if ( ok1 ) {
661 in_out_xmin = output.x();
662 in_out_ymin = output.y();
663 }
664 bool ok2 = transform( osg::Vec3d(lonmax, latmax, 0), to_srs, output );
665 if ( ok2 ) {
666 in_out_xmax = output.x();
667 in_out_ymax = output.y();
668 }
669
670 //bool ok1 = transform2D( lonmin, latmin, to_srs, in_out_xmin, in_out_ymin, context );
671 //bool ok2 = transform2D( lonmax, latmax, to_srs, in_out_xmax, in_out_ymax, context );
672 ok = ok1 && ok2;
673 }
674 }
675 }
676
677 return ok;
678 }
679
680 // --------------------------------------------------------------------------
681
UnifiedCubeProfile()682 UnifiedCubeProfile::UnifiedCubeProfile() :
683 Profile(SpatialReference::create( "unified-cube" ),
684 0.0, 0.0, 6.0, 1.0,
685 -180.0, -90.0, 180.0, 90.0,
686 6, 1 )
687
688 {
689 const SpatialReference* srs = getSRS()->getGeographicSRS();
690
691 // set up some constant extents
692 _faceExtent_gcs[0] = GeoExtent( srs, -180, -45, -90, 45 );
693 _faceExtent_gcs[1] = GeoExtent( srs, -90, -45, 0, 45 );
694 _faceExtent_gcs[2] = GeoExtent( srs, 0, -45, 90, 45 );
695 _faceExtent_gcs[3] = GeoExtent( srs, 90, -45, 180, 45 );
696 _faceExtent_gcs[4] = GeoExtent( srs, -180, 45, 180, 90 ); // north polar
697 _faceExtent_gcs[5] = GeoExtent( srs, -180, -90, 180, -45 ); // south polar
698 }
699
700 int
getFace(const TileKey & key)701 UnifiedCubeProfile::getFace( const TileKey& key )
702 {
703 return key.getTileX() >> key.getLevelOfDetail();
704 }
705
706 GeoExtent
transformGcsExtentOnFace(const GeoExtent & gcsExtent,int face) const707 UnifiedCubeProfile::transformGcsExtentOnFace( const GeoExtent& gcsExtent, int face ) const
708 {
709 if ( face < 4 )
710 {
711 const GeoExtent& fex = _faceExtent_gcs[face];
712
713 return GeoExtent(
714 getSRS(),
715 (double)face + (gcsExtent.xMin()-fex.xMin()) / fex.width(),
716 (gcsExtent.yMin()-fex.yMin()) / fex.height(),
717 (double)face + (gcsExtent.xMax()-fex.xMin()) / fex.width(),
718 (gcsExtent.yMax()-fex.yMin()) / fex.height() );
719 }
720 else
721 {
722 // transform all 4 corners; then do the min/max for x/y.
723 double lon[4] = { gcsExtent.xMin(), gcsExtent.xMax(), gcsExtent.xMax(), gcsExtent.xMin() };
724 double lat[4] = { gcsExtent.yMin(), gcsExtent.yMin(), gcsExtent.yMax(), gcsExtent.yMax() };
725 double x[4], y[4];
726
727 for( int i=0; i<4; ++i )
728 {
729 int dummy;
730 if ( ! CubeUtils::latLonToFaceCoords( lat[i], lon[i], x[i], y[i], dummy, face ) )
731 {
732 OE_WARN << LC << "transformGcsExtentOnFace, ll2fc failed" << std::endl;
733 }
734 }
735
736 double xmin = SMALLEST( x[0], x[1], x[2], x[3] );
737 double xmax = LARGEST( x[0], x[1], x[2], x[3] );
738 double ymin = SMALLEST( y[0], y[1], y[2], y[3] );
739 double ymax = LARGEST( y[0], y[1], y[2], y[3] );
740
741 CubeUtils::faceToCube( xmin, ymin, face );
742 CubeUtils::faceToCube( xmax, ymax, face );
743
744 return GeoExtent( getSRS(), xmin, ymin, xmax, ymax );
745 }
746 }
747
748 void
getIntersectingTiles(const GeoExtent & remoteExtent,unsigned localLOD,std::vector<TileKey> & out_intersectingKeys) const749 UnifiedCubeProfile::getIntersectingTiles(const GeoExtent& remoteExtent,
750 unsigned localLOD,
751 std::vector<TileKey>& out_intersectingKeys ) const
752 {
753 if ( getSRS()->isHorizEquivalentTo( remoteExtent.getSRS() ) )
754 {
755 addIntersectingTiles( remoteExtent, localLOD, out_intersectingKeys );
756 }
757 else
758 {
759 // the cube profile is non-contiguous. so there may be multiple local extents required
760 // to fully intersect the remote extent.
761
762 // first transform the remote extent to lat/long.
763 GeoExtent remoteExtent_gcs = remoteExtent.getSRS()->isGeographic()
764 ? remoteExtent
765 : remoteExtent.transform( remoteExtent.getSRS()->getGeographicSRS() );
766
767 // Chop the input extent into three separate extents: for the equatorial, north polar,
768 // and south polar tile regions.
769 for( int face=0; face<6; ++face )
770 {
771 GeoExtent partExtent_gcs = _faceExtent_gcs[face].intersectionSameSRS( remoteExtent_gcs );
772 if ( partExtent_gcs.isValid() )
773 {
774 GeoExtent partExtent = transformGcsExtentOnFace( partExtent_gcs, face );
775 addIntersectingTiles( partExtent, localLOD, out_intersectingKeys );
776 }
777 }
778 }
779 }
780
781 unsigned
getEquivalentLOD(const Profile * rhsProfile,unsigned rhsLOD) const782 UnifiedCubeProfile::getEquivalentLOD(const Profile* rhsProfile, unsigned rhsLOD) const
783 {
784 return rhsLOD;
785 }
786
787
~UnifiedCubeProfile()788 UnifiedCubeProfile::~UnifiedCubeProfile()
789 {
790 }
791