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