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