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