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