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