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