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