/**
  * @RDE_HEADER@
  * adapting GeISiR source to SituationMap
  * $Id$
  * $HeadURL$
  * $Author$
  */


#include "elevation_map.h"

#include <QtConcurrent>
#include <QPolygonF>
#include <MathUtils/math_const.h>
#include <QElapsedTimer>

#include <gdal.h>

#include <iostream>
#include <limits>
#include <memory>
#include <math.h>

namespace SituationMap
{

class BoundingRect
{
  double m_top = 0.0;
  double m_bottom = 0.0;
  double m_left = 0.0;
  double m_right = 0.0;

  void swapBounds()
  {
    if (m_top < m_bottom)
      std::swap(m_bottom, m_top);
    if (m_right < m_left)
      std::swap(m_left, m_right);
  } // void swapBounds()

public:

  BoundingRect() = default;

  BoundingRect(double left, double bottom, double right, double top)
    : m_top(top), m_bottom(bottom), m_left(left), m_right(right)
  {
    swapBounds();
  } // BoundingRect( double left, double bottom, double right, double top)

  BoundingRect(const BoundingRect &) = default;
  BoundingRect &operator =(const BoundingRect &) = default;

  ~BoundingRect() = default;

  inline void operator()(double left, double bottom, double right, double top)
  {
    m_top    = top;
    m_bottom = bottom;
    m_left   = left;
    m_right  = right;

    swapBounds();
  } // inline void operator() (double left, double bottom, double right, double top)

  double top()    const { return m_top;    }
  double bottom() const { return m_bottom; }
  double left()   const { return m_left;   }
  double right()  const { return m_right;  }

  bool contains(const QPointF &p) const { return ((p.x() >= m_left) && (p.x() <= m_right) && (p.y() >= m_bottom) && (p.y() <= m_top)); }

  QRectF toQRectF() const
  {
    return QRectF(QPointF(m_left, m_top), QPointF(m_right, m_bottom));
  } // toQRectF
}; // BoundingRect



struct ElevationMap::Internal
{
  BoundingRect     _location;         // Lage der Karte in Lat/Lon

  struct PointDouble
  {
    void setX(double x) { _x = x;}
    void setY(double y) { _y = y;}
    double x() { return _x;}
    double y() { return _y;}
    double _x{0.0};
    double _y{0.0};
  };
  PointDouble      _elexelSize;      // LatLon Size of 1 elevation sample, redundat

  int              _totalColumns = 0;
  int              _totalRows = 0;

  double           _toMeterFactor = 1.0;
  std::vector<int> _elevationField;

  template<class BufferType> bool copyBufferData(GDALRasterBandH hBand, int xBlkSize, int yBlkSize)
  {
    // NLDEVCTRL-804: Konsistenz d. Z-Werte zwischen 2D/3D (also ElevationMap vs. ESRI::BaseSurface)
    // This ElevationMap was using an approach similar to following.
    // Now it´s more or less identical and the difference using 'Altmark_dt1.tif'
    // are gone. AHTC is also fine :-)
    /*
    CPLErr GetHistogram( GDALRasterBand *poBand, GUIntBig *panHistogram )
    {
    memset( panHistogram, 0, sizeof(GUIntBig) * 256 );

    CPLAssert( poBand->GetRasterDataType() == GDT_Byte );

    int nXBlockSize, nYBlockSize;

    poBand->GetBlockSize( &nXBlockSize, &nYBlockSize );
    int nXBlocks = (poBand->GetXSize() + nXBlockSize - 1) / nXBlockSize;
    int nYBlocks = (poBand->GetYSize() + nYBlockSize - 1) / nYBlockSize;

    GByte *pabyData = (GByte *) CPLMalloc(nXBlockSize * nYBlockSize);

    for( int iYBlock = 0; iYBlock < nYBlocks; iYBlock++ )
    {
        for( int iXBlock = 0; iXBlock < nXBlocks; iXBlock++ )
        {
            int        nXValid, nYValid;

            poBand->ReadBlock( iXBlock, iYBlock, pabyData );

            // Compute the portion of the block that is valid
            // for partial edge blocks.
            poBand->GetActualBlockSize(iXBlock, iYBlock, &nXValid, &nYValid)

            // Collect the histogram counts.
            for( int iY = 0; iY < nYValid; iY++ )
            {
                for( int iX = 0; iX < nXValid; iX++ )
                {
                    panHistogram[pabyData[iX + iY * nXBlockSize]] += 1;
                }
            }
        }
    }
    }
    Source https://gdal.org/api/gdalrasterband_cpp.html
    */

    const size_t  bufSize = xBlkSize * yBlkSize;
    std::unique_ptr<BufferType[]> bufferPtr(new BufferType[bufSize]);

    int nXBlocks = (GDALGetRasterBandXSize(hBand) + xBlkSize - 1) / xBlkSize;
    int nYBlocks = (GDALGetRasterBandYSize(hBand) + yBlkSize - 1) / yBlkSize;

    for (int iYBlock = 0; iYBlock < nYBlocks; ++iYBlock)
    {
      for (int iXBlock = 0; iXBlock < nXBlocks; ++iXBlock)
      {

        CPLErr err = GDALReadBlock( hBand, iXBlock, iYBlock, bufferPtr.get() );  // Gdal liest spaltenweise !! x <-> y
        if (err != CE_None)
        {
          return false;
        }

        int nXValid, nYValid;

        // Unfortunately GDALGetActualBlockSize is not available at Suse 42.3 (gdal 2.1.3)
        // so check bounds using _totalColumns, _totalRows:
#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(2, 2, 0)
#  define HAS_ACTUAL_BLOCKSIZE
#endif

#ifdef HAS_ACTUAL_BLOCKSIZE
        GDALGetActualBlockSize(hBand, iXBlock, iYBlock, &nXValid, &nYValid);
#else
        nXValid = xBlkSize;
        nYValid = yBlkSize;
#endif // HAS_ACTUAL_BLOCKSIZE


        const int curYOffset = (iYBlock * yBlkSize);
        const int curXOffset = (iXBlock * xBlkSize);

        for (int iY = 0; iY < nYValid; ++iY)
        {
          const int row    = iY + curYOffset; // Line

#ifndef HAS_ACTUAL_BLOCKSIZE
          if (row >= _totalRows)
            continue;
#endif // HAS_ACTUAL_BLOCKSIZE

          for (int iX = 0; iX < nXValid; ++iX)
          {
            // identical to gdallocationinfo
            const int column = iX + curXOffset; // Pixel

#ifndef HAS_ACTUAL_BLOCKSIZE
            if (column >= _totalColumns)
              continue;
#endif // HAS_ACTUAL_BLOCKSIZE

            int bi = iX + (iY * xBlkSize);
            const int z = static_cast<int>( static_cast<double>(bufferPtr[bi]) / _toMeterFactor );

            _elevationField[column + row * _totalColumns] = z;
          } // for
        } // for

      } // for
    } // for

    return true;

  } // template<class BufferType> void copyBufferData(BufferType * bufferPtr, int xBlkSize, int yBlkSize)

}; // struct ElevationMap::Internal


ElevationMap::ElevationMap(const QFileInfo &fileInfo)
  : m_internal(new Internal())
  , m_aproximation(true)
  , m_valid(false)
  , m_logger(getLogger("ElevationMap"))
{
  GDALAllRegister();

  m_fileName = fileInfo.absoluteFilePath();

  INFO_LOG(m_logger) << " ctor ElevationMap: map of file '" << m_fileName << "'"
                       << ", this=" << this;

  if (!fileInfo.exists())
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " elevation data file \'" << m_fileName << "\' doesn't exists.";
    return;
  }

  if (!fileInfo.isReadable())
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " elevation data file \'" << m_fileName << "\' isn't readable.";
    return;
  }

  auto cb = [this]()
  {
    return readFromFile();
  };
  auto cbFinished = [this]()
  {
    NOTICE_LOG(m_logger) << " Finished reading elevation map of file '" << m_fileName << "'"
                         << ", this=" << this;
  };
  QObject::connect(&m_gdalReadWatcher, &QFutureWatcher<bool>::finished,
                   &m_gdalReadWatcher, cbFinished);
  m_gdalReadWatcher.setFuture(QtConcurrent::run(cb));
  // sync init, takes ~800ms during sitmap loading 72MB file
  //  QElapsedTimer t;
  //  t.start();
  //  m_gdalReadWatcher.waitForFinished();
  //  NOTICE_LOG(m_logger) << "Time loading elevation map: " << t.elapsed();
} // ElevationMap::ElevationMap( const QFile &file )


ElevationMap::~ElevationMap()
{
  TRACE_LOG(m_logger) << LOGFUNCTION << " START";
  m_gdalReadWatcher.cancel();
  m_gdalReadWatcher.waitForFinished();

  delete m_internal;
  m_internal = nullptr;
  TRACE_LOG(m_logger) << LOGFUNCTION << " END";
} // ElevationMap::~ElevationMap

bool ElevationMap::waitForInitialization()
{
  m_gdalReadWatcher.waitForFinished();
  return m_valid;
}

bool ElevationMap::isValid() const
{
  // N. Thoden, 27.04.2020
  // isFinished seems to be unreliable:
  // if mainThread request ::isValid() diretly after ::waitForInitialization()
  // you get false (by 'if (!m_gdalReadWatcher.isFinished())'
  // but m_valid is true. So skip the additional check.
//  if (!m_gdalReadWatcher.isFinished())
//    return false;

  return m_valid;
}

QString ElevationMap::filename() const
{
  return m_fileName;
}

bool ElevationMap::readFromFile()
{
  m_valid = false;

  GDALDatasetH hDataset = GDALOpen(m_fileName.toLocal8Bit().constData(), GA_ReadOnly);

  if (nullptr == hDataset)
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " GDALOpen failed, may be calling of GDALAllRegister is missing";
    return false;
  }

  // Detect basic parameters
  //
  GDALDriverH hDriver             = GDALGetDatasetDriver(hDataset);
  int         nRasterXSize        = GDALGetRasterXSize(hDataset);
  int         nRasterYSize        = GDALGetRasterYSize(hDataset);
  double      adfGeoTransform[6];
  CPLErr      err                 = GDALGetGeoTransform(hDataset, adfGeoTransform); // XY Spatial Referenz WGS84 (4326) expected!!!

  if (err != CE_None)
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to get geo raster transform.";
    GDALClose(hDataset);
    return false;
  }


  if (nullptr == hDriver) { GDALClose(hDataset); return false; }
  if (nRasterXSize <= 0)  { GDALClose(hDataset); return false; }
  if (nRasterYSize <= 0)  { GDALClose(hDataset); return false; }

  double pixWidth  = adfGeoTransform[1];
  double pixHeight = adfGeoTransform[5];

  auto toRAD = [](double x) -> double { return x * MathConst::pi / 180.0; };


  m_internal->_location(toRAD(adfGeoTransform[0])       // adfGeoTransform is in degrees
                        , toRAD(adfGeoTransform[3])
                        , toRAD(adfGeoTransform[0] + double(nRasterXSize) * pixWidth + double(nRasterYSize) * adfGeoTransform[2])
                        , toRAD(adfGeoTransform[3] + double(nRasterYSize) * pixHeight + double(nRasterXSize) * adfGeoTransform[4]));

  m_internal->_elexelSize.setX(fabs(toRAD(pixWidth)));
  m_internal->_elexelSize.setY(fabs(toRAD(pixHeight)));

  GDALRasterBandH hBand = GDALGetRasterBand(hDataset, 1);

  if (nullptr == hBand)
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to read raster band.";
    GDALClose(hDataset);
    return false;
  }
 
  int xBlkSize;
  int yBlkSize;
  GDALGetBlockSize(hBand, &xBlkSize, &yBlkSize);

  if (xBlkSize <0 || yBlkSize < 0)
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " no buffer, next operations will possible fail....";
    return false;
  }

  m_internal->_totalColumns   = nRasterXSize;
  m_internal->_totalRows      = nRasterYSize;
  m_internal->_elevationField.resize(nRasterXSize * nRasterYSize);
 
  const GDALDataType type = GDALGetRasterDataType(hBand);
  switch (type)
  {
  case GDT_UInt16:
    m_internal->_toMeterFactor = 1.0;
    m_valid = m_internal->copyBufferData<unsigned short int> ( hBand, xBlkSize, yBlkSize );
    break;
  case GDT_Int16:
    m_internal->_toMeterFactor = 1.0;
    m_valid = m_internal->copyBufferData<short int>          ( hBand, xBlkSize, yBlkSize );
    break;
  case GDT_Float32:
    m_internal->_toMeterFactor = 0.01;
    m_valid = m_internal->copyBufferData<float>              ( hBand, xBlkSize, yBlkSize );
    break;
  case GDT_Float64:
    m_internal->_toMeterFactor = 0.01;
    m_valid = m_internal->copyBufferData<double>             ( hBand, xBlkSize, yBlkSize );
    break;
  default:
  {
    ERROR_LOG(m_logger) << LOGFUNCTION << " data type isn't supported type or a floating point type"
                        << " (type=" << type << ")";
    GDALClose(hDataset);
    return false;
  }
  } // switch

  if (!m_valid)
  {
    CPLErrorNum errNum = CPLGetLastErrorNo();
    ERROR_LOG(m_logger) << LOGFUNCTION << " GDAL failed to read data block (error no " << int(errNum) << ")";
  }
 
  GDALClose(hDataset);
  
  INFO_LOG(m_logger) << LOGFUNCTION << " loaded, filename =" << m_fileName << "'"
                       << ", m_valid=" << m_valid;
  return m_valid;

} // bool ElevationMap::readFromFile(const QString & fileName)




QRectF ElevationMap::extent() const
{
  if (!isValid())
    return QRectF();

  return m_internal->_location.toQRectF();
}

QPair<double, double> ElevationMap::elexelSize() const
{
  return {m_internal->_elexelSize._x, m_internal->_elexelSize._y};
}


int ElevationMap::getElevation(int column, int row) const
{
  row = m_internal->_totalRows - row;
  const size_t index = column + row * m_internal->_totalColumns;

  if (m_internal->_elevationField.empty() || m_internal->_elevationField.size() < index)
    return std::numeric_limits<int>::min();

  return m_internal->_elevationField.at(index);
} // int ElevationMap::getElevation( int column, int row)


double ElevationMap::elevation(double lon, double lat) const
{
  if (!isValid())
    return 0;

  return getElevation(lon, lat) * m_internal->_toMeterFactor;
} // double ElevationMap::elevation( double lon, double lat)


double ElevationMap:: getElevation(double lon, double lat) const
{
  if (!m_internal->_location.contains(QPointF(lon, lat)))
  {
    return std::numeric_limits<int>::min();  // ausserhalb der Höhenkarte !!
  }
  double x = (lon - m_internal->_location.left()) / m_internal->_elexelSize.x();
  int column = floor(x);

  double y = (lat - m_internal->_location.bottom()) / m_internal->_elexelSize.y();
  int row = floor(y);

  int elevation = getElevation(column, row);
  if (elevation == std::numeric_limits<int>::min())
    return elevation;

  if (!m_aproximation)
    return elevation;

  /*
  ** Approximation: e0, e1, e2, e3 Hoehenwerte der 4 benachbarten Stützpunkte
  **
  **        + e3        + e2
  **     y
  **     ^        + er                er = e0 + x*(e1-e0) + y*(e3-e0) + x*y(e2-e3+e0-e1)
  **     |
  **     |
  **     |
  **     |  + e0        + e1
  **       -----> x
  **
  ** Besonderheiten: die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Rasters,
  ** so dass je nach Groesse von x, y andere Nachbarpunkte zu berücksichtigen sind, und x, sowie y entsprechend
  ** angepasst werden muesen.
  ** Bei Randrasterpunkten kann keine Interpolation erfolgen, so dass hier entsprechend extrapoliert wird.
  */

  x = (x - (double)column); // Residuum normiert auf Stützpunktabstand
  y = (y - (double)row);    // Residuum normiert auf Stützpunktabstand

  x -= 0.5;       // Mitte berücksichtigen
  y -= 0.5;       // Mitte berücksichtigen

  if (x < 0.0) { column--; x += 1.0;}  // linke Spalte zur Interpolation nutzen

  if (y < 0.0) { row--; y += 1.0;}     // untere Reihe zur Interpolation nutzen

  // Randrasterpunkte berücksichtigen: extrapolieren statt interpolieren
  if (column <= -1)  { x = -x; column = 0;}
  else if (column >= (m_internal->_totalColumns - 1)) { x += 1.0; column = m_internal->_totalColumns - 2;}

  if (row == -1)      { y = -y; row = 0;}
  else if (row >= (m_internal->_totalRows - 1)) { y += 1.0; row = m_internal->_totalRows - 2;}

  // Die Approximation:
  double substitute = elevation;
  double e0 = getOrSubstituteElevation(column, row, substitute);
  double e1 = getOrSubstituteElevation(column + 1, row, substitute);
  double e2 = getOrSubstituteElevation(column + 1, row + 1, substitute);
  double e3 = getOrSubstituteElevation(column, row + 1, substitute);

  double er = e0 + x * (e1 - e0) + y * (e3 - e0) + x * y * (e2 - e3 + e0 - e1);

  return er;   // floor(er+0.5);; runden herausgenommen 19.8.2008 JT
}  // double ElevationMap:: getElevation( double lon, double lat)


double ElevationMap::getOrSubstituteElevation(int column, int row, double substitute) const
{
  int elevation = getElevation(column, row);

  if (elevation != std::numeric_limits<int>::min())
    return elevation;

  int substCount = 0;
  double calcSubstitute = 0.0;

  elevation = getElevation(column + 1, row);
  if (elevation != std::numeric_limits<int>::min()) {substCount++;  calcSubstitute += elevation; }

  elevation = getElevation(column - 1, row);
  if (elevation != std::numeric_limits<int>::min()) {substCount++;  calcSubstitute += elevation; }

  elevation = getElevation(column, row + 1);
  if (elevation != std::numeric_limits<int>::min()) {substCount++;  calcSubstitute += elevation; }

  elevation = getElevation(column, row - 1);
  if (elevation != std::numeric_limits<int>::min()) {substCount++;  calcSubstitute += elevation; }

  if (substCount)
    substitute = calcSubstitute / substCount;

  return substitute;
} // double ElevationMap::getOrSubstituteElevation( int column, int row, double substitute)

double ElevationMap::distance(const QPointF &start, const QPointF &end) const
{
  if (!isValid())
    return 0;

  if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
  {
    return -1.0; // ausserhalb{ der Höhenkarte !!
  }

  double start_z = getElevation(start.x(), start.y());
  if (start_z <= std::numeric_limits<int>::min())    return -1.0;

  double  end_z = getElevation(end.x(), end.y());
  if (end_z <= std::numeric_limits<int>::min())    return -1.0;

  const QPointF &dist = distance_in_meter_WGS84(start, end);

  double dz = end_z - start_z;
  double d  = sqrt(dist.x() * dist.x() + dist.y() * dist.y() + dz * dz) ;

  return d;
} // double ElevationMap::distance( const QPointF &start, const QPointF &end )


double ElevationMap::profileDistance(const QPointF &start, const QPointF &end) const
{
  if (!isValid())
    return 0;

  if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
  {
    return -1.0; // ausserhalb der Höhenkarte !!
  }

  Vec3dVector elevProfile;
  getSamplePoints(start, end, elevProfile);

  for (Vec3dVector::iterator itr = elevProfile.begin(); itr != elevProfile.end(); ++itr)
    itr->setZ(getElevation(itr->x(), itr->y()));

  // Einzellaengen aufsummieren, Stützstellen mit undefinierten Hoehen werden ausgelassen!
  Point3D curPos(*elevProfile.begin());
  double sumLng = 0.0;
  for (Vec3dVector::iterator itr = elevProfile.begin() + 1; itr != elevProfile.end(); ++itr)
  {
    if ((*itr).z() > std::numeric_limits<int>::min())
    {
      if (curPos.z() > std::numeric_limits<int>::min())
      {
        const QPointF &dist = distance_in_meter_WGS84(QPointF(curPos.x(), curPos.y()), QPointF(itr->x(), itr->y()));
        double dz = itr->z() - curPos.z();
        sumLng += sqrt(dist.x() * dist.x() + dist.y() * dist.y() + dz * dz) ;
      }
      else
      {
        double dx = itr->x() - curPos.x();
        double dy = itr->y() - curPos.y();
        sumLng += sqrt(dx * dx + dy * dy);
      }
      curPos = *itr;
    }
  }
  return sumLng;
} // double ElevationMap::profileDistance( const QPointF &start, const QPointF &end)


unsigned int ElevationMap::elevationProfile(const QPointF &start, const QPointF &end, Vec3dVector &elevProfile) const
{
  unsigned int startSize = elevProfile.size();

  getSamplePoints(start, end, elevProfile);
  for (Vec3dVector::iterator itr = elevProfile.begin() + startSize; itr != elevProfile.end(); ++itr)
    itr->setZ(getElevation(itr->x(), itr->y()));

  return (elevProfile.size() - startSize);
} // unsigned int ElevationMap::elevationProfile( const QPointF &start, const QPointF &end, Vec3dVector &elevProfile)

bool ElevationMap::contains(const QPointF &point) const
{
  return m_internal->_location.contains(point);
}

bool ElevationMap::intersects(const QPolygonF &polygon) const
{
  const BoundingRect &bbox = m_internal->_location;

  // maybe QPointF uses float, but we can accept double to be reduced to float:
  const QPolygonF &ourBorder = QPolygonF({QPointF(bbox.left() , bbox.bottom()),
                                          QPointF(bbox.left() , bbox.top()),
                                          QPointF(bbox.right(), bbox.top()),
                                          QPointF(bbox.right(), bbox.bottom())});

#if QT_VERSION < QT_VERSION_CHECK(5,10,0)
  return ourBorder.intersected(polygon).size() > 0;
#else
  return ourBorder.intersects(polygon);
#endif
}

void ElevationMap::getSamplePoints(const QPointF &start, const QPointF &end, Vec3dVector &sPnts) const
{
  sPnts.clear();
  sPnts.push_back(Point3D::fromQPointF(start));

  bool endIsClipped = false;
  QPointF cStart(start);
  QPointF cEnd(end);

  if (!m_internal->_location.contains(start) || !m_internal->_location.contains(end))
  {
    int res = clipp2Map(start, end, cStart, cEnd);

    if (res == 0)  // no intersection
    {
      sPnts.push_back(Point3D::fromQPointF(end));
      return;
    }
    if (res & 1) // start clipped
      sPnts.push_back(Point3D::fromQPointF(cStart));
    if (res & 2) // end clipped
      endIsClipped = true;
  }

  QPointF delta(cEnd - cStart);

  // Berechne alle Schnittpunkte des Vektors( cStart, cEnd) mit den x-Grenzen der Hoehenstuetzrechtecke.
  // Die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Raster-Rechtecks
  int startCol = floor((cStart.x() - m_internal->_location.left()) / m_internal->_elexelSize.x());
  if (startCol < 0) startCol = 0;
  else if (startCol >= m_internal->_totalColumns) startCol = m_internal->_totalColumns - 1;

  int endCol   = floor((cEnd.x() - m_internal->_location.left()) / m_internal->_elexelSize.x());
  if (endCol < 0) endCol = 0;
  else if (endCol >= m_internal->_totalColumns) endCol = m_internal->_totalColumns - 1;

  int incCol = 0;
  if (endCol > startCol) incCol = 1;
  else if (endCol < startCol) incCol = -1;

  double m = delta.y() / delta.x();
  Vec2dVector xCuts;
  xCuts.reserve(std::abs(endCol - startCol));
  for (int col = startCol; col != endCol; col += incCol)
  {
    double x = m_internal->_location.left() + ((col + 0.5) + (incCol * 0.5)) * m_internal->_elexelSize.x();
    double y = cStart.y() + (x - cStart.x()) * m;
    xCuts.push_back(QPointF(x, y));
  }

  // Berechne alle Schnittpunkte des Vektors( cStart, cEnd) mit den y-Grenzen der Hoehenstuetzrechtecke.
  // Die Stützpunkte liegen definitiongemaess jeweils in der Mitte des Raster-Rechtecks
  int startRow = floor((cStart.y() - m_internal->_location.bottom()) / m_internal->_elexelSize.y());
  if (startRow < 0) startRow = 0;
  else if (startRow >= m_internal->_totalRows) startRow = m_internal->_totalRows - 1;

  int endRow   = floor((end.y()   - m_internal->_location.bottom()) / m_internal->_elexelSize.y());
  if (endRow < 0) endRow = 0;
  else if (endRow >= m_internal->_totalRows) endRow = m_internal->_totalRows - 1;

  int incRow = 0;
  if (endRow > startRow) incRow = 1;
  else if (endRow < startRow) incRow = -1;

  m = delta.x() / delta.y();
  Vec2dVector yCuts;
  yCuts.reserve(std::abs(endRow - startRow));
  for (int row = startRow; row != endRow; row += incRow)
  {
    double y = m_internal->_location.bottom() + ((row + 0.5) + (incRow * 0.5)) * m_internal->_elexelSize.y();
    double x = cStart.x() + (y - cStart.y()) * m;
    yCuts.push_back(QPointF(x, y));
  }

  // xCuts und yCuts zusammenführen, unter Berücksichtigung der Reihenfolge
  Vec2dVector mergeCuts;
  Vec2dVector::iterator xCutItr = xCuts.begin();
  Vec2dVector::iterator yCutItr = yCuts.begin();

  while ((xCutItr != xCuts.end()) && (yCutItr != yCuts.end()))
  {
    if (incCol > 0)
    {
      if ((*xCutItr).x() < (*yCutItr).x()) { mergeCuts.push_back(*xCutItr);  ++xCutItr; }
      else { mergeCuts.push_back(*yCutItr);  ++yCutItr; }
    }
    else
    {
      if ((*xCutItr).x() > (*yCutItr).x()) { mergeCuts.push_back(*xCutItr);  ++xCutItr; }
      else { mergeCuts.push_back(*yCutItr);  ++yCutItr; }
    }
  }
  while ((xCutItr != xCuts.end())) { mergeCuts.push_back(*xCutItr);  ++xCutItr; }
  while ((yCutItr != yCuts.end())) { mergeCuts.push_back(*yCutItr);  ++yCutItr; }

  // mergeCuts enthält sequenziell alle Schnittpunkte des Vektors mit den Hoehengridkanten
  // Die 2 Schnittpunkte pro Hoehenwert werden gemittelt zu einem zusammengefasst
  if (mergeCuts.size() > 1)
  {
    for (Vec2dVector::iterator itr = mergeCuts.begin() + 1; itr != mergeCuts.end(); ++itr)
      sPnts.push_back(Point3D::fromQPointF((*(itr - 1) + *itr) * 0.5));
  }
  sPnts.push_back(Point3D::fromQPointF(cEnd));

//  DEBUG_LOG(m_logger) << " xCuts.si " << xCuts.size()
//                      << ", yCuts.si " << yCuts.size()
//                      << ", mergeCuts.si " << mergeCuts.size();


  if (endIsClipped) sPnts.push_back(Point3D::fromQPointF(end));
} // void ElevationMap::getSamplePoints( const QPointF &start, const QPointF &end, Vec3dVector & sPnts)


// start und/oder end liegt ausserhalb !! in Lat/Lon !!
int ElevationMap::clipp2Map(const QPointF &start, const QPointF &end, QPointF &cStart, QPointF &cEnd) const
{
  Vec2dVector clippVec;

  if (end.x() != start.x())
  {
    double x = m_internal->_location.right();
    double m = (end.y() - start.y()) / (end.x() - start.x());
    if (isInIntervall(x, start.x(), end.x()))
    {
      double yc = start.y() + (x - start.x()) * m;
      if ((yc <= m_internal->_location.top()) && (yc >= m_internal->_location.bottom()))
        clippVec.push_back(QPointF(x, yc));
    }
    x = m_internal->_location.left();
    if (isInIntervall(x, start.x(), end.x()))
    {
      double yc = start.y() + (x - start.x()) * (end.y() - start.y()) / (end.x() - start.x());
      if ((yc <= m_internal->_location.top()) && (yc >= m_internal->_location.bottom()))
        clippVec.push_back(QPointF(x, yc));
    }
  }

  if (end.y() != start.y())
  {
    double y = m_internal->_location.top();
    double m = (end.x() - start.x()) / (end.y() - start.y());
    if (isInIntervall(y, start.y(), end.y()))
    {
      double xc = start.x() + (y - start.y()) * m;
      if ((xc < m_internal->_location.right()) && (xc > m_internal->_location.left()))
        clippVec.push_back(QPointF(xc, y));
    }
    y = m_internal->_location.bottom();
    if (isInIntervall(y, start.y(), end.y()))
    {
      double xc = start.x() + (y - start.y()) * m;
      if ((xc < m_internal->_location.right()) && (xc > m_internal->_location.left()))
        clippVec.push_back(QPointF(xc, y));
    }
  }
  int size = (int)clippVec.size();
  if (size == 0)
  {
    return 0;  // komlett ausserhalb ( oder komplett innerhalb)
  }

  if (size == 1)
  {
    if (m_internal->_location.contains(start))
    {
      cEnd = clippVec[0]; return 2;
    }
    else
    {
      cStart = clippVec[0]; return 1;
    }
  }

  if (size == 2)
  {
    if (squareLength(clippVec[0] - start) < squareLength(clippVec[1] - start))
    {
      cStart = clippVec[0];
      cEnd   = clippVec[1];
    }
    else
    {
      cStart = clippVec[1];
      cEnd   = clippVec[0];
    }
    return 3;
  }

  ERROR_LOG(m_logger) << LOGFUNCTION << "Internal Clipp Error!";

  return 0;
} // int ElevationMap::clipp2Map( const QPointF &start, const QPointF &end, QPointF &cStart, QPointF &cEnd)


// Line of Sight Profile: berücksichtigt Erdkrümmung und WGS84 Erdellipsoid
void ElevationMap::LOSProfile(const QPointF &start, const QPointF &end, LOSVector &losProfile)
{
  m_gdalReadWatcher.waitForFinished();

  if (!m_internal->_location.contains(start))    return ;  // ausserhalb der Höhenkarte !!

  bool clipped = false;
  double totalLng;
  QPointF startG = start;
  QPointF endG = end;

  if (!m_internal->_location.contains(end))  // Clippen
  {
    clipped = true;
    totalLng = length(endG - startG);
    if (endG.x() < m_internal->_location.left())
    {
      double m = (endG.y() - startG.y()) / (endG.x() - startG.x());
      double x = m_internal->_location.left();
      endG.setX(x);
      endG.setY(startG.y() + m * (x - startG.x()));
    }
    else if (endG.x() > m_internal->_location.right())
    {
      double m = (endG.y() - startG.y()) / (endG.x() - startG.x());
      double x = m_internal->_location.right();
      endG.setX(x);
      endG.setY(startG.y() + m * (x - startG.x()));
    }
    if (endG.y() < m_internal->_location.bottom())
    {
      double m = (endG.x() - startG.x()) / (endG.y() - startG.y());
      double y = m_internal->_location.bottom();
      endG.setY(y);
      endG.setX(startG.x() + m * (y - startG.y()));
    }
    else if (endG.y() > m_internal->_location.top())
    {
      double m = (endG.x() - startG.x()) / (endG.y() - startG.y());
      double y = m_internal->_location.top();
      endG.setY(y);
      endG.setX(startG.x() + m * (y - startG.y()));
    }
  }

  // Korrektursteigung aufgrund WGS84 Erdellipsoid, Lineare Approximation, Fehler bis 200 km: < 1.0m
  double startRadius = radiusWGS84(startG.y());
  double corrM = (radiusWGS84(endG.y()) - startRadius) / length(endG - startG);

  Vec3dVector elevProfile;
  getSamplePoints(start, end, elevProfile);
  for (Vec3dVector::iterator itr = elevProfile.begin(); itr != elevProfile.end(); ++itr)
  {
    double dist  = length(distance_in_meter_WGS84(start, itr->toPointF()));
    double height = elevation(itr->x(), itr->y());
    height += corrM * dist; // WGS84 Erdellipsoid Korrektur
    height -= (startRadius * (1 - cos(atan(dist / startRadius))));        // Einbeziehen der Erdkruemmung
    // alternativ als Naeherung: height -= ( dist * dist) / (startRadius * 2.0);
    losProfile.push_back(LOSentry(dist, height));
  }

  if (clipped)
    losProfile.push_back(LOSentry(totalLng, -NAN));

  qCritical() << " losProfile.isze() = " << losProfile.size();

} // void ElevationMap::LOSProfile( const QPointF &start, const QPointF &end, LOSVector &losProfile)



} // SituationMap



/* -----------------------------------------------------------------------------------------------
** CVS-Id/Release ...: $Header$
** -----------------------------------------------------------------------------------------------
**
** Historie:
**
**    $Log: not supported by cvs2svn $
**    Revision 1.17.1.1  2008/06/16 06:05:37  joerg
**    Transparenz und Bias einheitlich für alle Objekte implementiert
**
**    Revision 1.17  2006/03/24 09:05:01  klaus
**    Elevationprofile: Begrenzung auf  Kartengröße aufgehoben ( Clipping)
**
**    Revision 1.16  2005/11/09 13:44:36  klaus
**    Umsetzung auf Linux;
**
**    Revision 1.15  2005/04/15 13:05:45  joerg
**    Resource Control
**
**    Revision 1.14  2005/03/04 14:29:26  klaus
**    Elevationmap mit CachedTileList versehen und automatischem Entladen von Tiles
**
**    Revision 1.13  2005/03/03 12:37:53  klaus
**    dynamischen Laden und Entladen der ElevationTiles für Reliefmap
**
**    Revision 1.12  2005/03/01 12:36:59  klaus
**    Reliefmap eingeführt, aus Elevationmap entfernt
**
**    Revision 1.11  2005/02/25 15:32:56  klaus
**    Ergänzung um Reliefkarte
**
**    Revision 1.10  2005/02/10 14:15:19  klaus
**    Bias in Basisklasse MapObject überführt
**
**    Revision 1.9  2005/02/10 13:50:33  klaus
**    Bias für SightFans eingeführt
**
**    Revision 1.8  2005/02/09 16:45:26  klaus
**    1. Version der SightFans
**
**    Revision 1.7  2005/01/26 12:26:23  klaus
**    todos bearbeitet, Fehlerrückmeldnugen
**
**    Revision 1.6  2004/12/14 10:29:48  klaus
**    7. Erweiterung für Höhenkarten
**
**    Revision 1.5  2004/12/07 11:08:28  klaus
**    5. Erweiterung für Höhenkarten
**
**    Revision 1.4  2004/12/03 16:25:14  klaus
**    4. Erweiterung für Höhenkarten
**
**    Revision 1.3  2004/11/30 13:18:13  klaus
**    3. Erweiterung für Höhenkarten
**
**    Revision 1.2  2004/11/23 14:36:40  klaus
**    2. Erweiterung für Höhenkarten
**
**
** ----------------------------------------------------------------------------------------------- */