Implement RFC108 heatmap support (#4857)

also adds:
 - rfc86 scaletoken substitutions for PROCESSING entries
 - proj fastpaths for 3857->4326 reprojections

Author: tbonfort

CMakeLists.txt

@@ -245,7 +245,7 @@ mapogcfiltercommon.c maprendering.c mapwcs20.c mapogcsld.c
 mapresample.c mapwfs.c mapgdal.c mapogcsos.c mapscale.c mapwfs11.c mapwfs20.c
 mapgeomtransform.c mapogroutput.c mapsde.c mapwfslayer.c mapagg.cpp mapkml.cpp
 mapgeomutil.cpp mapkmlrenderer.cpp fontcache.c textlayout.c maputfgrid.cpp
-mapogr.cpp mapcontour.c mapsmoothing.c mapv8.cpp ${REGEX_SOURCES})
+mapogr.cpp mapcontour.c mapsmoothing.c mapv8.cpp ${REGEX_SOURCES} kerneldensity.c)
 
 if(BUILD_DYNAMIC)
   add_library(mapserver SHARED ${mapserver_SOURCES} ${agg_SOURCES} ${v8_SOURCES})

HISTORY.TXT

@@ -15,6 +15,10 @@ For a complete change history, please see the Git log comments.
 7.0 release (TBD)
 -----------------
 
+- Apply RFC86 scaletoken substitutions to layer->PROCESSING entries
+
+- RFC108 Heatmap / Kernel-Density Layers
+
 - RFC104 Bitmap Label removal, replaced with inlined truetype font
 
 - RFC103 Layer Level Character Encoding

kerneldensity.c

@@ -0,0 +1,345 @@
+/******************************************************************************
+ *
+ * Project:  MapServer
+ * Purpose:  KernelDensity layer implementation and related functions.
+ * Author:   Thomas Bonfort and the MapServer team.
+ *
+ ******************************************************************************
+ * Copyright (c) 2014 Regents of the University of Minnesota.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies of this Software or works derived from this Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#include "mapserver.h"
+#include <float.h>
+#ifdef USE_GDAL
+
+#include "gdal.h"
+
+
+void gaussian_blur(float *values, int width, int height, int radius) {
+  float *tmp = (float*)msSmallMalloc(width*height*sizeof(float));
+  int length = radius*2+1;
+  float *kernel = (float*)msSmallMalloc(length*sizeof(float));
+  float sigma=radius/3.0;
+	float a=1.0/ sqrt(2.0*M_PI*sigma*sigma);
+	float den=2.0*sigma*sigma;
+	int i,x,y;
+
+	for (i=0; i<length; i++) {
+	  float x=i - radius;
+	  float v=a * exp(-(x*x) / den);
+	  kernel[i]=v;
+	}
+	memset(tmp,0,width*height*sizeof(float));
+
+	for(y=0; y<height; y++) {
+		float* src_row=values + width*y;
+		float* dst_row=tmp + width*y;
+
+		for(x=radius; x<width-radius; x++) {
+			float accum=0;
+			for(i=0; i<length; i++) {
+				accum+=src_row[x+i-radius] * kernel[i];
+			}
+			dst_row[x]=accum;
+		}
+	}
+
+	for(x=0; x<width; x++) {
+		float* src_col=tmp+x;
+		float* dst_col=values+x;
+
+		for(y=radius; y<height-radius; y++) {
+			float accum=0;
+			for (i=0; i<length; i++) {
+				accum+=src_col[width*(y+i-radius)] * kernel[i];
+			}
+			dst_col[y*width]=accum;
+		}
+	}
+  free(tmp);
+  free(kernel);
+}
+
+
+int msComputeKernelDensityDataset(mapObj *map, imageObj *image, layerObj *kerneldensity_layer, void **hDSvoid, void **cleanup_ptr) {
+
+  int status,layer_idx, i,j, nclasses=0, have_sample=0;
+  rectObj searchrect;
+  shapeObj shape;
+  layerObj *layer;
+  float *values;
+  int radius = 10, im_width = image->width, im_height = image->height;
+  int expand_searchrect=1;
+  float normalization_scale=0.0;
+  double invcellsize = 1.0 / map->cellsize, georadius=0;
+  float valmax=FLT_MIN, valmin=FLT_MAX;
+  unsigned char *iValues;
+  GDALDatasetH hDS;
+  const char *pszProcessing;
+  int *classgroup = NULL;
+  
+  assert(kerneldensity_layer->connectiontype == MS_KERNELDENSITY);
+  *cleanup_ptr = NULL;
+
+  if(!kerneldensity_layer->connection || !*kerneldensity_layer->connection) {
+    msSetError(MS_MISCERR, "msComputeKernelDensityDataset()", "KernelDensity layer has no CONNECTION defined");
+    return MS_FAILURE;
+  }
+  pszProcessing = msLayerGetProcessingKey( kerneldensity_layer, "KERNELDENSITY_RADIUS" );
+  if(pszProcessing)
+    radius = atoi(pszProcessing);
+  else
+    radius = 10;
+
+  pszProcessing = msLayerGetProcessingKey( kerneldensity_layer, "KERNELDENSITY_COMPUTE_BORDERS" );
+  if(pszProcessing && strcasecmp(pszProcessing,"OFF"))
+    expand_searchrect = 1;
+  else
+    expand_searchrect = 0;
+
+  pszProcessing = msLayerGetProcessingKey( kerneldensity_layer, "KERNELDENSITY_NORMALIZATION" );
+  if(!pszProcessing || !strcasecmp(pszProcessing,"AUTO"))
+    normalization_scale = 0.0;
+  else {
+    normalization_scale = atof(pszProcessing);
+    if(normalization_scale != 0) {
+      normalization_scale = 1.0 / normalization_scale;
+    } else {
+      normalization_scale = 1.0;
+    }
+  }
+
+  layer_idx = msGetLayerIndex(map,kerneldensity_layer->connection);
+  if(layer_idx == -1) {
+    int nLayers, *aLayers;
+    aLayers = msGetLayersIndexByGroup(map, kerneldensity_layer->connection, &nLayers);
+    if(!aLayers || !nLayers) {
+      msSetError(MS_MISCERR, "KernelDensity layer (%s) references unknown layer (%s)", "msComputeKernelDensityDataset()",
+                 kerneldensity_layer->name,kerneldensity_layer->connection);
+      return (MS_FAILURE);
+    }
+    for(i=0; i<nLayers; i++) {
+      layer_idx = aLayers[i];
+      layer = GET_LAYER(map, layer_idx);
+      if(msScaleInBounds(map->scaledenom, layer->minscaledenom, layer->maxscaledenom))
+        break;
+    }
+    free(aLayers);
+    if(i == nLayers) {
+      msSetError(MS_MISCERR, "KernelDensity layer (%s) references no layer for current scale", "msComputeKernelDensityDataset()",
+                 kerneldensity_layer->name);
+      return (MS_FAILURE);
+    }
+  } else {
+    layer = GET_LAYER(map, layer_idx);
+  }
+  /* open the linked layer */
+  status = msLayerOpen(layer);
+  if(status != MS_SUCCESS) return MS_FAILURE;
+
+  status = msLayerWhichItems(layer, MS_FALSE, NULL);
+  if(status != MS_SUCCESS) {
+    msLayerClose(layer);
+    return MS_FAILURE;
+  }
+
+  /* identify target shapes */
+  if(layer->transform == MS_TRUE) {
+    searchrect = map->extent;
+    if(expand_searchrect) {
+      georadius = radius * map->cellsize;
+      searchrect.minx -= georadius;
+      searchrect.miny -= georadius;
+      searchrect.maxx += georadius;
+      searchrect.maxy += georadius;
+      im_width += 2 * radius;
+      im_height += 2 * radius;
+    }
+  }
+  else {
+    searchrect.minx = searchrect.miny = 0;
+    searchrect.maxx = map->width-1;
+    searchrect.maxy = map->height-1;
+  }
+
+#ifdef USE_PROJ
+  layer->project = msProjectionsDiffer(&(layer->projection), &(map->projection));
+  if(layer->project)
+    msProjectRect(&map->projection, &layer->projection, &searchrect); /* project the searchrect to source coords */
+#endif
+
+  status = msLayerWhichShapes(layer, searchrect, MS_FALSE);
+  if(status == MS_DONE) { /* no overlap */
+    msLayerClose(layer);
+    return MS_SUCCESS;
+  } else if(status != MS_SUCCESS) {
+    msLayerClose(layer);
+    return MS_FAILURE;
+  }
+  values = (float*) msSmallCalloc(im_width * im_height, sizeof(float));
+  
+  if(layer->classgroup && layer->numclasses > 0)
+    classgroup = msAllocateValidClassGroups(layer, &nclasses);
+
+  msInitShape(&shape);
+  while((status = msLayerNextShape(layer, &shape)) == MS_SUCCESS) {
+    int l,p,s,c;
+    double weight = 1.0;
+#ifdef USE_PROJ
+    if(layer->project)
+      msProjectShape(&layer->projection, &map->projection, &shape);
+#endif
+
+    /* the weight for the sample is set to 1.0 by default. If the
+     * layer has some classes defined, we will read the weight from
+     * the class->style->size (which can be binded to an attribute)
+     */
+    if(layer->numclasses > 0) {
+      c = msShapeGetClass(layer, map, &shape, classgroup, nclasses);
+      if((c == -1) || (layer->class[c]->status == MS_OFF)) {
+        goto nextshape; /* no class matched, skip */
+      }
+      for (s = 0; s < layer->class[c]->numstyles; s++) {
+        if (msScaleInBounds(map->scaledenom,
+            layer->class[c]->styles[s]->minscaledenom,
+            layer->class[c]->styles[s]->maxscaledenom)) {
+          if(layer->class[c]->styles[s]->bindings[MS_STYLE_BINDING_SIZE].index != -1) {
+            weight = atof(shape.values[layer->class[c]->styles[s]->bindings[MS_STYLE_BINDING_SIZE].index]);
+          } else {
+            weight = layer->class[c]->styles[s]->size;
+          }
+          break;
+        }
+      }
+      if(s == layer->class[c]->numstyles) {
+        /* no style in scale bounds */
+        goto nextshape;
+      }
+    }
+    for(l=0; l<shape.numlines; l++) {
+      for(p=0; p<shape.line[l].numpoints; p++) {
+        int x = MS_MAP2IMAGE_XCELL_IC(shape.line[l].point[p].x, map->extent.minx - georadius, invcellsize);
+        int y = MS_MAP2IMAGE_YCELL_IC(shape.line[l].point[p].y, map->extent.maxy + georadius, invcellsize);
+        if(x>=0 && y>=0 && x<im_width && y<im_height) {
+          float *value = values + y * im_width + x;
+          (*value) += weight;
+          have_sample = 1;
+        }
+      }
+    }
+
+nextshape:
+    msFreeShape(&shape);
+  }
+  msLayerClose(layer);
+  if(status == MS_DONE) {
+    status = MS_SUCCESS;
+  } else {
+    status = MS_FAILURE;
+  }
+
+  if(have_sample) { /* no use applying the filtering kernel if we have no samples */
+    gaussian_blur(values,im_width, im_height, radius);
+
+    if(normalization_scale == 0.0) {   /* auto normalization */
+      for (j=radius; j<im_height-radius; j++) {
+        for (i=radius; i<im_width-radius; i++) {
+          float val = values[j*im_width + i];
+          if(val >0 && val>valmax) {
+            valmax = val;
+          }
+          if(val>0 && val<valmin) {
+            valmin = val;
+          }
+        }
+      }
+    } else {
+      valmin = 0;
+      valmax = normalization_scale;
+    }
+  }
+
+
+  if(have_sample && expand_searchrect) {
+    iValues = msSmallMalloc(image->width*image->height*sizeof(unsigned char));
+    for (j=0; j<image->height; j++) {
+      for (i=0; i<image->width; i++) {
+        float norm=(values[(j+radius)*im_width + i + radius] - valmin) / valmax;
+        int v=255 * norm;
+        if (v<0) v=0;
+        else if (v>255) v=255;
+        iValues[j*image->width + i] = v;
+      }
+    }
+  } else {
+    iValues = msSmallCalloc(1,image->width*image->height*sizeof(unsigned char));
+    if(have_sample) {
+       for (j=radius; j<image->height-radius; j++) {
+         for (i=radius; i<image->width-radius; i++) {
+           float norm=(values[j*im_width + i] - valmin) / valmax;
+           int v=255 * norm;
+           if (v<0) v=0;
+           else if (v>255) v=255;
+           iValues[j*image->width + i]=v;
+         }
+       }
+    }
+  }
+
+  free(values);
+  
+  {
+    char ds_string [1024];
+    double adfGeoTransform[6];
+    snprintf(ds_string,1024,"MEM:::DATAPOINTER=%p,PIXELS=%u,LINES=%u,BANDS=1,DATATYPE=Byte,PIXELOFFSET=1,LINEOFFSET=%u",
+        iValues,image->width,image->height,image->width);
+    hDS = GDALOpenShared( ds_string, GA_ReadOnly );
+    if(hDS==NULL) {
+      msSetError(MS_MISCERR,"msComputeKernelDensityDataset()","failed to create in-memory gdal dataset for interpolated data");
+      status = MS_FAILURE;
+      free(iValues);
+    }
+    adfGeoTransform[0] = map->extent.minx - map->cellsize * 0.5; /* top left x */
+    adfGeoTransform[1] = map->cellsize;/* w-e pixel resolution */
+    adfGeoTransform[2] = 0; /* 0 */
+    adfGeoTransform[3] = map->extent.maxy + map->cellsize * 0.5;/* top left y */
+    adfGeoTransform[4] =0; /* 0 */
+    adfGeoTransform[5] = -map->cellsize;/* n-s pixel resolution (negative value) */
+    GDALSetGeoTransform(hDS,adfGeoTransform);
+    *hDSvoid = hDS;
+    *cleanup_ptr = (void*)iValues;
+  }
+  return status;
+}
+#else
+
+
+int msComputeKernelDensityDataset(mapObj *map, imageObj *image, layerObj *layer, void **hDSvoid, void **cleanup_ptr) {
+    msSetError(MS_MISCERR,"msComputeKernelDensityDataset()", "KernelDensity layers require GDAL support, however GDAL support is not compiled in this build");
+    return MS_FAILURE;
+}
+
+#endif
+
+int msCleanupKernelDensityDataset(mapObj *map, imageObj *image, layerObj *layer, void *cleanup_ptr) {
+  free(cleanup_ptr);
+  return MS_SUCCESS;
+}