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   state->samplesLon = maxSamples;
 100   state->samplesLat = maxSamples;
 101   state->sampleDelta = (neCorner - swCorner)/(maxSamples-1);
 102   state->swSample = swCorner + state->sampleDelta;
 103 }
 104
 105 /// Get the altitude at a sample in a render state.
 106 double tcSrtmModel::altitudeAt(const RenderState& state, int x, int y, tcGeo* outCoord, bool corrected, bool* accurate)
 107 {
 108   double lon;
 109   if (x <= 0)
 110   {
 111     lon = state.swCorner.lon();
 112   }
 113   else if (x >= state.samplesLon - 1)
 114   {
 115     lon = state.neCorner.lon();
 116   }
 117   else
 118   {
 119     lon = state.swSample.lon() + state.sampleDelta.lon() * (x-1);
 120   }
 121   double lat;
 122   if (y <= 0)
 123   {
 124     lat = state.swCorner.lat();
 125   }
 126   else if (y >= state.samplesLat - 1)
 127   {
 128     lat = state.neCorner.lat();
 129   }
 130   else
 131   {
 132     lat = state.swSample.lat() + state.sampleDelta.lat() * (y-1);
 133   }
 134   return altitudeAt(*outCoord = tcGeo(lon, lat), corrected, accurate);
 135 }
 136
 137 /// Get the altitude at a coordinate.
 138 double tcSrtmModel::altitudeAt(const tcGeo& coord, bool corrected, bool* accurate)
 139 {
 140   // Round down to nearest degree
 141   tcSrtmData* data = 0;
 142   if (0 != m_cache)
 143   {
 144     data = m_cache->loadData(coord);
 145   }
 146   if (0 != data)
 147   {
 148     return data->sampleAltitude(coord, corrected, accurate);
 149   }
 150   return 0.0;
 151 }
 152
 153 /// Get the normal at a coordinate.
 154 maths::Vector<3,float> tcSrtmModel::normalAt(const tcGeo& coord, bool corrected, bool* accurate)
 155 {
 156   if (0 != m_cache)
 157   {
 158     tcSrtmData* data = m_cache->loadData(coord);
 159     if (0 != data)
 160     {
 161       // Sample altitudes a fixed angle in each direction.
 162       const double radiusEarth = 6378.137e3;
 163       tcGeo angRes = data->sampleResolution();
 164       bool xpa, xma, ypa, yma;
 165       double xp = altitudeAt(coord + tcGeo(angRes.lon()*0.5, 0.0), corrected, &xpa);
 166       double xm = altitudeAt(coord - tcGeo(angRes.lon()*0.5, 0.0), corrected, &xma);
 167       double yp = altitudeAt(coord + tcGeo(0.0, angRes.lat()*0.5), corrected, &ypa);
 168       double ym = altitudeAt(coord - tcGeo(0.0, angRes.lat()*0.5), corrected, &yma);
 169       double dx = xp-xm;
 170       double dy = yp-ym;
 171       if (0 != accurate)
 172       {
 173         *accurate = xpa && xma && ypa && yma;
 174       }
 175       maths::Vector<2,float> res(radiusEarth * angRes.lon() * cos(coord.lat()),
 176                                  radiusEarth * angRes.lat());
 177       return maths::Vector<3,float>(
 178         -dx / res[0],
 179         -dy / res[1],
 180         2.0f
 181       ).normalized();
 182     }
 183   }
 184   if (0 != accurate)
 185   {
 186     *accurate = false;
 187   }
 188   return maths::Vector<3,float>(0.0f, 0.0f, 1.0f);
 189 }
 190
 191 /// Use raytracing to find how lit a point is.
 192 float tcSrtmModel::litAt(const tcGeo& coord, const maths::Vector<3,float>& light, float angularRadius, bool corrected, bool* accurate)
 193 {
 194   if (0 != m_cache)
 195   {
 196     tcSrtmData* data = m_cache->loadData(coord);
 197     if (0 != data)
 198     {
 199       // Sample altitudes a fixed angle in each direction.
 200       const float radiusEarth = 6378.137e3;
 201       tcGeo angRes = data->sampleResolution();
 202       float distance = 0.0f;
 203       float step = (angRes.lon() + angRes.lat()) / 4.0;
 204       bool ac = true;
 205       float result = 1.0f;
 206       bool xa;
 207       tcGeo cur = coord;
 208       float alt = altitudeAt(cur, corrected, &xa);
 209       ac = ac && xa;
 210       Q_ASSERT(light[2] > 0.0);
 211       while (alt < 8000.0)
 212       {
 213         // Move towards light source
 214         cur.setLon(cur.lon() + step*light[0] / cos(cur.lat()));
 215         cur.setLat(cur.lat() + step*light[1]);
 216         alt += radiusEarth*step*light[2];
 217         distance += radiusEarth*step;
 218         // Check for collision
 219         float x = altitudeAt(cur, corrected, &xa);
 220         ac = ac && xa;
 221         float flex = distance*angularRadius;
 222         if (x >= alt + flex)
 223         {
 224           result = 0.0f;
 225           break;
 226         }
 227         else if (x >= alt - flex)
 228         {
 229           float portion = (x - alt + flex) / (2.0f*flex);
 230           if (result > portion)
 231           {
 232             result = portion;
 233           }
 234         }
 235         // Can take bigger steps the further we are from the ground
 236         // here we are assuming that the terrain isn't extremely steep
 237         step = 0.5f * (alt - x + flex) / radiusEarth;
 238       }
 239       if (0 != accurate)
 240       {
 241         *accurate = ac;
 242       }
 243       return result;
 244     }
 245   }
 246   if (0 != accurate)
 247   {
 248     *accurate = false;
 249   }
 250   return 1.0f;
 251 }
 252
 253 /// Set the data set to use.
 254 void tcSrtmModel::setDataSet(const QString& name)
 255 {
 256   delete m_cache;
 257   m_cache = new tcSrtmCache(name);
 258   if (0 != m_cache)
 259   {
 260     connect(m_cache, SIGNAL(progressing(const QString&)), this, SIGNAL(progressing(const QString&)));
 261   }
 262   emit dataSetChanged();
 263 }
 264
 265 /// Update any samples affected by the elevation field.
 266 void tcSrtmModel::updateFromElevationData(const tcElevationData* elevation)
 267 {
 268   m_cache->updateFromElevationData(elevation);
 269   emit dataSetChanged();
 270 }