[FEATURE] Add triweight and Epanechnikov kernels. Make decay disabled…

… for all kernels expect Triangular kernel. Fix math for triangular kernels.
qgis · Mar 18, 2013 · 80c3981 · 80c3981
1 parent bf45cb2
commit 80c3981
Show file tree

Hide file tree

Showing 5 changed files with 99 additions and 25 deletions.
diff --git a/src/plugins/heatmap/heatmap.cpp b/src/plugins/heatmap/heatmap.cpp
@@ -111,7 +111,7 @@ void Heatmap::run()
     int columns = d.columns();
     int rows = d.rows();
     double cellsize = d.cellSizeX(); // or d.cellSizeY();  both have the same value
-    double myDecay = d.decayRatio();
+    mDecay = d.decayRatio();
     int kernelShape = d.kernelShape();
 
     // Start working on the input vector
@@ -262,7 +262,7 @@ void Heatmap::run()
           {
             continue;
           }
-          
+
           double pixelValue = weight * calculateKernelValue( distance, myBuffer, kernelShape );
 
           // clearing anamolies along the axes
@@ -332,41 +332,84 @@ int Heatmap::bufferSize( double radius, double cellsize )
     ++buffer;
   }
   return buffer;
-}    
-    
+}
+
 double Heatmap::calculateKernelValue( double distance, int bandwidth, int kernelShape )
 {
-  switch (kernelShape) {
+  switch ( kernelShape )
+  {
     case HeatmapGui::Triangular:
-      return ( 1 - ( distance / bandwidth ) );
-      
+      return triangularKernel( distance , bandwidth );
+
     case HeatmapGui::Uniform:
       return uniformKernel( distance, bandwidth );
-      
+
     case HeatmapGui::Quartic:
       return quarticKernel( distance, bandwidth );
+
+    case HeatmapGui::Triweight:
+      return triweightKernel( distance, bandwidth );
+
+    case HeatmapGui::Epanechnikov:
+      return epanechnikovKernel( distance, bandwidth );
   }
   return 0;
-  
+
 }
 
+/* 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 Heatmap::uniformKernel( double distance, int bandwidth )
 {
-  // Normalizing constant. Calculated by polar double integrating the kernel function
-  // with radius of 0 to bandwidth and equating area to 1.
-  double k = 2. / (M_PI * (double)bandwidth);
-
+  // Normalizing constant
+  double k = 2. / ( M_PI * ( double )bandwidth );
+
   // Derived from Wand and Jones (1995), p. 175
-  return k * ( 0.5 / (double)bandwidth);
+  return k * ( 0.5 / ( double )bandwidth );
 }
 
 double Heatmap::quarticKernel( double distance, int bandwidth )
 {
   // Normalizing constant
-  double k = 16. / (5. * M_PI * pow((double)bandwidth, 2));
-
+  double k = 16. / ( 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 );
+}
+
+double Heatmap::triweightKernel( double distance, int bandwidth )
+{
+  // 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 );
+}
+
+double Heatmap::epanechnikovKernel( double distance, int bandwidth )
+{
+  // Normalizing constant
+  double k = 8. / ( 3. * 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);
+  return k * ( 3. / 4. ) * ( 1. - pow( distance / ( double )bandwidth, 2 ) );
+}
+
+double Heatmap::triangularKernel( double distance, int bandwidth )
+{
+  // Normalizing constant. In this case it's calculated a little different
+  // due to the inclusion of the non-standard "decay" parameter
+
+  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 ) );
 }
 
 // Unload the plugin by cleaning up the GUI

diff --git a/src/plugins/heatmap/heatmap.h b/src/plugins/heatmap/heatmap.h
@@ -81,6 +81,8 @@ class Heatmap: public QObject, public QgisPlugin
     void help();
 
   private:
+    double mDecay;
+
     //! Worker to convert meters to map units
     double mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs );
     //! Worker to calculate buffer size in pixels
@@ -89,10 +91,15 @@ class Heatmap: public QObject, public QgisPlugin
     double calculateKernelValue( double distance, int bandwidth, int kernelShape );
     //! Uniform kernel function
     double uniformKernel( double distance, int bandwidth );
-    //! Quartic kernel function    
+    //! Quartic kernel function
     double quarticKernel( double distance, int bandwidth );
-
-
+    //! Triweight kernel function
+    double triweightKernel( double distance, int bandwidth );
+    //! Epanechnikov kernel function
+    double epanechnikovKernel( double distance, int bandwidth );
+    //! Triangular kernel function
+    double triangularKernel( double distance, int bandwidth );
+
     // MANDATORY PLUGIN PROPERTY DECLARATIONS  .....
 
     int mPluginType;

diff --git a/src/plugins/heatmap/heatmapgui.cpp b/src/plugins/heatmap/heatmapgui.cpp
@@ -225,6 +225,14 @@ void HeatmapGui::on_mBufferLineEdit_editingFinished()
 {
   updateBBox();
 }
+
+void HeatmapGui::on_kernelShapeCombo_currentIndexChanged( int index )
+{
+  Q_UNUSED( index );
+  // Only enable the decay edit if the kernel shape is set to triangular
+  mDecayLineEdit->setEnabled( index == HeatmapGui::Triangular );
+}
+
 /*
  *
  * Private Functions

diff --git a/src/plugins/heatmap/heatmapgui.h b/src/plugins/heatmap/heatmapgui.h
@@ -35,12 +35,14 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
       Meters,
       MapUnits
     };
-    
+
     enum kernelShape
     {
+      Quartic,
       Triangular,
       Uniform,
-      Quartic
+      Triweight,
+      Epanechnikov
     };
 
     /** Returns whether to apply weighted heat */
@@ -54,7 +56,7 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
 
     /** Return the radius Unit (meters/map units) */
     int radiusUnit();
-    
+
     /** Return the selected kernel shape */
     int kernelShape();
 
@@ -130,6 +132,7 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
     void on_mRadiusUnitCombo_currentIndexChanged( int index );
     void on_mInputVectorCombo_currentIndexChanged( int index );
     void on_mBufferLineEdit_editingFinished();
+    void on_kernelShapeCombo_currentIndexChanged( int index );
 };
 
 #endif
diff --git a/src/plugins/heatmap/heatmapguibase.ui b/src/plugins/heatmap/heatmapguibase.ui
@@ -227,6 +227,9 @@
         </item>
         <item row="3" column="1" colspan="2">
          <widget class="QLineEdit" name="mDecayLineEdit">
+          <property name="enabled">
+           <bool>false</bool>
+          </property>
           <property name="text">
            <string>0.0</string>
           </property>
@@ -241,6 +244,11 @@
         </item>
         <item row="0" column="1" colspan="2">
          <widget class="QComboBox" name="kernelShapeCombo">
+          <item>
+           <property name="text">
+            <string>Quartic (Biweight)</string>
+           </property>
+          </item>
           <item>
            <property name="text">
             <string>Triangular</string>
@@ -253,9 +261,14 @@
           </item>
           <item>
            <property name="text">
-            <string>Quartic</string>
+            <string>Triweight</string>
+           </property>
+          </item>
+          <item>
+           <property name="text">
+            <string>Epanechnikov</string>
            </property>
-          </item>          
+          </item>
          </widget>
         </item>
         <item row="0" column="0">