geo/tcSrtmData.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 tcSrtmData.cpp
  22  * @brief Block of raw SRTM elevation data.
  23  */
  24
  25 #include "tcSrtmData.h"
  26 #include <SplineDefinitions.h>
  27
  28 #include <QFile>
  29 #include <QDataStream>
  30
  31 #include <cmath>
  32
  33 /*
  34  * Constructors + destructor
  35  */
  36
  37 /// Load SRTM data.
  38 tcSrtmData::tcSrtmData(int lon, int lat, QFile& file, QObject* progressObj, const char* progressMember)
  39 : m_lon(lon)
  40 , m_lat(lat)
  41 , m_lines(0)
  42 , m_samples(0)
  43 , m_data(0)
  44 {
  45   connect(this, SIGNAL(progressing(const QString&)), progressObj, progressMember);
  46
  47   qint64 elements = file.size() / sizeof(uint16_t);
  48   double side = sqrt(elements);
  49   m_lines = m_samples = (int)side;
  50   Q_ASSERT(m_lines*m_samples == elements);
  51   m_data = new uint16_t[2*m_lines*m_samples];
  52   Q_ASSERT(0 != m_data);
  53
  54   QDataStream in(&file);
  55   in.setByteOrder(QDataStream::BigEndian);
  56   int lastPercent = -1;
  57   for (int i = 0; i < m_lines; ++i)
  58   {
  59     int percent = 20*i/(m_lines-1);
  60     if (percent != lastPercent)
  61     {
  62       lastPercent = percent;
  63       progressing(tr("Loading SRTM (%1,%2): %3%").arg(lon).arg(lat).arg(percent*5));
  64     }
  65     for (int j = 0; j < m_samples; ++j)
  66     {
  67       Q_ASSERT(!in.atEnd());
  68       uint16_t sample;
  69       in >> sample;
  70       m_data[indexData(0, i, j)] = sample;
  71       m_data[indexData(1, i, j)] = sample;
  72     }
  73   }
  74   /// @todo Handle larger files!
  75   Q_ASSERT(in.atEnd());
  76
  77   progressing(tr("Filling SRTM voids (%1,%2)").arg(lon).arg(lat));
  78   fillVoidsBicubic();
  79
  80   progressing(tr("SRTM (%1,%2): Complete").arg(lon).arg(lat));
  81 }
  82
  83 /// Destructor.
  84 tcSrtmData::~tcSrtmData()
  85 {
  86   delete [] m_data;
  87 }
  88
  89 /*
  90  * Accessors
  91  */
  92
  93 /// Get the longitude in degrees.
  94 int tcSrtmData::lon() const
  95 {
  96   return m_lon;
  97 }
  98
  99 /// Get the latitude in degrees.
 100 int tcSrtmData::lat() const
 101 {
 102   return m_lat;
 103 }
 104
 105 /// Get the number of scan lines.
 106 int tcSrtmData::lines() const
 107 {
 108   return m_lines;
 109 }
 110
 111 /// Get the number of samples in each scan line.
 112 int tcSrtmData::samples() const
 113 {
 114   return m_samples;
 115 }
 116
 117 /// Find the geographical angles between samples.
 118 tcGeo tcSrtmData::sampleResolution() const
 119 {
 120   return tcGeo(M_PI/180 / (m_samples - 1),
 121                M_PI/180 / (m_lines   - 1));
 122 }
 123
 124 /// Get the geographical coordinate of a sample.
 125 tcGeo tcSrtmData::sampleCoordinate(int line, int sample) const
 126 {
 127   return tcGeo(M_PI/180*((double) m_lon    + (double)sample/(m_samples-1)),
 128                M_PI/180*((double)(m_lat+1) - (double)line  /(m_lines  -1)));
 129 }
 130
 131 /// Does data exist for a sample?
 132 bool tcSrtmData::sampleExists(int line, int sample) const
 133 {
 134   Q_ASSERT(line >= 0 && line < m_lines);
 135   Q_ASSERT(sample >= 0 && sample < m_samples);
 136   return 0x0000 == (0x8000 & m_data[indexData(0, line, sample)]);
 137 }
 138
 139 /// Get the altitude at a sample (interpolating gaps).
 140 uint16_t tcSrtmData::sampleAltitude(int line, int sample, bool corrected, bool* exists) const
 141 {
 142   Q_ASSERT(line >= 0 && line < m_lines);
 143   Q_ASSERT(sample >= 0 && sample < m_samples);
 144   if (0 != exists)
 145   {
 146     *exists = (0x0000 == (0x8000 & m_data[indexData(corrected ? 1 : 0, line, sample)]));
 147   }
 148   return 0x7FFF & m_data[indexData(corrected ? 1 : 0, line, sample)];
 149 }
 150
 151 /// Get the altitude at a coordinate.
 152 float tcSrtmData::sampleAltitude(const tcGeo& coord, bool corrected, bool* exists) const
 153 {
 154   double dlon = coord.lon()*180.0/M_PI - m_lon;
 155   double dlat = coord.lat()*180.0/M_PI - m_lat;
 156   Q_ASSERT(dlon >= 0.0 && dlon < 1.0);
 157   Q_ASSERT(dlat >= 0.0 && dlat < 1.0);
 158   int line = (int)floor(dlat*(m_lines-1));
 159   float dy = dlat*(m_lines-1) - line;
 160   line = m_lines - 1 - line;
 161   int sample = (int)floor(dlon*(m_samples-1));
 162   float dx = dlon*(m_samples-1) - sample;
 163   Q_ASSERT(line > 0 && line < m_lines);
 164   Q_ASSERT(sample >= 0 && sample < m_samples-1);
 165   Q_ASSERT(dx >= 0.0f && dx <= 1.0f);
 166   Q_ASSERT(dy >= 0.0f && dy <= 1.0f);
 167   if (0 != exists)
 168   {
 169     *exists = (0x0000 == (0x8000 & m_data[indexData(corrected ? 1 : 0, line, sample)]));
 170   }
 171   uint16_t x00 = 0x7FFF & m_data[indexData(corrected ? 1 : 0, line,   sample  )];
 172   uint16_t x10 = 0x7FFF & m_data[indexData(corrected ? 1 : 0, line,   sample+1)];
 173   uint16_t x01 = 0x7FFF & m_data[indexData(corrected ? 1 : 0, line-1, sample  )];
 174   uint16_t x11 = 0x7FFF & m_data[indexData(corrected ? 1 : 0, line-1, sample+1)];
 175   float mdx = 1.0f-dx;
 176   float mdy = 1.0f-dy;
 177   return ((float)x00*mdx*mdy + (float)x01*mdx*dy + (float)x10*dx*mdy + (float)x11*dx*dy);
 178 }
 179
 180 /*
 181  * Mutators
 182  */
 183
 184 /// Fill voids using bilinear interpolation.
 185 void tcSrtmData::fillVoidsBilinear()
 186 {
 187   for (int line = 0; line < m_lines; ++line)
 188   {
 189     int lastNonVoid = -1;
 190     uint16_t lastNonVoidValue;
 191     for (int sample = 0; sample < m_samples; ++sample)
 192     {
 193       int index = indexData(0, line, sample);
 194       uint16_t value = m_data[index];
 195       bool isVoid = (0 != (value & 0x8000));
 196       if (!isVoid)
 197       {
 198         if (lastNonVoid != -1)
 199         {
 200           int gap = sample - lastNonVoid;
 201           float startAlt = lastNonVoidValue;
 202           float dAlt = (float)(value - lastNonVoidValue) / gap;
 203           // Lovely, between lastNonVoid and sample is a void
 204           for (int i = lastNonVoid+1; i < sample; ++i)
 205           {
 206             int iIndex = indexData(0, line, i);
 207             uint16_t alt = startAlt + dAlt * (i - lastNonVoid);
 208             alt |= 0x8000;
 209             m_data[iIndex] = alt;
 210             // Keep track of the gap
 211             m_data[indexData(1, line, i)] = gap;
 212           }
 213         }
 214         lastNonVoid = sample;
 215         lastNonVoidValue = value;
 216       }
 217       else
 218       {
 219         // Clear any previous data
 220         m_data[index] = 0x8000;
 221       }
 222     }
 223   }
 224   for (int sample = 0; sample < m_samples; ++sample)
 225   {
 226     int lastNonVoid = -1;
 227     uint16_t lastNonVoidValue;
 228     for (int line = 0; line < m_lines; ++line)
 229     {
 230       int index = indexData(0, line, sample);
 231       uint16_t value = m_data[index];
 232       bool isVoid = (0 != (value & 0x8000));
 233       if (!isVoid)
 234       {
 235         if (lastNonVoid != -1)
 236         {
 237           int gap = line - lastNonVoid;
 238           float startAlt = lastNonVoidValue;
 239           float dAlt = (float)(value - lastNonVoidValue) / gap;
 240           // Lovely, between lastNonVoid and line is a void
 241           for (int i = lastNonVoid+1; i < line; ++i)
 242           {
 243             // Average with previous value if non zero, otherwise overwrite
 244             int iIndex = indexData(0, i, sample);
 245             int iIndexC = indexData(1, i, sample);
 246             uint16_t alt = startAlt + dAlt * (i - lastNonVoid);
 247             uint16_t prevAlt = 0x7fff & m_data[iIndex];
 248             if (prevAlt != 0)
 249             {
 250               int otherGap = m_data[iIndexC];
 251               // Prefer the value of the smaller gap
 252               alt = (float)((int)alt * otherGap + (int)prevAlt * gap) / (int)(gap+otherGap);
 253             }
 254             alt |= 0x8000;
 255             m_data[iIndex] = alt;
 256             m_data[iIndexC] = 0x8000; //alt;
 257           }
 258         }
 259         else
 260         {
 261           for (int i = lastNonVoid+1; i < line; ++i)
 262           {
 263             int iIndex = indexData(0, i, sample);
 264             int iIndexC = indexData(1, i, sample);
 265             // The corrected value may contain a gap, which isn't needed anymore
 266             m_data[iIndexC] = 0x8000; //m_data[iIndex];
 267           }
 268         }
 269         lastNonVoid = line;
 270         lastNonVoidValue = value;
 271       }
 272     }
 273   }
 274 }
 275
 276 /// Fill voids using bicubic interpolation.
 277 void tcSrtmData::fillVoidsBicubic()
 278 {
 279   for (int line = 0; line < m_lines; ++line)
 280   {
 281     int lastNonVoid = -1;
 282     uint16_t lastNonVoidValue;
 283     for (int sample = 0; sample < m_samples; ++sample)
 284     {
 285       int index = indexData(0, line, sample);
 286       uint16_t value = m_data[index];
 287       bool isVoid = (0 != (value & 0x8000));
 288       if (!isVoid)
 289       {
 290         if (lastNonVoid != -1)
 291         {
 292           int gap = sample - lastNonVoid;
 293           // Obtain some control points for the splines
 294           // use the two samples either side of the gap as control points 1 and 2
 295           // use the ones before and after these samples, extended to the size of
 296           // the gap as control points 0 and 3
 297           float startAlt = lastNonVoidValue;
 298           float beforeStartAlt = startAlt;
 299           if (lastNonVoid > 0)
 300           {
 301             uint16_t beforeStart = m_data[indexData(0, line, lastNonVoid-1)];
 302             if (0 == (beforeStart & 0x8000))
 303             {
 304               beforeStartAlt = startAlt + ((float)beforeStart - startAlt) * gap;
 305             }
 306           }
 307           float endAlt = value;
 308           float afterEndAlt = endAlt;
 309           if (sample < m_samples-1)
 310           {
 311             uint16_t afterEnd = m_data[indexData(0, line, sample+1)];
 312             if (0 == (afterEnd & 0x8000))
 313             {
 314               afterEndAlt = endAlt + ((float)afterEnd - endAlt) * gap;
 315             }
 316           }
 317           // Lovely, between lastNonVoid and sample is a void
 318           for (int i = lastNonVoid+1; i < sample; ++i)
 319           {
 320             int iIndex = indexData(0, line, i);
 321             float u = (float)(i-lastNonVoid) / gap;
 322             float u2 = u*u;
 323             float u3 = u2*u;
 324             uint16_t alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 325             alt |= 0x8000;
 326             m_data[iIndex] = alt;
 327             // Keep track of the gap
 328             m_data[indexData(1, line, i)] = gap;
 329           }
 330         }
 331         lastNonVoid = sample;
 332         lastNonVoidValue = value;
 333       }
 334       else
 335       {
 336         // Clear any previous data
 337         m_data[index] = 0x8000;
 338       }
 339     }
 340   }
 341   for (int sample = 0; sample < m_samples; ++sample)
 342   {
 343     int lastNonVoid = -1;
 344     uint16_t lastNonVoidValue;
 345     for (int line = 0; line < m_lines; ++line)
 346     {
 347       int index = indexData(0, line, sample);
 348       uint16_t value = m_data[index];
 349       bool isVoid = (0 != (value & 0x8000));
 350       if (!isVoid)
 351       {
 352         if (lastNonVoid != -1)
 353         {
 354           int gap = line - lastNonVoid;
 355           // Obtain some control points for the splines
 356           // use the two samples either side of the gap as control points 1 and 2
 357           // use the ones before and after these samples, extended to the size of
 358           // the gap as control points 0 and 3
 359           float startAlt = lastNonVoidValue;
 360           float beforeStartAlt = startAlt;
 361           if (lastNonVoid > 0)
 362           {
 363             uint16_t beforeStart = m_data[indexData(0, lastNonVoid-1, sample)];
 364             if (0 == (beforeStart & 0x8000))
 365             {
 366               beforeStartAlt = startAlt + ((float)beforeStart - startAlt) * gap;
 367             }
 368           }
 369           float endAlt = value;
 370           float afterEndAlt = endAlt;
 371           if (line < m_lines-1)
 372           {
 373             uint16_t afterEnd = m_data[indexData(0, line+1, sample)];
 374             if (0 == (afterEnd & 0x8000))
 375             {
 376               afterEndAlt = endAlt + ((float)afterEnd - endAlt) * gap;
 377             }
 378           }
 379           // Lovely, between lastNonVoid and line is a void
 380           for (int i = lastNonVoid+1; i < line; ++i)
 381           {
 382             // Average with previous value if non zero, otherwise overwrite
 383             int iIndex = indexData(0, i, sample);
 384             int iIndexC = indexData(1, i, sample);
 385             float u = (float)(i-lastNonVoid) / gap;
 386             float u2 = u*u;
 387             float u3 = u2*u;
 388             uint16_t alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 389             uint16_t prevAlt = 0x7fff & m_data[iIndex];
 390             if (prevAlt != 0)
 391             {
 392               int otherGap = m_data[iIndexC];
 393               // Prefer the value of the smaller gap
 394               alt = (float)((int)alt * otherGap + (int)prevAlt * gap) / (int)(gap+otherGap);
 395             }
 396             alt |= 0x8000;
 397             m_data[iIndex] = alt;
 398             m_data[iIndexC] = 0x8000; //alt;
 399           }
 400         }
 401         else
 402         {
 403           for (int i = lastNonVoid+1; i < line; ++i)
 404           {
 405             int iIndex = indexData(0, i, sample);
 406             int iIndexC = indexData(1, i, sample);
 407             // The corrected value may contain a gap, which isn't needed anymore
 408             m_data[iIndexC] = 0x8000; //m_data[iIndex];
 409           }
 410         }
 411         lastNonVoid = line;
 412         lastNonVoidValue = value;
 413       }
 414     }
 415   }
 416 }
 417
 418 /*
 419  * Private functions
 420  */
 421
 422 /// Index into the data array.
 423 unsigned int tcSrtmData::indexData(int corrected, int line, int sample) const
 424 {
 425   return (corrected*m_lines + line)*m_samples + sample;
 426 }