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