geo/tcElevationData.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 tcElevationData.cpp
  22  * @brief A block of pixel data representing elevation.
  23  */
  24
  25 #include "tcElevationData.h"
  26
  27 /*
  28  * Constructors + destructor
  29  */
  30
  31 /// Default constructor.
  32 tcElevationData::tcElevationData()
  33 : Parent()
  34 , m_geoToPixels()
  35 {
  36 }
  37
  38 /// Primary constructor.
  39 tcElevationData::tcElevationData(int width, int height)
  40 : Parent(width, height)
  41 , m_geoToPixels()
  42 {
  43 }
  44
  45 /// Destructor.
  46 tcElevationData::~tcElevationData()
  47 {
  48 }
  49
  50 /*
  51  * Elevation interface
  52  */
  53
  54 /// Sample elevation at a coordinate.
  55 float tcElevationData::elevationAt(const tcGeo& coord) const
  56 {
  57   maths::Vector<2,double> pt = m_geoToPixels * coord;
  58   return sampleFloat(pt[0]/(width()-1), pt[1]/(height()-1));
  59 }
  60
  61 /// Get the max elevation over a range of pixels.
  62 float tcElevationData::maxElevation(const tcGeo& swCorner, const tcGeo& neCorner, bool* inRange) const
  63 {
  64   maths::Vector<2,double> corner1 = m_geoToPixels * swCorner;
  65   maths::Vector<2,double> corner2 = m_geoToPixels * neCorner;
  66
  67   maths::Vector<2,double> minCorner(qMin(corner1[0], corner2[0]),
  68                                     qMin(corner1[1], corner2[1]));
  69   maths::Vector<2,double> maxCorner(qMax(corner1[0], corner2[0]),
  70                                     qMax(corner1[1], corner2[1]));
  71
  72   // Go through pixels in this range.
  73   int count = 0;
  74   float max = 0.0f;
  75   int minj = qMax((int)ceil(minCorner[1]), 0);
  76   int maxj = qMin((int)ceil(maxCorner[1]), height());
  77   int mini = qMax((int)ceil(minCorner[0]), 0);
  78   int maxi = qMin((int)ceil(maxCorner[0]), width());
  79   for (int j = minj; j < maxj; ++j)
  80   {
  81     for (int i = mini; i < maxi; ++i)
  82     {
  83       int index = j*width() + i;
  84       ++count;
  85       max = qMax(max, buffer()[index]);
  86     }
  87   }
  88
  89   if (0 != inRange)
  90   {
  91     *inRange = true;
  92   }
  93
  94   if (count > 0)
  95   {
  96     // Just average the samples
  97     return max;
  98   }
  99   else
 100   {
 101     // No samples? resort to interpolation
 102     maths::Vector<2,float> mid(corner1 + corner2);
 103     mid /= 2;
 104     mid[0] /= width() -1;
 105     mid[1] /= height()-1;
 106     if (mid[0] >= 0.0f && mid[0] <= 1.0f &&
 107         mid[1] >= 0.0f && mid[1] <= 1.0f)
 108     {
 109       return Parent::sampleFloat(mid[0], mid[1]);
 110     }
 111     else
 112     {
 113       if (0 != inRange)
 114       {
 115         *inRange = false;
 116       }
 117       return 0.0f;
 118     }
 119   }
 120 }
 121
 122 /// Get the mean elevation over a range of pixels.
 123 float tcElevationData::meanElevation(const tcGeo& swCorner, const tcGeo& neCorner, bool* inRange) const
 124 {
 125   maths::Vector<2,double> corner1 = m_geoToPixels * swCorner;
 126   maths::Vector<2,double> corner2 = m_geoToPixels * neCorner;
 127
 128   maths::Vector<2,double> minCorner(qMin(corner1[0], corner2[0]),
 129                                     qMin(corner1[1], corner2[1]));
 130   maths::Vector<2,double> maxCorner(qMax(corner1[0], corner2[0]),
 131                                     qMax(corner1[1], corner2[1]));
 132
 133   // Go through pixels in this range.
 134   int count = 0;
 135   float total = 0.0f;
 136   int minj = qMax((int)ceil(minCorner[1]), 0);
 137   int maxj = qMin((int)ceil(maxCorner[1]), height());
 138   int mini = qMax((int)ceil(minCorner[0]), 0);
 139   int maxi = qMin((int)ceil(maxCorner[0]), width());
 140   for (int j = minj; j < maxj; ++j)
 141   {
 142     for (int i = mini; i < maxi; ++i)
 143     {
 144       int index = j*width() + i;
 145       ++count;
 146       total += buffer()[index];
 147     }
 148   }
 149
 150   if (0 != inRange)
 151   {
 152     *inRange = true;
 153   }
 154
 155   if (count > 0)
 156   {
 157     // Just average the samples
 158     return total / count;
 159   }
 160   else
 161   {
 162     // No samples? resort to interpolation
 163     maths::Vector<2,float> mid(corner1 + corner2);
 164     mid /= 2;
 165     mid[0] /= width() -1;
 166     mid[1] /= height()-1;
 167     if (mid[0] >= 0.0f && mid[0] <= 1.0f &&
 168         mid[1] >= 0.0f && mid[1] <= 1.0f)
 169     {
 170       return Parent::sampleFloat(mid[0], mid[1]);
 171     }
 172     else
 173     {
 174       if (0 != inRange)
 175       {
 176         *inRange = false;
 177       }
 178       return 0.0f;
 179     }
 180   }
 181 }
 182
 183 /// Set the geographical to pixels transformation.
 184 void tcElevationData::setGeoToPixels(const tcAffineTransform2<double>& geoToPixels)
 185 {
 186   m_geoToPixels = geoToPixels;
 187 }
 188
 189 /*
 190  * Accessors
 191  */
 192
 193 /// Get the coordinate of the most south westerly corner.
 194 tcGeo tcElevationData::swCorner() const
 195 {
 196   tcAffineTransform2<double> pixelsToGeo = m_geoToPixels.inverse();
 197   return (tcGeo)(pixelsToGeo * maths::Vector<2,double>(0.0, height()-1));
 198 }
 199
 200 /// Get the coordinate of the most north easterly corner.
 201 tcGeo tcElevationData::neCorner() const
 202 {
 203   tcAffineTransform2<double> pixelsToGeo = m_geoToPixels.inverse();
 204   return (tcGeo)(pixelsToGeo * maths::Vector<2,double>(width()-1, 0.0));
 205 }