geo/tcSrtmModel.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 tcSrtmModel.cpp
  22  * @brief SRTM elevation model.
  23  */
  24
  25 #include "tcSrtmModel.h"
  26 #include "tcSrtmCache.h"
  27 #include "tcSrtmData.h"
  28
  29 #include <inttypes.h>
  30
  31 /*
  32  * Constructors + destructor
  33  */
  34
  35 /// Default constructor.
  36 tcSrtmModel::tcSrtmModel()
  37 : m_cache(0)
  38 {
  39 }
  40
  41 /// Destructor.
  42 tcSrtmModel::~tcSrtmModel()
  43 {
  44   delete m_cache;
  45 }
  46
  47 /*
  48  * Main interface
  49  */
  50
  51 /// Get information to allow alignment with actual samples.
  52 void tcSrtmModel::sampleAlign(const tcGeo& swCorner, const tcGeo& neCorner, RenderState* state, int maxSamples)
  53 {
  54   state->swCorner = swCorner;
  55   state->neCorner = neCorner;
  56   state->moreAvailableLon = false;
  57   state->moreAvailableLat = false;
  58
  59   tcSrtmData* data = 0;
  60   if (0 != m_cache)
  61   {
  62     data = m_cache->findData(swCorner, neCorner);
  63     if (data)
  64     {
  65       tcGeo resolution = data->sampleResolution();
  66       maths::Vector<2,double> swSampleSpace = swCorner / resolution;
  67       maths::Vector<2,double> neSampleSpace = neCorner / resolution;
  68       maths::Vector<2,double> first(floor(swSampleSpace[0] + 1.0),
  69                                     floor(swSampleSpace[1] + 1.0));
  70       maths::Vector<2,double> last( ceil( neSampleSpace[0] - 1.0),
  71                                     ceil( neSampleSpace[1] - 1.0));
  72       state->samplesLon = 3 + (int)(last[0] - first[0]);
  73       state->samplesLat = 3 + (int)(last[1] - first[1]);
  74       state->swSample = swCorner + resolution*(first - swSampleSpace);
  75       state->sampleDelta = resolution;
  76
  77       // Sample evenly
  78       if (maxSamples >= 2)
  79       {
  80         if (state->samplesLon > maxSamples)
  81         {
  82           state->sampleDelta.setLon((neCorner.lon()-swCorner.lon()) / (maxSamples-1));
  83           state->swSample.setLon(swCorner.lon() + state->sampleDelta.lon());
  84           state->samplesLon = maxSamples;
  85           state->moreAvailableLon = true;
  86         }
  87         if (state->samplesLat > maxSamples)
  88         {
  89           state->sampleDelta.setLat((neCorner.lat()-swCorner.lat()) / (maxSamples-1));
  90           state->swSample.setLat(swCorner.lat() + state->sampleDelta.lat());
  91           state->samplesLat = maxSamples;
  92           state->moreAvailableLat = true;
  93         }
  94       }
  95       return;
  96     }
  97   }
  98
  99   int samplesLon = (int)((neCorner.lon() - swCorner.lon())/(M_PI/180.0/60.0));
 100   int samplesLat = (int)((neCorner.lat() - swCorner.lat())/(M_PI/180.0/60.0));
 101   if (samplesLon > maxSamples)
 102   {
 103     state->moreAvailableLon = true;
 104   }
 105   if (samplesLat > maxSamples)
 106   {
 107     state->moreAvailableLat = true;
 108   }
 109
 110
 111   state->samplesLon = maxSamples;
 112   state->samplesLat = maxSamples;
 113   state->sampleDelta = (neCorner - swCorner)/(maxSamples-1);
 114   state->swSample = swCorner + state->sampleDelta;
 115 }
 116
 117 /// Get the altitude at a sample in a render state.
 118 double tcSrtmModel::altitudeAt(const RenderState& state, int x, int y, tcGeo* outCoord, bool corrected, bool* accurate)
 119 {
 120   double lon;
 121   if (x <= 0)
 122   {
 123     lon = state.swCorner.lon();
 124   }
 125   else if (x >= state.samplesLon - 1)
 126   {
 127     lon = state.neCorner.lon();
 128   }
 129   else
 130   {
 131     lon = state.swSample.lon() + state.sampleDelta.lon() * (x-1);
 132   }
 133   double lat;
 134   if (y <= 0)
 135   {
 136     lat = state.swCorner.lat();
 137   }
 138   else if (y >= state.samplesLat - 1)
 139   {
 140     lat = state.neCorner.lat();
 141   }
 142   else
 143   {
 144     lat = state.swSample.lat() + state.sampleDelta.lat() * (y-1);
 145   }
 146   return altitudeAt(*outCoord = tcGeo(lon, lat), corrected, accurate);
 147 }
 148
 149 /// Get the altitude at a coordinate.
 150 double tcSrtmModel::altitudeAt(const tcGeo& coord, bool corrected, bool* accurate)
 151 {
 152   // Round down to nearest degree
 153   tcSrtmData* data = 0;
 154   if (0 != m_cache)
 155   {
 156     data = m_cache->loadData(coord);
 157   }
 158   if (0 != data)
 159   {
 160     return data->sampleAltitude(coord, corrected, accurate);
 161   }
 162   return 0.0;
 163 }
 164
 165 /// Get the normal at a coordinate.
 166 maths::Vector<3,float> tcSrtmModel::normalAt(const tcGeo& coord, bool corrected, bool* accurate)
 167 {
 168   if (0 != m_cache)
 169   {
 170     tcSrtmData* data = m_cache->loadData(coord);
 171     if (0 != data)
 172     {
 173       // Sample altitudes a fixed angle in each direction.
 174       const double radiusEarth = 6378.137e3;
 175       tcGeo angRes = data->sampleResolution();
 176       bool alta, xpa, xma, ypa, yma;
 177       double alt = altitudeAt(coord, corrected, &alta);
 178       double xp = altitudeAt(coord + tcGeo(angRes.lon()*0.5, 0.0), corrected, &xpa);
 179       double xm = altitudeAt(coord - tcGeo(angRes.lon()*0.5, 0.0), corrected, &xma);
 180       double yp = altitudeAt(coord + tcGeo(0.0, angRes.lat()*0.5), corrected, &ypa);
 181       double ym = altitudeAt(coord - tcGeo(0.0, angRes.lat()*0.5), corrected, &yma);
 182       double dx = xp-xm;
 183       double dy = yp-ym;
 184       if (0 != accurate)
 185       {
 186         *accurate = xpa && xma && ypa && yma;
 187       }
 188       maths::Vector<2,float> res((radiusEarth+alt) * angRes.lon() * cos(coord.lat()),
 189                                  (radiusEarth+alt) * angRes.lat());
 190       // (1.0f, 0.0f, dx/res[0]) x (0.0f, -1.0f, dy/res[1])
 191       maths::Vector<3,float> norm(-dx/res[0], -dy/res[1], 1.0f);
 192       norm.normalize();
 193       return norm;
 194     }
 195   }
 196   if (0 != accurate)
 197   {
 198     *accurate = false;
 199   }
 200   return maths::Vector<3,float>(0.0f, 0.0f, 1.0f);
 201 }
 202
 203 /// Use raytracing to find how lit a point is.
 204 float tcSrtmModel::litAt(const tcGeo& coord, const maths::Vector<3,float>& light, float angularRadius, bool corrected, bool* accurate)
 205 {
 206   if (0 != m_cache)
 207   {
 208     tcSrtmData* data = m_cache->loadData(coord);
 209     if (0 != data)
 210     {
 211       // Sample altitudes a fixed angle in each direction.
 212       const float radiusEarth = 6378.137e3;
 213       tcGeo angRes = data->sampleResolution();
 214       float distance = 0.0f;
 215       float step = (angRes.lon() + angRes.lat()) / 4.0;
 216       bool ac = true;
 217       float result = 1.0f;
 218       bool xa;
 219       tcGeo cur = coord;
 220       float alt = altitudeAt(cur, corrected, &xa);
 221       ac = ac && xa;
 222       Q_ASSERT(light[2] > 0.0);
 223       while (alt < 8000.0)
 224       {
 225         // Move towards light source
 226         cur.setLon(cur.lon() + step*light[0] / cos(cur.lat()));
 227         cur.setLat(cur.lat() + step*light[1]);
 228         alt += radiusEarth*step*light[2];
 229         distance += radiusEarth*step;
 230         // Check for collision
 231         float x = altitudeAt(cur, corrected, &xa);
 232         ac = ac && xa;
 233         float flex = distance*angularRadius;
 234         if (x >= alt + flex)
 235         {
 236           result = 0.0f;
 237           break;
 238         }
 239         else if (x >= alt - flex)
 240         {
 241           float portion = (x - alt + flex) / (2.0f*flex);
 242           if (result > portion)
 243           {
 244             result = portion;
 245           }
 246         }
 247         // Can take bigger steps the further we are from the ground
 248         // here we are assuming that the terrain isn't extremely steep
 249         step = 0.5f * (alt - x + flex) / radiusEarth;
 250       }
 251       if (0 != accurate)
 252       {
 253         *accurate = ac;
 254       }
 255       return result;
 256     }
 257   }
 258   if (0 != accurate)
 259   {
 260     *accurate = false;
 261   }
 262   return 1.0f;
 263 }
 264
 265 /// Set the data set to use.
 266 void tcSrtmModel::setDataSet(const QString& name)
 267 {
 268   delete m_cache;
 269   if (!name.isEmpty())
 270   {
 271     m_cache = new tcSrtmCache(name);
 272     if (0 != m_cache)
 273     {
 274       connect(m_cache, SIGNAL(progressing(const QString&)), this, SIGNAL(progressing(const QString&)));
 275     }
 276   }
 277   else
 278   {
 279     m_cache = 0;
 280   }
 281   emit dataSetChanged();
 282 }
 283
 284 /// Update any samples affected by the elevation field.
 285 void tcSrtmModel::updateFromElevationData(const tcElevationData* elevation)
 286 {
 287   if (0 != m_cache)
 288   {
 289     m_cache->updateFromElevationData(elevation);
 290     emit dataSetChanged();
 291   }
 292 }