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