geo/tcGlobe.cpp

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 /**
  21  * @file tcGlobe.cpp
  22  * @brief Manages data for a globe.
  23  */
  24
  25 #include "tcGlobe.h"
  26 #include "tcGeoImageData.h"
  27 #include "tcLandsatData.h"
  28 #include "tcProcessingData.h"
  29 #include "tcCustomSunDirection.h"
  30 #include "tcSrtmModel.h"
  31 #include <glMatrix.h>
  32 #include <glVector.h>
  33
  34 #include <GL/gl.h>
  35
  36 #include <cmath>
  37
  38 /*
  39  * Constructors + destructor
  40  */
  41
  42 /// Primary constructor.
  43 tcGlobe::tcGlobe(double meanRadius)
  44 : m_meanRadius(meanRadius)
  45 , m_elevation(new tcSrtmModel[2])
  46 , m_imagery()
  47 , m_elevationMode()
  48 , m_elevationInterpolation(1.0f)
  49 , m_colourCoding(NoColourCoding)
  50 , m_colourMapping()
  51 {
  52   for (int i = 0; i < 2; ++i)
  53   {
  54     m_elevationMode[i] = CorrectedElevation;
  55   }
  56   for (int i = 0; i < 3; ++i)
  57   {
  58     m_colourMapping[i][0] = 0;
  59     m_colourMapping[i][1] = 2-i;
  60   }
  61   for (int i = 3; i < 6; ++i)
  62   {
  63     m_colourMapping[i][0] = -1;
  64     m_colourMapping[i][1] = -1;
  65   }
  66   QList<tcShadowClassifyingData*> datasets;
  67   // this one is a bit rubbish, the sun is too high
  68   //datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950281999206EDC01/", m_elevation);
  69
  70   // high azimuth, low elevation
  71   datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282001307EDC00/", m_elevation);
  72   datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282002006EDC00/", m_elevation);
  73   //datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282002310EDC00/", m_elevation);
  74
  75   datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282000081EDC00/", m_elevation);
  76   //datasets << new tcLandsatData("/home/james/cs/pro/data/etp195r28_5t19900910/",  m_elevation);
  77
  78 // Relative to 0N, 0E
  79 #define NEW_DATASET(AZ, EL, TILT) \
  80   datasets << new tcCustomSunDirection(m_elevation+1, tcGeo((180.0-AZ)*M_PI/180, (EL)*M_PI/180), \
  81                                                       tcGeo(7.0*M_PI/180, (TILT)*M_PI/180), true)
  82 // Relative to our region of interest
  83 #define NEW_DATASET_REL(AZ, EL) \
  84   datasets << new tcCustomSunDirection(m_elevation+1, tcGeo((180.0-AZ)*M_PI/180, (EL)*M_PI/180), \
  85                                                       tcGeo(7.0*M_PI/180,   46.0*M_PI/180))
  86
  87   // Mustn't be co-planar
  88   for (int i = -10; i <= 10; i += 20)
  89   {
  90     for (int j = 30; j <= 150; j += 10)
  91     {
  92       //NEW_DATASET((double)90, (double)j, (double)i);
  93     }
  94   }
  95
  96   foreach (tcShadowClassifyingData* dataset, datasets)
  97   {
  98     addImagery(dataset);
  99   }
 100   addImagery(new tcProcessingData(datasets, m_elevation));
 101 }
 102
 103 /// Destructor.
 104 tcGlobe::~tcGlobe()
 105 {
 106   foreach (tcGeoImageData* data, m_imagery)
 107   {
 108     delete data;
 109   }
 110   delete [] m_elevation;
 111 }
 112
 113 /*
 114  * Data sets
 115  */
 116
 117 /// Add some imagery data.
 118 void tcGlobe::addImagery(tcGeoImageData* data)
 119 {
 120   m_imagery.push_back(data);
 121 }
 122
 123 /// Get the imagery data.
 124 const QList<tcGeoImageData*>& tcGlobe::imagery() const
 125 {
 126   return m_imagery;
 127 }
 128
 129 /// Get the elevation model.
 130 tcSrtmModel* tcGlobe::dem() const
 131 {
 132   return m_elevation;
 133 }
 134
 135 /*
 136  * Rendering
 137  */
 138
 139 /// Draw a line of latitude.
 140 void tcGlobe::drawLineOfLatitude(double latitude) const
 141 {
 142   double z  = sin(latitude) * m_meanRadius;
 143   double xy = cos(latitude) * m_meanRadius;
 144   glBegin(GL_LINE_LOOP);
 145   {
 146     for (int lon = 0; lon < 360; ++lon)
 147     {
 148       glVertex3d(xy*sin(M_PI/180*lon), xy*cos(M_PI/180*lon), z);
 149     }
 150   }
 151   glEnd();
 152 }
 153
 154 /// Render a cell.
 155 void tcGlobe::renderCell(tcObserver* const observer, const tcGeo& swCorner, const tcGeo& neCorner, int samples, bool normals,
 156                          bool northEdge, bool eastEdge, bool southEdge, bool westEdge) const
 157 {
 158   // Sample at a sensible level
 159   tcSrtmModel::RenderState state;
 160   m_elevation->sampleAlign(swCorner, neCorner, &state, samples);
 161
 162   if (state.moreAvailableLon || state.moreAvailableLat)
 163   {
 164     // Find the square distances to each corner
 165     GLmat4d modelview;
 166     glGetv(GL_MODELVIEW_MATRIX, &modelview);
 167     tcGeo geoCorners[4] = {
 168       tcGeo(swCorner.lon(), swCorner.lat()),
 169       tcGeo(swCorner.lon(), neCorner.lat()),
 170       tcGeo(neCorner.lon(), neCorner.lat()),
 171       tcGeo(neCorner.lon(), swCorner.lat())
 172     };
 173     GLvec3d cartCorners[4];
 174     float toCorners[4];
 175     double altitudeMean = m_meanRadius + altitudeAt((geoCorners[0] + geoCorners[1] + geoCorners[2] + geoCorners[3])/4);
 176     for (int i = 0; i < 4; ++i)
 177     {
 178       cartCorners[i] = (GLvec3d)geoCorners[i] * altitudeMean;
 179       toCorners[i] = (modelview*(cartCorners[i], 1.0)).slice<0,3>().sqr();
 180       // Cull faces which are roughly backfacing
 181 // actually lets not
 182 #if 0
 183       if ((modelview*(cartCorners[i], 0.0))[2] <= 0.0)
 184       {
 185         return;
 186       }
 187 #endif
 188     }
 189
 190     // Decide whether to subdivide
 191     float diagonal = ( (cartCorners[0]-cartCorners[2]).sqr()
 192                      + (cartCorners[3]-cartCorners[1]).sqr())/2*4;
 193     // If it is disproportionately tall, only subdivide horizontally
 194     bool tall = (cartCorners[1] - cartCorners[0]).sqr() > (cartCorners[3] - cartCorners[0]).sqr()*4.0
 195              || (cartCorners[2] - cartCorners[3]).sqr() > (cartCorners[2] - cartCorners[1]).sqr()*4.0;
 196     // If it is disproportionately wide, only subdivide vertically
 197     bool wide = (cartCorners[3] - cartCorners[0]).sqr() > (cartCorners[1] - cartCorners[0]).sqr()*4.0
 198              || (cartCorners[2] - cartCorners[1]).sqr() > (cartCorners[2] - cartCorners[3]).sqr()*4.0;
 199     bool subdivide = true;
 200     for (int i = 0; i < 4; ++i)
 201     {
 202       if (toCorners[i] > diagonal)
 203       {
 204         subdivide = false;
 205         break;
 206       }
 207     }
 208
 209     if (subdivide)
 210     {
 211       if (tall && !wide)
 212       {
 213         // bottom
 214         renderCell(observer,  geoCorners[0],                        (geoCorners[3] + geoCorners[2]) * 0.5, samples, normals,
 215                    false, eastEdge, southEdge, westEdge);
 216         // top
 217         renderCell(observer, (geoCorners[0] + geoCorners[1]) * 0.5,  geoCorners[2], samples, normals,
 218                    northEdge, eastEdge, false, westEdge);
 219       }
 220       else if (wide && !tall)
 221       {
 222         // left
 223         renderCell(observer,  geoCorners[0],                        (geoCorners[1] + geoCorners[2]) * 0.5, samples, normals,
 224                    northEdge, false, southEdge, westEdge);
 225         // right
 226         renderCell(observer, (geoCorners[0] + geoCorners[3]) * 0.5,  geoCorners[2], samples, normals,
 227                    northEdge, eastEdge, southEdge, false);
 228       }
 229       else
 230       {
 231         // bottom left
 232         renderCell(observer,  geoCorners[0],                        (geoCorners[0] + geoCorners[2]) * 0.5, samples, normals,
 233                    false, false, southEdge, westEdge);
 234         // bottom right
 235         renderCell(observer, (geoCorners[0] + geoCorners[3]) * 0.5, (geoCorners[3] + geoCorners[2]) * 0.5, samples, normals,
 236                    false, eastEdge, southEdge, false);
 237         // top left
 238         renderCell(observer, (geoCorners[0] + geoCorners[1]) * 0.5, (geoCorners[1] + geoCorners[2]) * 0.5, samples, normals,
 239                    northEdge, false, false, westEdge);
 240         // top right
 241         renderCell(observer, (geoCorners[0] + geoCorners[2]) * 0.5,  geoCorners[2], samples, normals,
 242                    northEdge, eastEdge, false, false);
 243       }
 244       return;
 245     }
 246   }
 247
 248 #if 0
 249   if (subdivide)
 250   {
 251     glColor3f(0.0f, 0.0f, 1.0f);
 252   }
 253   else
 254   {
 255     glColor3f(1.0f, 0.5f, 0.0f);
 256   }
 257
 258   glBegin(GL_LINE_LOOP);
 259   for (int i = 0; i < 4; ++i)
 260   {
 261     glVertex3(cartCorners[i]);
 262   }
 263   glEnd();
 264 #endif
 265
 266   #define EDGE_KERNEL_1 \
 267           glColor4f(0.5f, 0.3f, 0.2f, 1.0f); \
 268           glVertex3(dir * (m_meanRadius+alt)); \
 269           glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
 270           glVertex3(dir * m_meanRadius);
 271   #define EDGE_KERNEL_2 \
 272           glColor4f(0.5f, 0.5f, 1.0f, 1.0f); \
 273           glVertex3(dir * m_meanRadius); \
 274           glColor4f(0.5f, 0.5f, 1.0f, 1.0f); \
 275           glVertex3(dir * (m_meanRadius-(accurate ? 100.0 : 1000.0)));
 276   #define EDGE_KERNEL_3 \
 277           glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
 278           glVertex3(dir * (m_meanRadius-(accurate ? 100.0 : 1000.0))); \
 279           glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
 280           glVertex3(dir * (m_meanRadius-5000.0));
 281   #define EDGE_KERNEL_4 \
 282           glColor4f(0.5f, 0.5f, 0.5f, 1.0f); \
 283           glVertex3(dir * (m_meanRadius-5000.0)); \
 284           glColor4f(1.0f, 0.0f, 0.0f, 0.5f); \
 285           glVertex3(dir * (m_meanRadius-8000.0));
 286   #define EDGE_KERNEL(STAGE, EDGE, SAMPLES, LON, LAT) \
 287       if (EDGE##Edge) \
 288       { \
 289         glBegin(GL_TRIANGLE_STRIP); \
 290         for (int i = 0; i < SAMPLES; ++i) \
 291         { \
 292           tcGeo coord; \
 293           bool accurate = true; \
 294           double alt = altitudeAt(state, (LON), (LAT), &coord, &accurate); \
 295           GLvec3d dir = coord; \
 296           EDGE_KERNEL_##STAGE \
 297         } \
 298         glEnd(); \
 299       }
 300   #define EDGE(STAGE) \
 301       EDGE_KERNEL(STAGE, north, state.samplesLon, i,                  state.samplesLat-1); \
 302       EDGE_KERNEL(STAGE, east,  state.samplesLat, state.samplesLon-1, i); \
 303       EDGE_KERNEL(STAGE, south, state.samplesLon, i,                  0); \
 304       EDGE_KERNEL(STAGE, west,  state.samplesLat, 0,                  i);
 305
 306   if (normals)
 307   {
 308     // Draw normals
 309     float normColours[3][2][3] = {
 310       { { 1.0f, 1.0f, 1.0f }, { 1.0f, 0.0f, 0.0f } },
 311       { { 1.0f, 1.0f, 1.0f }, { 0.0f, 1.0f, 0.0f } },
 312       { { 1.0f, 1.0f, 1.0f }, { 0.0f, 0.0f, 1.0f } }
 313     };
 314     glBegin(GL_LINES);
 315     for (int i = 0; i < state.samplesLon; ++i)
 316     {
 317       for (int j = 0; j < state.samplesLat; ++j)
 318       {
 319         tcGeo coord;
 320         bool accurate = true;
 321         double alt = altitudeAt(state, i, j, &coord, &accurate);
 322         for (int n = 0; n < 3; ++n)
 323         {
 324           GLvec3f norm = normalAt(n, coord);
 325           if (!norm.zero())
 326           {
 327             norm = coord * norm;
 328             GLvec3d dir = coord;
 329             double rad = m_meanRadius + alt;
 330             glColor3fv(normColours[n][0]);
 331             glVertex3(dir * rad);
 332             glColor3fv(normColours[n][1]);
 333             glVertex3(dir * rad + norm*50.0f);
 334           }
 335         }
 336       }
 337     }
 338     glEnd();
 339   }
 340   else
 341   {
 342 #if 0
 343     // Render the solid rock walls
 344     glDisable(GL_CULL_FACE);
 345     EDGE(1)
 346     EDGE(2)
 347     EDGE(3)
 348     EDGE(4)
 349     glEnable(GL_CULL_FACE);
 350 #endif
 351
 352     // Render this cell
 353     /// @todo cache one edge of strip to save time on other
 354     for (int i = 0; i < state.samplesLon-1; ++i)
 355     {
 356       glBegin(GL_TRIANGLE_STRIP);
 357       {
 358         for (int j = 0; j < state.samplesLat; ++j)
 359         {
 360           for (int k = 0; k < 2; ++k)
 361           {
 362             tcGeo coord;
 363             bool accurate = true;
 364             double alt = altitudeAt(state, i+k, j, &coord, &accurate);
 365
 366             // Get colour
 367             foreach (tcGeoImageData* imagery, m_imagery)
 368             {
 369               imagery->texCoord(coord);
 370             }
 371
 372             if (!accurate)
 373             {
 374               glColor4f(0.5f, 0.5f, 1.0f, 1.0f);
 375             }
 376             else
 377             {
 378               glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
 379               //glColor4f(1.0f, 1.0f-(float)alt/3278.0f, 0.0f, 1.0f);
 380             }
 381             // Colour code if applicable
 382             if (m_colourCoding == ElevationSampleAlignment)
 383             {
 384               if (state.moreAvailableLon && state.moreAvailableLat)
 385               {
 386                 glColor3f(1.0f, 0.0f, 1.0f);
 387               }
 388               else if (state.moreAvailableLon)
 389               {
 390                 glColor3f(1.0f, 0.0f, 0.0f);
 391               }
 392               else if (state.moreAvailableLat)
 393               {
 394                 glColor3f(0.0f, 0.0f, 1.0f);
 395               }
 396               else
 397               {
 398                 glColor3f(0.0f, 1.0f, 0.0f);
 399               }
 400             }
 401             GLvec3d dir = coord;
 402             double rad = m_meanRadius + alt;
 403             glVertex3(dir * rad);
 404           }
 405         }
 406       }
 407       glEnd();
 408     }
 409   }
 410 }
 411
 412 /// Render from the POV of an observer.
 413 void tcGlobe::render(tcObserver* const observer, bool adaptive, const tcGeo* extent)
 414 {
 415   /// @todo use a really simple fragment shader to cull backfacing lines
 416 #if 1
 417   // Lines of latitude
 418   glColor3f(0.0f, 1.0f, 0.0f);
 419   for (int lat = -75; lat <= 75; lat += 15)
 420   {
 421     if (lat != 0)
 422     {
 423       drawLineOfLatitude(M_PI/180*lat);
 424     }
 425   }
 426
 427   double tropic = (23.0 + 26.0/60 + 22.0/3600) * M_PI/180;
 428   // Equator
 429   glColor3f(1.0f, 0.0f, 0.0f);
 430   drawLineOfLatitude(0.0);
 431   // Tropics (Capricorn and Cancer)
 432   glColor3f(1.0f, 0.0f, 1.0f);
 433   drawLineOfLatitude(-tropic);
 434   glColor3f(1.0f, 1.0f, 0.0f);
 435   drawLineOfLatitude(+tropic);
 436   // Arctic and Antarctic Circles
 437   glColor3f(1.0f, 1.0f, 1.0f);
 438   drawLineOfLatitude(+M_PI/2 - tropic);
 439   drawLineOfLatitude(-M_PI/2 + tropic);
 440
 441   // Lines of longitude
 442   for (int lon = 0; lon < 360; lon += 15)
 443   {
 444     double x = sin(M_PI/180*lon) * m_meanRadius;
 445     double y = -cos(M_PI/180*lon) * m_meanRadius;
 446     int minLat = 15;
 447     int maxLat = 165;
 448     if (lon % 180 == 0)
 449     {
 450       glLineWidth(2.0f);
 451       glColor3f(1.0f, lon/180, 0.0f);
 452     }
 453     if (lon % 90 == 0)
 454     {
 455       minLat = 0;
 456       maxLat = 180;
 457     }
 458     glBegin(GL_LINE_STRIP);
 459     {
 460       for (int lat = minLat; lat <= maxLat; ++lat)
 461       {
 462         double z = cos(M_PI/180*lat) * m_meanRadius;
 463         double xy = sin(M_PI/180*lat);
 464         glVertex3d(xy*x, xy*y, z);
 465       }
 466     }
 467     glEnd();
 468     if (lon % 180 == 0)
 469     {
 470       glLineWidth(1.0f);
 471       glColor3f(0.0f, 1.0f, 0.0f);
 472     }
 473   }
 474 #endif
 475
 476   // Draw data diagramatically
 477   foreach (tcGeoImageData* imagery, m_imagery)
 478   {
 479     imagery->renderSchematic(m_meanRadius, observer);
 480   }
 481
 482   if (0 == extent)
 483   {
 484     int colourMapping[6];
 485     for (int band = 0; band < m_imagery.size(); ++band)
 486     {
 487       for (int i = 0; i < 6; ++i)
 488       {
 489         colourMapping[i] = (m_colourMapping[i][0] == band ? m_colourMapping[i][1] : -1);
 490       }
 491       tcGeoImageData* imagery = m_imagery[band];
 492       imagery->setupThumbnailRendering(6, colourMapping);
 493     }
 494     // Go through cells
 495     const int dlon = 30;
 496     const int dlat = 30;
 497     for (int lon = -180; lon < 180; lon += dlon)
 498     {
 499       for (int lat = -90; lat < 90; lat += dlat)
 500       {
 501         tcGeo sw(M_PI/180 * (lon), M_PI/180 * (lat));
 502         tcGeo ne(M_PI/180 * (lon+dlon), M_PI/180 * (lat+dlat));
 503         renderCell(observer, sw, ne, adaptive ? 5 : 0, false);
 504       }
 505     }
 506     foreach (tcGeoImageData* imagery, m_imagery)
 507     {
 508       imagery->setupNormalRendering();
 509     }
 510     for (int lon = -180; lon < 180; lon += dlon)
 511     {
 512       for (int lat = -90; lat < 90; lat += dlat)
 513       {
 514         tcGeo sw(M_PI/180 * (lon), M_PI/180 * (lat));
 515         tcGeo ne(M_PI/180 * (lon+dlon), M_PI/180 * (lat+dlat));
 516         renderCell(observer, sw, ne, adaptive ? 5 : 0, true);
 517       }
 518     }
 519   }
 520   else
 521   {
 522     tcGeo sw = extent[0];
 523     tcGeo ne = extent[1];
 524     if (sw.lon() > ne.lon())
 525     {
 526       sw.setLon(ne.lon());
 527       ne.setLon(extent[0].lon());
 528     }
 529     if (sw.lat() > ne.lat())
 530     {
 531       sw.setLat(ne.lat());
 532       ne.setLat(extent[0].lat());
 533     }
 534     int colourMapping[6];
 535     for (int band = 0; band < m_imagery.size(); ++band)
 536     {
 537       for (int i = 0; i < 6; ++i)
 538       {
 539         colourMapping[i] = (m_colourMapping[i][0] == band ? m_colourMapping[i][1] : -1);
 540       }
 541       tcGeoImageData* imagery = m_imagery[band];
 542       imagery->setupDetailedRendering(6, colourMapping, sw, ne);
 543     }
 544     /// @todo If it is really big, split it
 545     renderCell(observer, sw, ne, adaptive ? 16 : 0, false, true, true, true, true);
 546     foreach (tcGeoImageData* imagery, m_imagery)
 547     {
 548       imagery->setupNormalRendering();
 549     }
 550     renderCell(observer, sw, ne, adaptive ? 16 : 0, true, true, true, true, true);
 551   }
 552   foreach (tcGeoImageData* imagery, m_imagery)
 553   {
 554     imagery->finishRendering();
 555   }
 556 }
 557
 558 /// Set the elevation mode to render in.
 559 void tcGlobe::setElevationMode(int dem, ElevationMode mode)
 560 {
 561   Q_ASSERT(dem >= 0 && dem < 2);
 562   m_elevationMode[dem] = mode;
 563 }
 564
 565 /// Set the level of interpolation.
 566 void tcGlobe::setElevationInterpolation(float interpolation)
 567 {
 568   m_elevationInterpolation = interpolation;
 569 }
 570
 571 /// Set the elevation data set name.
 572 void tcGlobe::setElevationDataSet(int dem, const QString& name)
 573 {
 574   Q_ASSERT(dem >= 0 && dem < 2);
 575   m_elevation[dem].setDataSet(name);
 576 }
 577
 578 /// Set colour coding method.
 579 void tcGlobe::setColourCoding(ColourCoding colourCoding)
 580 {
 581   m_colourCoding = colourCoding;
 582 }
 583
 584 /// Adjust the mapping between bands and colour channels.
 585 void tcGlobe::setColourMapping(int outputChannel, int inputBand, int inputGroup)
 586 {
 587   if (outputChannel < 6)
 588   {
 589     m_colourMapping[outputChannel][0] = inputGroup;
 590     m_colourMapping[outputChannel][1] = inputBand;
 591   }
 592 }
 593
 594 /*
 595  * Accessors
 596  */
 597
 598 /// Get the mean radius.
 599 double tcGlobe::meanRadius() const
 600 {
 601   return m_meanRadius;
 602 }
 603
 604 /// Get the altitude above sea level at a sample in a render state.
 605 double tcGlobe::altitudeAt(const tcSrtmModel::RenderState& state, int x, int y, tcGeo* outCoord, bool* isAccurate) const
 606 {
 607   bool accurate[2] = {
 608     true,
 609     true
 610   };
 611   double alt[2] = {
 612     0.0,
 613     0.0
 614   };
 615   bool needElev[2] = {
 616     m_elevationInterpolation < 1.0f,
 617     m_elevationInterpolation > 0.0f
 618   };
 619   for (int i = 0; i < 2; ++i)
 620   {
 621     switch (m_elevationMode[i])
 622     {
 623       case NoElevation:
 624         break;
 625       case RawElevation:
 626         {
 627           alt[i] = m_elevation[i].altitudeAt(state, x, y, outCoord, false, &accurate[i]);
 628           if (!accurate[i])
 629           {
 630             alt[i] = 0.0;
 631           }
 632         }
 633         break;
 634       case CorrectedElevation:
 635         {
 636           // get accurate set
 637           alt[i] = m_elevation[i].altitudeAt(state, x, y, outCoord, true, &accurate[i]);
 638         }
 639         break;
 640       default:
 641         Q_ASSERT(false);
 642         break;
 643     }
 644   }
 645   double result = alt[0]*(1.0-m_elevationInterpolation) + alt[1]*m_elevationInterpolation;
 646   if (0 != isAccurate)
 647   {
 648     *isAccurate = accurate[0] && accurate[1];
 649   }
 650   return result;
 651 }
 652
 653 /// Get the altitude above sea level at a coordinate.
 654 double tcGlobe::altitudeAt(const tcGeo& coord, bool* isAccurate) const
 655 {
 656   bool accurate[2] = {
 657     true,
 658     true
 659   };
 660   double alt[2] = {
 661     0.0,
 662     0.0
 663   };
 664   bool needElev[2] = {
 665     m_elevationInterpolation < 1.0f,
 666     m_elevationInterpolation > 0.0f
 667   };
 668   for (int i = 0; i < 2; ++i)
 669   {
 670     switch (m_elevationMode[i])
 671     {
 672       case NoElevation:
 673         break;
 674       case RawElevation:
 675         {
 676           alt[i] = m_elevation[i].altitudeAt(coord, false, &accurate[i]);
 677           if (!accurate[i])
 678           {
 679             alt[i] = 0.0;
 680           }
 681         }
 682         break;
 683       case CorrectedElevation:
 684         {
 685           // get accurate set
 686           alt[i] = m_elevation[i].altitudeAt(coord, true, &accurate[i]);
 687         }
 688         break;
 689       default:
 690         Q_ASSERT(false);
 691         break;
 692     }
 693   }
 694   double result = alt[0]*(1.0-m_elevationInterpolation) + alt[1]*m_elevationInterpolation;
 695   if (0 != isAccurate)
 696   {
 697     *isAccurate = accurate[0] && accurate[1];
 698   }
 699   return result;
 700 }
 701
 702 /// Get the radius at a coordinate.
 703 double tcGlobe::radiusAt(const tcGeo& coord) const
 704 {
 705   return m_meanRadius + altitudeAt(coord);
 706 }
 707
 708 /// Get the texture coordinate of the effective texture at a geographical coordinate.
 709 maths::Vector<2,double> tcGlobe::textureCoordOfGeo(const tcGeo& coord) const
 710 {
 711   /// @todo Reimplement tcGlobe::textureCoordOfGeo with multiple imagery
 712   Q_ASSERT(0 && "Reimplement tcGlobe::TextureCoordOfGeo with multiple imagery");
 713   return m_imagery[0]->geoToEffectiveTex() * coord;
 714 }
 715
 716 /// Get the current normal at a coordinate.
 717 maths::Vector<3,float> tcGlobe::normalAt(int norm, const tcGeo& coord) const
 718 {
 719   Q_ASSERT(norm >= 0 && norm < 3);
 720   if (m_colourMapping[3+norm][0] >= 0)
 721   {
 722     return m_imagery[m_colourMapping[3+norm][0]]->normalAt(norm, coord);
 723   }
 724   else
 725   {
 726     return maths::Vector<3,float>(0.0f);
 727   }
 728 }