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