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 = 32;
 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
 179           // Get colour
 180           float intensity = 0.0;
 181           foreach (tcPhotographyData* data, m_data)
 182           {
 183             intensity += data->sampleIntensity(coord);
 184           }
 185
 186           if (!accurate)
 187           {
 188             glColor4f(0.0f, 0.0f, 0.5f, 0.5f);
 189           }
 190           else
 191           {
 192             glColor4f(intensity, intensity, intensity, 1.0f);
 193             //glColor4f(1.0f, 1.0f-(float)alt/3278.0f, 0.0f, 1.0f);
 194           }
 195           GLvec3d dir = coord;
 196           double rad = m_meanRadius + alt;
 197           glVertex3(dir * rad);
 198         }
 199       }
 200     }
 201     glEnd();
 202   }
 203 }
 204
 205 /// Render from the POV of an observer.
 206 void tcGlobe::render(tcObserver* const observer)
 207 {
 208   /// @todo use a really simple fragment shader to cull backfacing lines
 209 #if 1
 210   // Lines of latitude
 211   glColor3f(0.0f, 1.0f, 0.0f);
 212   for (int lat = -75; lat <= 75; lat += 15)
 213   {
 214     if (lat != 0)
 215     {
 216       drawLineOfLatitude(M_PI/180*lat);
 217     }
 218   }
 219
 220   double tropic = (23.0 + 26.0/60 + 22.0/3600) * M_PI/180;
 221   // Equator
 222   glColor3f(1.0f, 0.0f, 0.0f);
 223   drawLineOfLatitude(0.0);
 224   // Tropics (Capricorn and Cancer)
 225   glColor3f(1.0f, 0.0f, 1.0f);
 226   drawLineOfLatitude(-tropic);
 227   glColor3f(1.0f, 1.0f, 0.0f);
 228   drawLineOfLatitude(+tropic);
 229   // Arctic and Antarctic Circles
 230   glColor3f(1.0f, 1.0f, 1.0f);
 231   drawLineOfLatitude(+M_PI/2 - tropic);
 232   drawLineOfLatitude(-M_PI/2 + tropic);
 233
 234   // Lines of longitude
 235   for (int lon = 0; lon < 360; lon += 15)
 236   {
 237     double x = sin(M_PI/180*lon) * m_meanRadius;
 238     double y = -cos(M_PI/180*lon) * m_meanRadius;
 239     int minLat = 15;
 240     int maxLat = 165;
 241     if (lon % 180 == 0)
 242     {
 243       glLineWidth(2.0f);
 244       glColor3f(1.0f, lon/180, 0.0f);
 245     }
 246     if (lon % 90 == 0)
 247     {
 248       minLat = 0;
 249       maxLat = 180;
 250     }
 251     glBegin(GL_LINE_STRIP);
 252     {
 253       for (int lat = minLat; lat <= maxLat; ++lat)
 254       {
 255         double z = cos(M_PI/180*lat) * m_meanRadius;
 256         double xy = sin(M_PI/180*lat);
 257         glVertex3d(xy*x, xy*y, z);
 258       }
 259     }
 260     glEnd();
 261     if (lon % 180 == 0)
 262     {
 263       glLineWidth(1.0f);
 264       glColor3f(0.0f, 1.0f, 0.0f);
 265     }
 266   }
 267 #endif
 268
 269   // Draw data diagramatically
 270   foreach (tcPhotographyData* data, m_data)
 271   {
 272     data->renderSchematic(m_meanRadius, observer);
 273   }
 274
 275   // Go through cells
 276   const int dlon = 30;
 277   const int dlat = 30;
 278   for (int lon = -180; lon < 180; lon += dlon)
 279   {
 280     for (int lat = -90; lat < 90; lat += dlat)
 281     {
 282       tcGeo sw(M_PI/180 * (lon), M_PI/180 * (lat));
 283       tcGeo ne(M_PI/180 * (lon+dlon), M_PI/180 * (lat+dlat));
 284       renderCell(observer, sw, ne);
 285     }
 286   }
 287 }
 288
 289 /*
 290  * Accessors
 291  */
 292
 293 /// Get the mean radius.
 294 double tcGlobe::meanRadius() const
 295 {
 296   return m_meanRadius;
 297 }