Move heatmap generation code to analysis lib

And clean it up a lot

python/analysis/analysis.sip

@@ -44,6 +44,7 @@
 %Include raster/qgsderivativefilter.sip
 %Include raster/qgsaspectfilter.sip
 %Include raster/qgshillshadefilter.sip
+%Include raster/qgskde.sip
 %Include raster/qgsninecellfilter.sip
 %Include raster/qgsrastercalcnode.sip
 %Include raster/qgsrastercalculator.sip

python/analysis/raster/qgskde.sip

@@ -0,0 +1,104 @@
+/**
+ * \class QgsKernelDensityEstimation
+ * \ingroup analysis
+ * Performs Kernel Density Estimation ("heatmap") calculations on a vector layer.
+ * @note added in QGIS 3.0
+ */
+class QgsKernelDensityEstimation
+{
+%TypeHeaderCode
+#include <qgskde.h>
+%End
+
+  public:
+
+    //! Kernel shape type
+    enum KernelShape
+    {
+      KernelQuartic, //!< Quartic kernel
+      KernelTriangular, //!< Triangular kernel
+      KernelUniform, //!< Uniform (flat) kernel
+      KernelTriweight, //!< Triweight kernel
+      KernelEpanechnikov, //!< Epanechnikov kernel
+    };
+
+    //! Output values type
+    enum OutputValues
+    {
+      OutputRaw, //!< Output the raw KDE values
+      OutputScaled, //!< Output mathematically correct scaled values
+    };
+
+    //! Result of operation
+    enum Result
+    {
+      Success, //!< Operation completed successfully
+      DriverError, //!< Could not open the driver for the specified format
+      InvalidParameters, //!< Input parameters were not valid
+      FileCreationError, //!< Error creating output file
+      RasterIoError, //!< Error writing to raster
+    };
+
+    //! KDE parameters
+    struct Parameters
+    {
+      //! Vector point layer
+      QgsVectorLayer* vectorLayer;
+
+      //! Fixed radius, in map units
+      double radius;
+
+      //! Field for radius, or empty if using a fixed radius
+      QString radiusField;
+
+      //! Field name for weighting field, or empty if not using weights
+      QString weightField;
+
+      //! Size of pixel in output file
+      double pixelSize;
+
+      //! Kernel shape
+      QgsKernelDensityEstimation::KernelShape shape;
+
+      //! Decay ratio (Triangular kernels only)
+      double decayRatio;
+
+      //! Type of output value
+      QgsKernelDensityEstimation::OutputValues outputValues;
+    };
+
+    /**
+     * Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use
+     * to generate the surface. The output path and file format are also required.
+     */
+    QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat );
+
+    /**
+     * Runs the KDE calculation across the whole layer at once. Either call this method, or manually
+     * call run(), addFeature() and finalise() separately.
+     */
+    Result run();
+
+    /**
+     * Prepares the output file for writing and setups up the surface calculation. This must be called
+     * before adding features via addFeature().
+     * @see addFeature()
+     * @see finalise()
+     */
+    Result prepare();
+
+    /**
+     * Adds a single feature to the KDE surface. prepare() must be called before adding features.
+     * @see prepare()
+     * @see finalise()
+     */
+    Result addFeature( const QgsFeature& feature );
+
+    /**
+     * Finalises the output file. Must be called after adding all features via addFeature().
+     * @see prepare()
+     * @see addFeature()
+     */
+    Result finalise();
+
+};

src/analysis/CMakeLists.txt

@@ -27,6 +27,7 @@ SET(QGIS_ANALYSIS_SRCS
   raster/qgsruggednessfilter.cpp
   raster/qgsderivativefilter.cpp
   raster/qgshillshadefilter.cpp
+  raster/qgskde.cpp
   raster/qgsslopefilter.cpp
   raster/qgsaspectfilter.cpp
   raster/qgstotalcurvaturefilter.cpp
@@ -106,6 +107,7 @@ SET(QGIS_ANALYSIS_HDRS
   raster/qgsaspectfilter.h
   raster/qgsderivativefilter.h
   raster/qgshillshadefilter.h
+  raster/qgskde.h
   raster/qgsninecellfilter.h
   raster/qgsrastercalculator.h
   raster/qgsrelief.h

src/analysis/raster/qgskde.cpp

@@ -0,0 +1,448 @@
+/***************************************************************************
+  qgskde.cpp
+  ----------
+  Date                 : October 2016
+  Copyright            : (C) 2016 by Nyall Dawson
+  Email                : nyall dot dawson at gmail dot com
+ ***************************************************************************
+ *                                                                         *
+ *   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 "qgskde.h"
+#include "qgsvectorlayer.h"
+#include "qgsgeometry.h"
+
+#define NO_DATA -9999
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#if defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1800
+#define TO8F(x) (x).toUtf8().constData()
+#else
+#define TO8F(x) QFile::encodeName( x ).constData()
+#endif
+
+QgsKernelDensityEstimation::QgsKernelDensityEstimation( const QgsKernelDensityEstimation::Parameters& parameters, const QString& outputFile, const QString& outputFormat )
+    : mInputLayer( parameters.vectorLayer )
+    , mOutputFile( outputFile )
+    , mOutputFormat( outputFormat )
+    , mRadiusField( -1 )
+    , mWeightField( -1 )
+    , mRadius( parameters.radius )
+    , mPixelSize( parameters.pixelSize )
+    , mShape( parameters.shape )
+    , mDecay( parameters.decayRatio )
+    , mOutputValues( parameters.outputValues )
+    , mBufferSize( -1 )
+    , mDatasetH( nullptr )
+    , mRasterBandH( nullptr )
+{
+  if ( !parameters.radiusField.isEmpty() )
+    mRadiusField = mInputLayer->fields().lookupField( parameters.radiusField );
+  if ( !parameters.weightField.isEmpty() )
+    mWeightField = mInputLayer->fields().lookupField( parameters.weightField );
+}
+
+QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::run()
+{
+  Result result = prepare();
+  if ( result != Success )
+    return result;
+
+  QgsAttributeList requiredAttributes;
+
+  if ( mRadiusField >= 0 )
+    requiredAttributes << mRadiusField;
+
+  if ( mWeightField >= 0 )
+    requiredAttributes << mWeightField;
+
+  QgsFeatureIterator fit = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( requiredAttributes ) );
+
+  QgsFeature f;
+  while ( fit.nextFeature( f ) )
+  {
+    addFeature( f );
+  }
+
+  return finalise();
+}
+
+QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::prepare()
+{
+  GDALAllRegister();
+
+  GDALDriverH driver = GDALGetDriverByName( mOutputFormat.toUtf8() );
+  if ( !driver )
+  {
+    return DriverError;
+  }
+
+  if ( !mInputLayer )
+    return InvalidParameters;
+
+  mBounds = calculateBounds();
+  if ( mBounds.isNull() )
+    return InvalidParameters;
+
+  int rows = qMax( qRound( mBounds.height() / mPixelSize ), 1 );
+  int cols = qMax( qRound( mBounds.width() / mPixelSize ), 1 );
+
+  if ( !createEmptyLayer( driver, mBounds, rows, cols ) )
+    return FileCreationError;
+
+  // open the raster in GA_Update mode
+  mDatasetH = GDALOpen( TO8F( mOutputFile ), GA_Update );
+  if ( !mDatasetH )
+    return FileCreationError;
+  mRasterBandH = GDALGetRasterBand( mDatasetH, 1 );
+  if ( !mRasterBandH )
+    return FileCreationError;
+
+  mBufferSize = -1;
+  if ( mRadiusField < 0 )
+    mBufferSize = radiusSizeInPixels( mRadius );
+
+  return Success;
+}
+
+QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::addFeature( const QgsFeature& feature )
+{
+  QgsGeometry featureGeometry = feature.geometry();
+  if ( featureGeometry.isEmpty() )
+  {
+    return Success;
+  }
+
+  // convert the geometry to multipoint
+  QgsMultiPoint multiPoints;
+  if ( !featureGeometry.isMultipart() )
+  {
+    QgsPoint p = featureGeometry.asPoint();
+    // avoiding any empty points or out of extent points
+    if ( !mBounds.contains( p ) )
+      return Success;
+
+    multiPoints << p;
+  }
+  else
+  {
+    multiPoints = featureGeometry.asMultiPoint();
+  }
+
+  // if radius is variable then fetch it and calculate new pixel buffer size
+  double radius = mRadius;
+  int buffer = mBufferSize;
+  if ( mRadiusField >= 0 )
+  {
+    radius = feature.attribute( mRadiusField ).toDouble();
+    buffer = radiusSizeInPixels( radius );
+  }
+  int blockSize = 2 * buffer + 1; //Block SIDE would be more appropriate
+
+  // calculate weight
+  double weight = 1.0;
+  if ( mWeightField >= 0 )
+  {
+    weight = feature.attribute( mWeightField ).toDouble();
+  }
+
+  Result result = Success;
+
+  //loop through all points in multipoint
+  for ( QgsMultiPoint::const_iterator pointIt = multiPoints.constBegin(); pointIt != multiPoints.constEnd(); ++pointIt )
+  {
+    // avoiding any empty points or out of extent points
+    if ( !mBounds.contains( *pointIt ) )
+    {
+      continue;
+    }
+
+    // calculate the pixel position
+    unsigned int xPosition = ((( *pointIt ).x() - mBounds.xMinimum() ) / mPixelSize ) - buffer;
+    unsigned int yPosition = ((( *pointIt ).y() - mBounds.yMinimum() ) / mPixelSize ) - buffer;
+
+    // get the data
+    float *dataBuffer = ( float * ) CPLMalloc( sizeof( float ) * blockSize * blockSize );
+    if ( GDALRasterIO( mRasterBandH, GF_Read, xPosition, yPosition, blockSize, blockSize,
+                       dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
+    {
+      result = RasterIoError;
+    }
+
+    for ( int xp = 0; xp <= buffer; xp++ )
+    {
+      for ( int yp = 0; yp <= buffer; yp++ )
+      {
+        double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
+
+        // is pixel outside search bandwidth of feature?
+        if ( distance > buffer )
+        {
+          continue;
+        }
+
+        double pixelValue = weight * calculateKernelValue( distance, buffer, mShape, mOutputValues );
+
+        // clearing anamolies along the axes
+        if ( xp == 0 && yp == 0 )
+        {
+          pixelValue /= 4;
+        }
+        else if ( xp == 0 || yp == 0 )
+        {
+          pixelValue /= 2;
+        }
+
+        int pos[4];
+        pos[0] = ( buffer + xp ) * blockSize + ( buffer + yp );
+        pos[1] = ( buffer + xp ) * blockSize + ( buffer - yp );
+        pos[2] = ( buffer - xp ) * blockSize + ( buffer + yp );
+        pos[3] = ( buffer - xp ) * blockSize + ( buffer - yp );
+        for ( int p = 0; p < 4; p++ )
+        {
+          if ( dataBuffer[ pos[p] ] == NO_DATA )
+          {
+            dataBuffer[ pos[p] ] = 0;
+          }
+          dataBuffer[ pos[p] ] += pixelValue;
+        }
+      }
+    }
+    if ( GDALRasterIO( mRasterBandH, GF_Write, xPosition, yPosition, blockSize, blockSize,
+                       dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 ) != CE_None )
+    {
+      result = RasterIoError;
+    }
+    CPLFree( dataBuffer );
+  }
+
+  return result;
+}
+
+QgsKernelDensityEstimation::Result QgsKernelDensityEstimation::finalise()
+{
+  GDALClose(( GDALDatasetH ) mDatasetH );
+  mDatasetH = nullptr;
+  mRasterBandH = nullptr;
+  return Success;
+}
+
+int QgsKernelDensityEstimation::radiusSizeInPixels( double radius ) const
+{
+  int buffer = radius / mPixelSize;
+  if ( radius - ( mPixelSize * buffer ) > 0.5 )
+  {
+    ++buffer;
+  }
+  return buffer;
+}
+
+bool QgsKernelDensityEstimation::createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const
+{
+  double geoTransform[6] = { bounds.xMinimum(), mPixelSize, 0, bounds.yMinimum(), 0, mPixelSize };
+  GDALDatasetH emptyDataset = GDALCreate( driver, mOutputFile.toUtf8(), columns, rows, 1, GDT_Float32, nullptr );
+  if ( !emptyDataset )
+    return false;
+
+  if ( GDALSetGeoTransform( emptyDataset, geoTransform ) != CE_None )
+    return false;
+
+  // Set the projection on the raster destination to match the input layer
+  if ( GDALSetProjection( emptyDataset, mInputLayer->crs().toWkt().toLocal8Bit().data() ) != CE_None )
+    return false;
+
+  GDALRasterBandH poBand = GDALGetRasterBand( emptyDataset, 1 );
+  if ( !poBand )
+    return false;
+
+  if ( GDALSetRasterNoDataValue( poBand, NO_DATA ) != CE_None )
+    return false;
+
+  float* line = static_cast< float * >( CPLMalloc( sizeof( float ) * columns ) );
+  for ( int i = 0; i < columns; i++ )
+  {
+    line[i] = NO_DATA;
+  }
+  // Write the empty raster
+  for ( int i = 0; i < rows ; i++ )
+  {
+    if ( GDALRasterIO( poBand, GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 ) != CE_None )
+    {
+      return false;
+    }
+  }
+
+  CPLFree( line );
+  //close the dataset
+  GDALClose( emptyDataset );
+  return true;
+}
+
+double QgsKernelDensityEstimation::calculateKernelValue( const double distance, const int bandwidth, const QgsKernelDensityEstimation::KernelShape shape, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  switch ( shape )
+  {
+    case KernelTriangular:
+      return triangularKernel( distance, bandwidth, outputType );
+
+    case KernelUniform:
+      return uniformKernel( distance, bandwidth, outputType );
+
+    case KernelQuartic:
+      return quarticKernel( distance, bandwidth, outputType );
+
+    case KernelTriweight:
+      return triweightKernel( distance, bandwidth, outputType );
+
+    case KernelEpanechnikov:
+      return epanechnikovKernel( distance, bandwidth, outputType );
+  }
+  return 0; //no warnings
+}
+
+/* The kernel functions below are taken from "Kernel Smoothing" by Wand and Jones (1995), p. 175
+ *
+ * Each kernel is multiplied by a normalizing constant "k", which normalizes the kernel area
+ * to 1 for a given bandwidth size.
+ *
+ * k is calculated by polar double integration of the kernel function
+ * between a radius of 0 to the specified bandwidth and equating the area to 1. */
+
+double QgsKernelDensityEstimation::uniformKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  Q_UNUSED( distance );
+  switch ( outputType )
+  {
+    case OutputScaled:
+    {
+      // Normalizing constant
+      double k = 2. / ( M_PI * ( double )bandwidth );
+
+      // Derived from Wand and Jones (1995), p. 175
+      return k * ( 0.5 / ( double )bandwidth );
+    }
+    case OutputRaw:
+      return 1.0;
+  }
+  return 0.0; // NO warnings!!!!!
+}
+
+double QgsKernelDensityEstimation::quarticKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  switch ( outputType )
+  {
+    case OutputScaled:
+    {
+      // Normalizing constant
+      double k = 116. / ( 5. * M_PI * pow(( double )bandwidth, 2 ) );
+
+      // Derived from Wand and Jones (1995), p. 175
+      return k * ( 15. / 16. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 );
+    }
+    case OutputRaw:
+      return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 2 );
+  }
+  return 0.0; //no, seriously, I told you NO WARNINGS!
+}
+
+double QgsKernelDensityEstimation::triweightKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  switch ( outputType )
+  {
+    case OutputScaled:
+    {
+      // Normalizing constant
+      double k = 128. / ( 35. * M_PI * pow(( double )bandwidth, 2 ) );
+
+      // Derived from Wand and Jones (1995), p. 175
+      return k * ( 35. / 32. ) * pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 );
+    }
+    case OutputRaw:
+      return pow( 1. - pow( distance / ( double )bandwidth, 2 ), 3 );
+  }
+  return 0.0; // this is getting ridiculous... don't you ever listen to a word I say?
+}
+
+double QgsKernelDensityEstimation::epanechnikovKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  switch ( outputType )
+  {
+    case OutputScaled:
+    {
+      // Normalizing constant
+      double k = 8. / ( 3. * M_PI * pow(( double )bandwidth, 2 ) );
+
+      // Derived from Wand and Jones (1995), p. 175
+      return k * ( 3. / 4. ) * ( 1. - pow( distance / ( double )bandwidth, 2 ) );
+    }
+    case OutputRaw:
+      return ( 1. - pow( distance / ( double )bandwidth, 2 ) );
+  }
+
+  return 0.0; // la la la i'm not listening
+}
+
+double QgsKernelDensityEstimation::triangularKernel( const double distance, const int bandwidth, const QgsKernelDensityEstimation::OutputValues outputType ) const
+{
+  switch ( outputType )
+  {
+    case OutputScaled:
+    {
+      // Normalizing constant. In this case it's calculated a little different
+      // due to the inclusion of the non-standard "decay" parameter
+
+      if ( mDecay >= 0 )
+      {
+        double k = 3. / (( 1. + 2. * mDecay ) * M_PI * pow(( double )bandwidth, 2 ) );
+
+        // Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter)
+        return k * ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
+      }
+      else
+      {
+        // Non-standard or mathematically valid negative decay ("coolmap")
+        return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
+      }
+    }
+    case OutputRaw:
+      return ( 1. - ( 1. - mDecay ) * ( distance / ( double )bandwidth ) );
+  }
+  return 0.0; // ....
+}
+
+QgsRectangle QgsKernelDensityEstimation::calculateBounds() const
+{
+  if ( !mInputLayer )
+    return QgsRectangle();
+
+  QgsRectangle bbox = mInputLayer->extent();
+
+  double radius = 0;
+  if ( mRadiusField >= 0 )
+  {
+    // if radius is using a field, find the max value
+    radius = mInputLayer->maximumValue( mRadiusField ).toDouble();
+  }
+  else
+  {
+    radius = mRadius;
+  }
+  // expand the bounds by the maximum search radius
+  bbox.setXMinimum( bbox.xMinimum() - radius );
+  bbox.setYMinimum( bbox.yMinimum() - radius );
+  bbox.setXMaximum( bbox.xMaximum() + radius );
+  bbox.setYMaximum( bbox.yMaximum() + radius );
+  return bbox;
+}
+
+
+
+

src/analysis/raster/qgskde.h

@@ -0,0 +1,174 @@
+/***************************************************************************
+  qgskde.h
+  --------
+  Date                 : October 2016
+  Copyright            : (C) 2016 by Nyall Dawson
+  Email                : nyall dot dawson at gmail dot com
+ ***************************************************************************
+ *                                                                         *
+ *   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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSKDE_H
+#define QGSKDE_H
+
+#include "qgsrectangle.h"
+#include <QString>
+
+// GDAL includes
+#include <gdal.h>
+#include <cpl_string.h>
+#include <cpl_conv.h>
+
+class QgsVectorLayer;
+class QProgressDialog;
+class QgsFeature;
+
+
+/**
+ * \class QgsKernelDensityEstimation
+ * \ingroup analysis
+ * Performs Kernel Density Estimation ("heatmap") calculations on a vector layer.
+ * @note added in QGIS 3.0
+ */
+class ANALYSIS_EXPORT QgsKernelDensityEstimation
+{
+  public:
+
+    //! Kernel shape type
+    enum KernelShape
+    {
+      KernelQuartic = 0, //!< Quartic kernel
+      KernelTriangular, //!< Triangular kernel
+      KernelUniform, //!< Uniform (flat) kernel
+      KernelTriweight, //!< Triweight kernel
+      KernelEpanechnikov, //!< Epanechnikov kernel
+    };
+
+    //! Output values type
+    enum OutputValues
+    {
+      OutputRaw = 0, //!< Output the raw KDE values
+      OutputScaled, //!< Output mathematically correct scaled values
+    };
+
+    //! Result of operation
+    enum Result
+    {
+      Success, //!< Operation completed successfully
+      DriverError, //!< Could not open the driver for the specified format
+      InvalidParameters, //!< Input parameters were not valid
+      FileCreationError, //!< Error creating output file
+      RasterIoError, //!< Error writing to raster
+    };
+
+    //! KDE parameters
+    struct Parameters
+    {
+      //! Vector point layer
+      QgsVectorLayer* vectorLayer;
+
+      //! Fixed radius, in map units
+      double radius;
+
+      //! Field for radius, or empty if using a fixed radius
+      QString radiusField;
+
+      //! Field name for weighting field, or empty if not using weights
+      QString weightField;
+
+      //! Size of pixel in output file
+      double pixelSize;
+
+      //! Kernel shape
+      KernelShape shape;
+
+      //! Decay ratio (Triangular kernels only)
+      double decayRatio;
+
+      //! Type of output value
+      OutputValues outputValues;
+    };
+
+    /**
+     * Constructor for QgsKernelDensityEstimation. Requires a Parameters object specifying the options to use
+     * to generate the surface. The output path and file format are also required.
+     */
+    QgsKernelDensityEstimation( const Parameters& parameters, const QString& outputFile, const QString& outputFormat );
+
+    /**
+     * Runs the KDE calculation across the whole layer at once. Either call this method, or manually
+     * call run(), addFeature() and finalise() separately.
+     */
+    Result run();
+
+    /**
+     * Prepares the output file for writing and setups up the surface calculation. This must be called
+     * before adding features via addFeature().
+     * @see addFeature()
+     * @see finalise()
+     */
+    Result prepare();
+
+    /**
+     * Adds a single feature to the KDE surface. prepare() must be called before adding features.
+     * @see prepare()
+     * @see finalise()
+     */
+    Result addFeature( const QgsFeature& feature );
+
+    /**
+     * Finalises the output file. Must be called after adding all features via addFeature().
+     * @see prepare()
+     * @see addFeature()
+     */
+    Result finalise();
+
+  private:
+
+    //! Calculate the value given to a point width a given distance for a specified kernel shape
+    double calculateKernelValue( const double distance, const int bandwidth, const KernelShape shape, const OutputValues outputType ) const;
+    //! Uniform kernel function
+    double uniformKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
+    //! Quartic kernel function
+    double quarticKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
+    //! Triweight kernel function
+    double triweightKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
+    //! Epanechnikov kernel function
+    double epanechnikovKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
+    //! Triangular kernel function
+    double triangularKernel( const double distance, const int bandwidth, const OutputValues outputType ) const;
+
+    QgsRectangle calculateBounds() const;
+
+    QgsVectorLayer* mInputLayer;
+
+    QString mOutputFile;
+    QString mOutputFormat;
+
+    int mRadiusField;
+    int mWeightField;
+    double mRadius;
+    double mPixelSize;
+    QgsRectangle mBounds;
+
+    KernelShape mShape;
+    double mDecay;
+    OutputValues mOutputValues;
+
+    int mBufferSize;
+
+    GDALDatasetH mDatasetH;
+    GDALRasterBandH mRasterBandH;
+
+    //! Creates a new raster layer and initialises it to the no data value
+    bool createEmptyLayer( GDALDriverH driver, const QgsRectangle& bounds, int rows, int columns ) const;
+    int radiusSizeInPixels( double radius ) const;
+};
+
+
+#endif // QGSKDE_H