Fixes #3,#11,#12,#15 resolves #2,#14,#23 couple of perfomance tweaks

heatmap/heatmap.c

@@ -3,6 +3,9 @@
 #include <math.h>
 #include <string.h>
 
+float constant = 50.0;
+float multiplier = 200.0;
+
 struct info
 {
 	float minX;
@@ -33,7 +36,7 @@ BOOL WINAPI DllMain(    HINSTANCE hinstDLL,  // handle to DLL module
 #endif
 
 //walk the list of points, get the boundary values    
-void getBounds(struct info *inf, float *points, unsigned int cPoints)
+void getBounds(struct info *inf, float *points, unsigned int cPoints, int weighted)
 {
     unsigned int i = 0;
 
@@ -43,8 +46,11 @@ void getBounds(struct info *inf, float *points, unsigned int cPoints)
     float maxX = points[i];
     float maxY = points[i+1];
 
+    int inc = 2;
+    if (weighted) inc = 3;
+
     //then iterate over the list and find the max/min values
-    for(i = 0; i < cPoints; i=i+2)
+    for(i = 0; i < cPoints; i=i+inc)
     {
         float x = points[i];
         float y = points[i+1];
@@ -86,7 +92,7 @@ struct point translate(struct info *inf, struct point pt)
     return pt;
 }
 
-unsigned char* calcDensity(struct info *inf, float *points, int cPoints)
+unsigned char* calcDensity(struct info *inf, float *points, int cPoints, int weighted)
 {
     int width = inf->width;
     int height = inf->height;
@@ -110,7 +116,10 @@ unsigned char* calcDensity(struct info *inf, float *points, int cPoints)
         pixels[i] = 0xff;
     }
 
-    for(i = 0; i < cPoints; i=i+2)
+    int inc = 2;
+    if (weighted) inc = 3;
+
+    for(i = 0; i < cPoints; i=i+inc)
     {
         pt.x = points[i];
         pt.y = points[i+1];
@@ -123,8 +132,17 @@ unsigned char* calcDensity(struct info *inf, float *points, int cPoints)
                 if (j < 0 || k < 0 || j >= width || k >= height) continue; 
 
                 dist = sqrt( (j-pt.x)*(j-pt.x) + (k-pt.y)*(k-pt.y) );
-
-                pixVal = (int)(200.0*(dist/radius)+50.0);
+
+                if(dist>radius) continue; // stop point contributing to pixels outside its radius
+
+                if (weighted)
+                {
+                  pixVal = (int)((multiplier*(dist/radius)+constant)/points[i+2]);
+                }
+                else
+                {
+                  pixVal = (int)(multiplier*(dist/radius)+constant);
+                }
                 if (pixVal > 255) pixVal = 255;
 
                 ndx = k*width + j;
@@ -147,7 +165,6 @@ unsigned char *colorize(struct info *inf, unsigned char* pixels_bw, int *scheme,
 {
     int width = inf->width;
     int height = inf->height;
-    int dotsize = inf->dotsize;
 
     int i = 0;
     int pix = 0;
@@ -190,7 +207,7 @@ unsigned char *tx(float *points,
                   unsigned char *pix_color, 
                   int opacity, 
                   int boundsOverride, 
-                  float minX, float minY, float maxX, float maxY)
+                  float minX, float minY, float maxX, float maxY, int weighted)
 {
     unsigned char *pixels_bw = NULL;
 
@@ -216,7 +233,7 @@ unsigned char *tx(float *points,
     }
     else
     {
-        getBounds(&inf, points, cPoints);
+        getBounds(&inf, points, cPoints, weighted);
     }
 
     #ifdef DEBUG
@@ -225,7 +242,7 @@ unsigned char *tx(float *points,
 
     //iterate through points, place a dot at each center point
     //and set pix value from 0 - 255 using multiply method for radius [dotsize].
-    pixels_bw = calcDensity(&inf, points, cPoints);
+    pixels_bw = calcDensity(&inf, points, cPoints, weighted);
 
     //using provided color scheme and opacity, update pixel value to RGBA values
     pix_color = colorize(&inf, pixels_bw, scheme, pix_color, opacity);

heatmap/heatmap.py

@@ -3,11 +3,16 @@
 import ctypes
 import platform
 import math
-
 import colorschemes
-
 from PIL import Image
 
+use_pyproj = False
+try:
+  import pyproj
+  use_pyproj = True
+except:
+  pass
+
 class Heatmap:
     """
     Create heatmaps from a list of 2D coordinates.
@@ -45,8 +50,6 @@ class Heatmap:
 </kml>"""
 
     def __init__(self, libpath=None):
-        self.minXY = ()
-        self.maxXY = ()
         self.img = None
         # if you're reading this, it's probably because this
         # hacktastic garbage failed.  sorry.  I deserve a jab or two via @jjguy.
@@ -75,7 +78,7 @@ def __init__(self, libpath=None):
         if not self._heatmap:
             raise Exception("Heatmap shared library not found in PYTHONPATH.")
 
-    def heatmap(self, points, dotsize=150, opacity=128, size=(1024, 1024), scheme="classic", area=None):
+    def heatmap(self, points, dotsize=150, opacity=128, size=(1024, 1024), scheme="classic", area=None, weighted=0, srcepsg=None, dstepsg='EPSG:3857'):
         """
         points  -> an iterable list of tuples, where the contents are the
                    x,y coordinates to plot. e.g., [(1, 1), (2, 2), (3, 3)]
@@ -90,11 +93,20 @@ def heatmap(self, points, dotsize=150, opacity=128, size=(1024, 1024), scheme="c
         area    -> Specify bounding coordinates of the output image. Tuple of
                    tuples: ((minX, minY), (maxX, maxY)).  If None or unspecified,
                    these values are calculated based on the input data.
+        weighted -> Is the data weighted
+        srcepsg    -> epsg code of the source, set to None to ignore. If outputting to KML for google earth client overlay and have WGS84 or World Equidistant Cylindrical coords leave as None to save processing.
+        dstepsg    -> epsg code of the destination, ignored if srcepsg is not set. Defaults to EPSG:3857 (Cylindrical Mercator). Due to linear interpolation in heatmap.c it only makes sense to use linear output projections. If outputting to KML for google earth client overlay use EPSG:4087 (World Equidistant Cylindrical).
         """
         self.dotsize = dotsize
         self.opacity = opacity
         self.size = size
         self.points = points
+        self.weighted = weighted
+        self.srcepsg = srcepsg
+        self.dstepsg = dstepsg
+
+        if self.srcepsg and not use_pyproj:
+          raise Exception('srcepsg entered but pyproj is not available')
 
         if area is not None:
             self.area = area
@@ -103,21 +115,28 @@ def heatmap(self, points, dotsize=150, opacity=128, size=(1024, 1024), scheme="c
             self.area = ((0, 0), (0, 0))
             self.override = 0
 
+        #convert area for heatmap.c if required
+        ((east, south), (west, north)) = self.area
+        if use_pyproj and self.srcepsg is not None and self.srcepsg != self.dstepsg:
+          source = pyproj.Proj(init=self.srcepsg)
+          dest = pyproj.Proj(init=self.dstepsg)
+          (east,south) = pyproj.transform(source,dest,east,south)
+          (west,north) = pyproj.transform(source,dest,west,north)
+
         if scheme not in self.schemes():
             tmp = "Unknown color scheme: %s.  Available schemes: %s" % (
                 scheme, self.schemes())
             raise Exception(tmp)
 
-        arrPoints = self._convertPoints(points)
+        arrPoints = self._convertPoints()
         arrScheme = self._convertScheme(scheme)
         arrFinalImage = self._allocOutputBuffer()
 
         ret = self._heatmap.tx(
-            arrPoints, len(points) * 2, size[0], size[1], dotsize,
+            arrPoints, len(arrPoints), size[0], size[1], dotsize,
             arrScheme, arrFinalImage, opacity, self.override,
-            ctypes.c_float(self.area[0][0]), ctypes.c_float(
-                self.area[0][1]),
-            ctypes.c_float(self.area[1][0]), ctypes.c_float(self.area[1][1]))
+            ctypes.c_float(east), ctypes.c_float(south),
+            ctypes.c_float(west), ctypes.c_float(north), weighted)
 
         if not ret:
             raise Exception("Unexpected error during processing.")
@@ -129,43 +148,52 @@ def heatmap(self, points, dotsize=150, opacity=128, size=(1024, 1024), scheme="c
     def _allocOutputBuffer(self):
         return (ctypes.c_ubyte * (self.size[0] * self.size[1] * 4))()
 
-    def _convertPoints(self, pts):
+    def _convertPoints(self):
         """ flatten the list of tuples, convert into ctypes array """
 
-        #TODO is there a better way to do this??
-        flat = []
-        for i, j in pts:
-            flat.append(i)
-            flat.append(j)
-        #build array of input points
-        arr_pts = (ctypes.c_float * (len(pts) * 2))(*flat)
+        if isinstance(self.points,tuple):
+          self.points = list(self.points)
+        if isinstance(self.points[0],tuple):
+          self.points = list(sum(self.points,()))
+        elif isinstance(self.points[0],list):
+          self.points = sum(self.points,[])
+
+        #convert if required, need to copy as may use points later for _range.
+        if use_pyproj and self.srcepsg is not None and self.srcepsg != self.dstepsg:
+          converted =list(self.points)
+          source = pyproj.Proj(init=self.srcepsg)
+          dest = pyproj.Proj(init=self.dstepsg) 
+          #nicer way? map/lambda will retun 2/3 tuple so need to flatten again
+          inc = 3 if self.weighted else 2
+          for i in range(0, len(self.points), inc):
+            (x,y) = pyproj.transform(source,dest,self.points[i],self.points[i+1])
+            converted[i] = x
+            converted[i+1] = y
+          arr_pts = (ctypes.c_float * (len(converted))) (*converted)
+        else:
+          arr_pts = (ctypes.c_float * (len(self.points))) (*self.points)
         return arr_pts
 
     def _convertScheme(self, scheme):
         """ flatten the list of RGB tuples, convert into ctypes array """
 
-        #TODO is there a better way to do this??
-        flat = []
-        for r, g, b in colorschemes.schemes[scheme]:
-            flat.append(r)
-            flat.append(g)
-            flat.append(b)
-        arr_cs = (
-            ctypes.c_int * (len(colorschemes.schemes[scheme]) * 3))(*flat)
+        flat = list(sum(colorschemes.schemes[scheme],()))
+        arr_cs = (ctypes.c_int * (len(flat)))(*flat)
         return arr_cs
 
-    def _ranges(self, points):
+    def _ranges(self):
         """ walks the list of points and finds the
         max/min x & y values in the set """
-        minX = points[0][0]
-        minY = points[0][1]
+        minX = self.points[0]
+        minY = self.points[1]
         maxX = minX
         maxY = minY
-        for x, y in points:
-            minX = min(x, minX)
-            minY = min(y, minY)
-            maxX = max(x, maxX)
-            maxY = max(y, maxY)
+        inc = 3 if self.weighted else 2
+        for i in range(0,len(self.points),inc):
+            minX = min(self.points[i], minX)
+            minY = min(self.points[i+1], minY)
+            maxX = max(self.points[i], maxX)
+            maxY = max(self.points[i+1], maxY)
 
         return ((minX, minY), (maxX, maxY))
 
@@ -184,9 +212,16 @@ def saveKML(self, kmlFile):
         self.img.save(tilePath)
 
         if self.override:
-            ((east, south), (west, north)) = self.area
+            ((west, south), (east, north)) = self.area
         else:
-            ((east, south), (west, north)) = self._ranges(self.points)
+            ((west, south), (east, north)) = self._ranges()
+
+        #convert overlay BBOX if required
+        if use_pyproj and self.srcepsg is not None and self.srcepsg != 'EPSG:4326':
+          source = pyproj.Proj(init=self.srcepsg)
+          dest = pyproj.Proj(init='EPSG:4326')
+          (east,south) = pyproj.transform(source,dest,east,south)
+          (west,north) = pyproj.transform(source,dest,west,north)
 
         bytes = self.KML % (tilePath, north, south, east, west)
         file(kmlFile, "w").write(bytes)

requirements.txt

@@ -1 +1,2 @@
 Pillow>=2.1.0
+pyproj

setup.py

@@ -1,9 +1,22 @@
 import os
 import glob
 
-from distutils.core import setup, Extension
-from distutils.command.install import install
-from distutils.command.build_ext import build_ext
+#here use a flag so don't automatically use setuptools if available, hard to test otherwise
+with_setuptools = False
+if 'USE_SETUPTOOLS' in os.environ or 'pip' in __file__ or 'easy_install' in __file__:
+  try:
+    from setuptools.command.install import install
+    from setuptools import setup
+    from setuptools import Extension
+    from setuptools.command.build_ext import build_ext
+    with_setuptools = True
+  except ImportError:
+    pass
+if not with_setuptools:
+    from distutils.command.install import install
+    from distutils.core import setup
+    from distutils.core import Extension
+    from distutils.command.build_ext import build_ext
 
 # sorry for this, welcome feedback on the "right" way.
 # shipping pre-compiled bainries on windows, have
@@ -16,7 +29,6 @@ def run(self):
             return
         build_ext.run(self)
 
-
 class post_install(install):
     def run(self):
         install.run(self)
@@ -32,19 +44,30 @@ def run(self):
 
 cHeatmap = Extension('cHeatmap', sources=['heatmap/heatmap.c', ])
 
-setup(name='heatmap',
-      version="2.2.1",
-      description='Module to create heatmaps',
-      author='Jeffrey J. Guy',
-      author_email='[email protected]',
-      url='http://jjguy.com/heatmap/',
-      license='MIT License',
-      packages=['heatmap', ],
-      py_modules=['heatmap.colorschemes', ],
-      ext_modules=[cHeatmap, ],
-      cmdclass={'install': post_install,
-                'build_ext': mybuild},
-      requires=["Pillow", ],
-      test_suite="test",
-      tests_require=["Pillow", ],
-      )
+#separate calls to remove errors
+basekw = {
+      'name' : 'heatmap',
+      'version' : "2.2.1",
+      'description' : 'Module to create heatmaps',
+      'author' : 'Jeffrey J. Guy',
+      'author_email' : '[email protected]',
+      'url' : 'http://jjguy.com/heatmap/',
+      'license' : 'MIT License',
+      'packages' : ['heatmap', ],
+      'py_modules' : ['heatmap.colorschemes', ],
+      'ext_modules' : [cHeatmap, ],
+      'cmdclass' : {'install': post_install,
+                'build_ext': mybuild}
+      }
+setuptoolskw = {      
+      'install_requires' : ['Pillow'],
+      'extras_require' : {'proj' : 'pyproj'},
+      'test_suite' : "test",
+      'tests_require' : ['pyproj']
+      }
+distutilskw = {
+      'requires' : ["Pillow"]
+      }
+
+basekw.update(setuptoolskw) if with_setuptools else basekw.update(distutilskw)
+setup(**basekw)