qgsmaptopixelgeometrysimplifier.cpp

Go to the documentation of this file.

/***************************************************************************
    qgsmaptopixelgeometrysimplifier.cpp
    ---------------------
    begin                : December 2013
    copyright            : (C) 2013 by Alvaro Huarte
    email                : http://wiki.osgeo.org/wiki/Alvaro_Huarte

 ***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <limits>
#include "qgsmaptopixelgeometrysimplifier.h"
#include "qgsapplication.h"

QgsMapToPixelSimplifier::QgsMapToPixelSimplifier( int simplifyFlags, double tolerance )
    : mSimplifyFlags( simplifyFlags )
    , mTolerance( tolerance )
{
}
QgsMapToPixelSimplifier::~QgsMapToPixelSimplifier()
{
}

// Helper simplification methods

float QgsMapToPixelSimplifier::calculateLengthSquared2D( double x1, double y1, double x2, double y2 )
{
  float vx = ( float )( x2 - x1 );
  float vy = ( float )( y2 - y1 );

  return vx*vx + vy*vy;
}

inline static QgsRectangle calculateBoundingBox( QGis::WkbType wkbType, unsigned char* wkb, size_t numPoints )
{
  double xmin =  std::numeric_limits<double>::max(), x, y;
  double ymin =  std::numeric_limits<double>::max();
  double xmax = -std::numeric_limits<double>::max();
  double ymax = -std::numeric_limits<double>::max();

  int sizeOfDoubleX = sizeof( double );
  int sizeOfDoubleY = QGis::wkbDimensions( wkbType ) == 3 /*hasZValue*/ ? 2 * sizeof( double ) : sizeof( double );

  for ( size_t i = 0; i < numPoints; ++i )
  {
    memcpy( &x, wkb, sizeof( double ) ); wkb += sizeOfDoubleX;
    memcpy( &y, wkb, sizeof( double ) ); wkb += sizeOfDoubleY;

    if ( xmin > x ) xmin = x;
    if ( ymin > y ) ymin = y;
    if ( xmax < x ) xmax = x;
    if ( ymax < y ) ymax = y;
  }

  return QgsRectangle( xmin, ymin, xmax, ymax );
}

inline static bool generalizeWkbGeometry( QGis::WkbType wkbType, unsigned char* sourceWkb, size_t sourceWkbSize, unsigned char* targetWkb, size_t& targetWkbSize, const QgsRectangle& envelope, bool writeHeader )
{
  Q_UNUSED( sourceWkb );
  unsigned char* wkb2 = targetWkb;
  unsigned int geometryType = QGis::singleType( QGis::flatType( wkbType ) );

  int sizeOfDoubleX = sizeof( double );
  int sizeOfDoubleY = QGis::wkbDimensions( wkbType ) == 3 /*hasZValue*/ ? 2 * sizeof( double ) : sizeof( double );

  // Skip the unnecesary generalization because of is a very single geometry
  size_t minimumSize = ( geometryType == QGis::WKBLineString ? 4 + 2 * ( sizeOfDoubleX + sizeOfDoubleY ) : 8 + 5 * ( sizeOfDoubleX + sizeOfDoubleY ) );
  if ( writeHeader ) minimumSize += 5;
  if ( sourceWkbSize <= minimumSize )
  {
    targetWkbSize = 0;
    return false;
  }

  double x1 = envelope.xMinimum();
  double y1 = envelope.yMinimum();
  double x2 = envelope.xMaximum();
  double y2 = envelope.yMaximum();

  // Write the main header of the geometry
  if ( writeHeader )
  {
    char byteOrder = QgsApplication::endian(); // byteOrder
    memcpy( targetWkb, &byteOrder, 1 );
    targetWkb += 1;

    memcpy( targetWkb, &geometryType, 4 ); // type
    targetWkb += 4;

    if ( geometryType == QGis::WKBPolygon ) // numRings
    {
      int numRings = 1;
      memcpy( targetWkb, &numRings, 4 );
      targetWkb += 4;
    }
  }

  // Write the generalized geometry
  if ( geometryType == QGis::WKBLineString )
  {
    int numPoints = 2;
    memcpy( targetWkb, &numPoints, 4 ); // numPoints;
    targetWkb += 4;

    memcpy( targetWkb, &x1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &x2, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y2, sizeof( double ) ); targetWkb += sizeof( double );
  }
  else
  {
    int numPoints = 5;
    memcpy( targetWkb, &numPoints, 4 ); // numPoints;
    targetWkb += 4;

    memcpy( targetWkb, &x1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &x2, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &x2, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y2, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &x1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y2, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &x1, sizeof( double ) ); targetWkb += sizeof( double );
    memcpy( targetWkb, &y1, sizeof( double ) ); targetWkb += sizeof( double );
  }
  targetWkbSize += targetWkb - wkb2;

  return true;
}

bool QgsMapToPixelSimplifier::simplifyWkbGeometry( int simplifyFlags, QGis::WkbType wkbType, unsigned char* sourceWkb, size_t sourceWkbSize, unsigned char* targetWkb, size_t& targetWkbSize, const QgsRectangle& envelope, double map2pixelTol, bool writeHeader, bool isaLinearRing )
{
  bool canbeGeneralizable = true;
  bool hasZValue = QGis::wkbDimensions( wkbType ) == 3;
  bool result = false;

  // Can replace the geometry by its BBOX ?
  if (( simplifyFlags & QgsMapToPixelSimplifier::SimplifyEnvelope ) && ( envelope.xMaximum() - envelope.xMinimum() ) < map2pixelTol && ( envelope.yMaximum() - envelope.yMinimum() ) < map2pixelTol )
  {
    canbeGeneralizable = generalizeWkbGeometry( wkbType, sourceWkb, sourceWkbSize, targetWkb, targetWkbSize, envelope, writeHeader );
    if ( canbeGeneralizable ) return true;
  }
  if ( !( simplifyFlags & QgsMapToPixelSimplifier::SimplifyGeometry ) ) canbeGeneralizable = false;

  // Write the main header of the geometry
  if ( writeHeader )
  {
    memcpy( targetWkb, sourceWkb, 1 ); // byteOrder
    sourceWkb += 1;
    targetWkb += 1;

    int geometryType;
    memcpy( &geometryType, sourceWkb, 4 );
    int flatType = QGis::flatType(( QGis::WkbType )geometryType );
    memcpy( targetWkb, &flatType, 4 ); // type
    sourceWkb += 4;
    targetWkb += 4;

    targetWkbSize += 5;
  }

  unsigned char* wkb1 = sourceWkb;
  unsigned char* wkb2 = targetWkb;
  unsigned int flatType = QGis::flatType( wkbType );

  // Write the geometry
  if ( flatType == QGis::WKBLineString || isaLinearRing )
  {
    double x, y, lastX = 0, lastY = 0;

    int sizeOfDoubleX = sizeof( double );
    int sizeOfDoubleY = QGis::wkbDimensions( wkbType ) == 3 /*hasZValue*/ ? 2 * sizeof( double ) : sizeof( double );

    int numPoints;
    memcpy( &numPoints, sourceWkb, 4 );
    sourceWkb += 4;
    if ( numPoints <= ( isaLinearRing ? 5 : 2 ) ) canbeGeneralizable = false;

    int numTargetPoints = 0;
    memcpy( targetWkb, &numTargetPoints, 4 );
    targetWkb += 4;
    targetWkbSize += 4;

    double* ptr = ( double* )targetWkb;
    map2pixelTol *= map2pixelTol; //-> Use mappixelTol for 'LengthSquare' calculations.

    // Process each vertex...
    for ( int i = 0, numPoints_i = ( isaLinearRing ? numPoints - 1 : numPoints ); i < numPoints_i; ++i )
    {
      memcpy( &x, sourceWkb, sizeof( double ) ); sourceWkb += sizeOfDoubleX;
      memcpy( &y, sourceWkb, sizeof( double ) ); sourceWkb += sizeOfDoubleY;

      if ( i == 0 || !canbeGeneralizable || QgsMapToPixelSimplifier::calculateLengthSquared2D( x, y, lastX, lastY ) > map2pixelTol || ( !isaLinearRing && ( i == 1 || i >= numPoints - 2 ) ) )
      {
        memcpy( ptr, &x, sizeof( double ) ); lastX = x; ptr++;
        memcpy( ptr, &y, sizeof( double ) ); lastY = y; ptr++;
        numTargetPoints++;
      }
    }
    targetWkb = wkb2 + 4;

    // Fix the topology of the geometry
    if ( isaLinearRing )
    {
      memcpy( &x, targetWkb + 0, sizeof( double ) );
      memcpy( &y, targetWkb + 8, sizeof( double ) );
      memcpy( ptr, &x, sizeof( double ) ); ptr++;
      memcpy( ptr, &y, sizeof( double ) ); ptr++;
      numTargetPoints++;
    }
    targetWkbSize += numTargetPoints * 16;
    targetWkb = wkb2;

    memcpy( targetWkb, &numTargetPoints, 4 );
    result = numPoints != numTargetPoints;
  }
  else
    if ( flatType == QGis::WKBPolygon )
    {
      int numRings;
      memcpy( &numRings, sourceWkb, 4 );
      sourceWkb += 4;

      memcpy( targetWkb, &numRings, 4 );
      targetWkb += 4;
      targetWkbSize += 4;

      for ( int i = 0; i < numRings; ++i )
      {
        int numPoints_i;
        memcpy( &numPoints_i, sourceWkb, 4 );
        QgsRectangle envelope_i = numRings == 1 ? envelope : calculateBoundingBox( wkbType, sourceWkb + 4, numPoints_i );

        size_t sourceWkbSize_i = 4 + numPoints_i * ( hasZValue ? 3 : 2 ) * sizeof( double );
        size_t targetWkbSize_i = 0;

        result |= simplifyWkbGeometry( simplifyFlags, wkbType, sourceWkb, sourceWkbSize_i, targetWkb, targetWkbSize_i, envelope_i, map2pixelTol, false, true );
        sourceWkb += sourceWkbSize_i;
        targetWkb += targetWkbSize_i;

        targetWkbSize += targetWkbSize_i;
      }
    }
    else
      if ( flatType == QGis::WKBMultiLineString || flatType == QGis::WKBMultiPolygon )
      {
        int numGeoms;
        memcpy( &numGeoms, sourceWkb, 4 );
        sourceWkb += 4;
        wkb1 += 4;

        memcpy( targetWkb, &numGeoms, 4 );
        targetWkb += 4;
        targetWkbSize += 4;

        for ( int i = 0; i < numGeoms; ++i )
        {
          size_t sourceWkbSize_i = 0;
          size_t targetWkbSize_i = 0;

          // ... calculate the wkb-size of the current child complex geometry
          if ( flatType == QGis::WKBMultiLineString )
          {
            int numPoints_i;
            memcpy( &numPoints_i, wkb1 + 5, 4 );
            int wkbSize_i = 4 + numPoints_i * ( hasZValue ? 3 : 2 ) * sizeof( double );

            sourceWkbSize_i += 5 + wkbSize_i;
            wkb1 += 5 + wkbSize_i;
          }
          else
          {
            int numPrings_i;
            memcpy( &numPrings_i, wkb1 + 5, 4 );
            sourceWkbSize_i = 9;
            wkb1 += 9;

            for ( int j = 0; j < numPrings_i; ++j )
            {
              int numPoints_i;
              memcpy( &numPoints_i, wkb1, 4 );
              int wkbSize_i = 4 + numPoints_i * ( hasZValue ? 3 : 2 ) * sizeof( double );

              sourceWkbSize_i += wkbSize_i;
              wkb1 += wkbSize_i;
            }
          }
          result |= simplifyWkbGeometry( simplifyFlags, QGis::singleType( wkbType ), sourceWkb, sourceWkbSize_i, targetWkb, targetWkbSize_i, envelope, map2pixelTol, true, false );
          sourceWkb += sourceWkbSize_i;
          targetWkb += targetWkbSize_i;

          targetWkbSize += targetWkbSize_i;
        }
      }
  return result;
}


bool QgsMapToPixelSimplifier::canbeGeneralizedByMapBoundingBox( const QgsRectangle& envelope, double map2pixelTol )
{
  // Can replace the geometry by its BBOX ?
  if (( envelope.xMaximum() - envelope.xMinimum() ) < map2pixelTol && ( envelope.yMaximum() - envelope.yMinimum() ) < map2pixelTol )
  {
    return true;
  }
  return false;
}

QgsGeometry* QgsMapToPixelSimplifier::simplify( QgsGeometry* geometry ) const
{
  QgsGeometry* g = new QgsGeometry();

  size_t wkbSize = geometry->wkbSize();
  unsigned char* wkb = ( unsigned char* )malloc( wkbSize );
  memcpy( wkb, geometry->asWkb(), wkbSize );
  g->fromWkb( wkb, wkbSize );
  simplifyGeometry( g, mSimplifyFlags, mTolerance );

  return g;
}

bool QgsMapToPixelSimplifier::simplifyGeometry( QgsGeometry* geometry, int simplifyFlags, double tolerance )
{
  size_t targetWkbSize = 0;

  // Check whether the geometry can be simplified using the map2pixel context
  QGis::GeometryType geometryType = geometry->type();
  if ( !( geometryType == QGis::Line || geometryType == QGis::Polygon ) ) return false;

  QgsRectangle envelope = geometry->boundingBox();
  QGis::WkbType wkbType = geometry->wkbType();

  unsigned char* wkb = ( unsigned char* )geometry->asWkb( );
  size_t wkbSize = geometry->wkbSize( );

  // Simplify the geometry rewriting temporally its WKB-stream for saving calloc's.
  if ( simplifyWkbGeometry( simplifyFlags, wkbType, wkb, wkbSize, wkb, targetWkbSize, envelope, tolerance ) )
  {
    unsigned char* targetWkb = ( unsigned char* )malloc( targetWkbSize );
    memcpy( targetWkb, wkb, targetWkbSize );
    geometry->fromWkb( targetWkb, targetWkbSize );
    return true;
  }
  return false;
}

bool QgsMapToPixelSimplifier::simplifyGeometry( QgsGeometry* geometry ) const
{
  return simplifyGeometry( geometry, mSimplifyFlags, mTolerance );
}