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 "tcPhotographyData.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 {
  48   addDataSet(new tcLandsatData("/home/james/cs/pro/data/LE71950281999206EDC01/"));
  49 }
  50
  51 /// Destructor.
  52 tcGlobe::~tcGlobe()
  53 {
  54 }
  55
  56 /*
  57  * Data sets
  58  */
  59
  60 /// Add a dataset to the globe.
  61 void tcGlobe::addDataSet(tcPhotographyData* data)
  62 {
  63   m_data.push_back(data);
  64 }
  65
  66 /*
  67  * Rendering
  68  */
  69
  70 /// Draw a line of latitude.
  71 void tcGlobe::drawLineOfLatitude(double latitude) const
  72 {
  73   double z  = sin(latitude) * m_meanRadius;
  74   double xy = cos(latitude) * m_meanRadius;
  75   glBegin(GL_LINE_LOOP);
  76   {
  77     for (int lon = 0; lon < 360; ++lon)
  78     {
  79       glVertex3d(xy*sin(M_PI/180*lon), xy*cos(M_PI/180*lon), z);
  80     }
  81   }
  82   glEnd();
  83 }
  84
  85 /// Render a cell.
  86 void tcGlobe::renderCell(tcObserver* const observer, const tcGeo& swCorner, const tcGeo& neCorner) const
  87 {
  88   // Find the square distances to each corner
  89   GLmat4d modelview;
  90   glGetv(GL_MODELVIEW_MATRIX, &modelview);
  91   tcGeo geoCorners[4] = {
  92     tcGeo(swCorner.lon(), swCorner.lat()),
  93     tcGeo(swCorner.lon(), neCorner.lat()),
  94     tcGeo(neCorner.lon(), neCorner.lat()),
  95     tcGeo(neCorner.lon(), swCorner.lat())
  96   };
  97   GLvec3d cartCorners[4];
  98   float toCorners[4];
  99   for (int i = 0; i < 4; ++i)
 100   {
 101     cartCorners[i] = (GLvec3d)geoCorners[i] * m_meanRadius;
 102     toCorners[i] = (modelview*(cartCorners[i], 1.0)).slice<0,3>().sqr();
 103     // Cull faces which are roughly backfacing
 104     if ((modelview*(cartCorners[i], 0.0))[2] <= 0.0)
 105     {
 106       return;
 107     }
 108   }
 109
 110   // Decide whether to subdivide
 111   float diagonal = (cartCorners[0]-cartCorners[2]).sqr()*4;
 112   bool subdivide = true;
 113   // If it is disproportionately tall, only subdivide horizontally
 114   bool tall = (cartCorners[1] - cartCorners[0]) > (cartCorners[3] - cartCorners[0]).sqr()*4.0;
 115   for (int i = 0; i < 4; ++i)
 116   {
 117     if (toCorners[i] > diagonal)
 118     {
 119       subdivide = false;
 120       break;
 121     }
 122   }
 123
 124   if (subdivide)
 125   {
 126     glColor3f(0.0f, 0.0f, 1.0f);
 127   }
 128   else
 129   {
 130     glColor3f(1.0f, 0.5f, 0.0f);
 131   }
 132
 133 #if 0
 134   glBegin(GL_LINE_LOOP);
 135   for (int i = 0; i < 4; ++i)
 136   {
 137     glVertex3(cartCorners[i]);
 138   }
 139   glEnd();
 140 #endif
 141
 142   if (subdivide)
 143   {
 144     if (tall)
 145     {
 146       // bottom
 147       renderCell(observer,  geoCorners[0],                        (geoCorners[3] + geoCorners[2]) * 0.5);
 148       // top
 149       renderCell(observer, (geoCorners[0] + geoCorners[1]) * 0.5,  geoCorners[2]);
 150     }
 151     else
 152     {
 153       // bottom left
 154       renderCell(observer,  geoCorners[0],                        (geoCorners[0] + geoCorners[2]) * 0.5);
 155       // bottom right
 156       renderCell(observer, (geoCorners[0] + geoCorners[3]) * 0.5, (geoCorners[3] + geoCorners[2]) * 0.5);
 157       // top right
 158       renderCell(observer, (geoCorners[0] + geoCorners[1]) * 0.5, (geoCorners[1] + geoCorners[2]) * 0.5);
 159       // top right
 160       renderCell(observer, (geoCorners[0] + geoCorners[2]) * 0.5,  geoCorners[2]);
 161     }
 162     return;
 163   }
 164
 165   // Render this cell
 166   const int resolution = 32;
 167   for (int i = 0; i < resolution; ++i)
 168   {
 169     glBegin(GL_TRIANGLE_STRIP);
 170     {
 171       for (int j = 0; j <= resolution; ++j)
 172       {
 173         for (int k = 0; k < 2; ++k)
 174         {
 175           double dlon = (neCorner.lon()-swCorner.lon()) / resolution;
 176           double dlat = (neCorner.lat()-swCorner.lat()) / resolution;
 177           tcGeo coord = swCorner + tcGeo((i+k) * dlon, j * dlat);
 178           bool accurate = true;
 179           double alt = 0.0;
 180           switch (m_elevationMode)
 181           {
 182             case NoElevation:
 183               {
 184                 // get accurate set
 185                 m_elevation->altitudeAt(coord, true, &accurate);
 186               }
 187               break;
 188             case RawElevation:
 189               {
 190                 alt = m_elevation->altitudeAt(coord, false, &accurate);
 191                 if (!accurate)
 192                 {
 193                   alt = 0.0;
 194                 }
 195               }
 196               break;
 197             case CorrectedElevation:
 198               {
 199                 // get accurate set
 200                 m_elevation->altitudeAt(coord, true, &accurate);
 201                 if (m_elevationCorrection <= 0.0f)
 202                 {
 203                   alt = m_elevation->altitudeAt(coord, false, 0);
 204                 }
 205                 else if (m_elevationCorrection >= 1.0f)
 206                 {
 207                   alt = m_elevation->altitudeAt(coord, true, 0);
 208                 }
 209                 else
 210                 {
 211                   double alt1 = m_elevation->altitudeAt(coord, false, 0);
 212                   double alt2 = m_elevation->altitudeAt(coord, true, 0);
 213                   alt = alt1*(1.0-m_elevationCorrection) + alt2*m_elevationCorrection;
 214                 }
 215               }
 216               break;
 217             default:
 218               break;
 219           }
 220
 221           // Get colour
 222           float intensity = 0.0;
 223           foreach (tcPhotographyData* data, m_data)
 224           {
 225             intensity += data->sampleIntensity(coord);
 226           }
 227
 228           if (!accurate)
 229           {
 230             glColor4f(intensity, intensity, intensity, 0.5f);
 231           }
 232           else
 233           {
 234             glColor4f(intensity, intensity, intensity, 1.0f);
 235             //glColor4f(1.0f, 1.0f-(float)alt/3278.0f, 0.0f, 1.0f);
 236           }
 237           GLvec3d dir = coord;
 238           double rad = m_meanRadius + alt;
 239           glVertex3(dir * rad);
 240         }
 241       }
 242     }
 243     glEnd();
 244   }
 245 }
 246
 247 /// Render from the POV of an observer.
 248 void tcGlobe::render(tcObserver* const observer)
 249 {
 250   /// @todo use a really simple fragment shader to cull backfacing lines
 251 #if 1
 252   // Lines of latitude
 253   glColor3f(0.0f, 1.0f, 0.0f);
 254   for (int lat = -75; lat <= 75; lat += 15)
 255   {
 256     if (lat != 0)
 257     {
 258       drawLineOfLatitude(M_PI/180*lat);
 259     }
 260   }
 261
 262   double tropic = (23.0 + 26.0/60 + 22.0/3600) * M_PI/180;
 263   // Equator
 264   glColor3f(1.0f, 0.0f, 0.0f);
 265   drawLineOfLatitude(0.0);
 266   // Tropics (Capricorn and Cancer)
 267   glColor3f(1.0f, 0.0f, 1.0f);
 268   drawLineOfLatitude(-tropic);
 269   glColor3f(1.0f, 1.0f, 0.0f);
 270   drawLineOfLatitude(+tropic);
 271   // Arctic and Antarctic Circles
 272   glColor3f(1.0f, 1.0f, 1.0f);
 273   drawLineOfLatitude(+M_PI/2 - tropic);
 274   drawLineOfLatitude(-M_PI/2 + tropic);
 275
 276   // Lines of longitude
 277   for (int lon = 0; lon < 360; lon += 15)
 278   {
 279     double x = sin(M_PI/180*lon) * m_meanRadius;
 280     double y = -cos(M_PI/180*lon) * m_meanRadius;
 281     int minLat = 15;
 282     int maxLat = 165;
 283     if (lon % 180 == 0)
 284     {
 285       glLineWidth(2.0f);
 286       glColor3f(1.0f, lon/180, 0.0f);
 287     }
 288     if (lon % 90 == 0)
 289     {
 290       minLat = 0;
 291       maxLat = 180;
 292     }
 293     glBegin(GL_LINE_STRIP);
 294     {
 295       for (int lat = minLat; lat <= maxLat; ++lat)
 296       {
 297         double z = cos(M_PI/180*lat) * m_meanRadius;
 298         double xy = sin(M_PI/180*lat);
 299         glVertex3d(xy*x, xy*y, z);
 300       }
 301     }
 302     glEnd();
 303     if (lon % 180 == 0)
 304     {
 305       glLineWidth(1.0f);
 306       glColor3f(0.0f, 1.0f, 0.0f);
 307     }
 308   }
 309 #endif
 310
 311   // Draw data diagramatically
 312   foreach (tcPhotographyData* data, m_data)
 313   {
 314     data->renderSchematic(m_meanRadius, observer);
 315   }
 316
 317   // Go through cells
 318   const int dlon = 30;
 319   const int dlat = 30;
 320   for (int lon = -180; lon < 180; lon += dlon)
 321   {
 322     for (int lat = -90; lat < 90; lat += dlat)
 323     {
 324       tcGeo sw(M_PI/180 * (lon), M_PI/180 * (lat));
 325       tcGeo ne(M_PI/180 * (lon+dlon), M_PI/180 * (lat+dlat));
 326       renderCell(observer, sw, ne);
 327     }
 328   }
 329 }
 330
 331 /// Set the elevation mode to render in.
 332 void tcGlobe::setElevationMode(ElevationMode mode)
 333 {
 334   m_elevationMode = mode;
 335 }
 336
 337 /// Set the level of correction to show.
 338 void tcGlobe::setElevationCorrection(float correction)
 339 {
 340   m_elevationCorrection = correction;
 341 }
 342
 343 /*
 344  * Accessors
 345  */
 346
 347 /// Get the mean radius.
 348 double tcGlobe::meanRadius() const
 349 {
 350   return m_meanRadius;
 351 }