kerneldensity.c

/******************************************************************************
 *
 * 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>

#include "gdal.h"

static
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);
}

void msKernelDensityProcessing(layerObj *layer, interpolationProcessingParams *interpParams) {
    const char *interpParamsProcessing = msLayerGetProcessingKey( layer, "KERNELDENSITY_RADIUS" );
    if(interpParamsProcessing) {
        interpParams->radius = atoi(interpParamsProcessing);
    } else {
        interpParams->radius = 10;
    }

    interpParamsProcessing = msLayerGetProcessingKey( layer, "KERNELDENSITY_COMPUTE_BORDERS" );
    if(interpParamsProcessing && strcasecmp(interpParamsProcessing,"OFF")) {
        interpParams->expand_searchrect = 1;
    } else {
        interpParams->expand_searchrect = 0;
    }

    interpParamsProcessing = msLayerGetProcessingKey( layer, "KERNELDENSITY_NORMALIZATION" );
    if(!interpParamsProcessing || !strcasecmp(interpParamsProcessing,"AUTO")) {
        interpParams->normalization_scale = 0.0;
    } else {
        interpParams->normalization_scale = atof(interpParamsProcessing);
        if(interpParams->normalization_scale != 0) {
            interpParams->normalization_scale = 1.0 / interpParams->normalization_scale;
        } else {
            interpParams->normalization_scale = 1.0;
        }
    }
}

void msKernelDensity(imageObj *image, float *values, int width, int height, int npoints,
                   interpolationProcessingParams *interpParams, unsigned char *iValues) {
    int i,j;
    float valmax=FLT_MIN, valmin=FLT_MAX;
    int radius = interpParams->radius;
    float normalization_scale = interpParams->normalization_scale;
    int expand_searchrect = interpParams->expand_searchrect;

    gaussian_blur(values, width, height, radius);

    if(normalization_scale == 0.0) {   /* auto normalization */
        for (j = radius; j < height-radius; j++) {
            for (i = radius; i < width-radius; i++) {
                float val = values[j*width + i];
                if(val >0 && val>valmax) {
                    valmax = val;
                }
                if(val>0 && val<valmin) {
                    valmin = val;
                }
            }
        }
    } else {
        valmin = 0;
        valmax = normalization_scale;
    }

    if(expand_searchrect) {
        for (j=0; j < image->height; j++) {
            for (i=0; i < image->width; i++) {
                float norm = (values[(j+radius)*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 {
        if(npoints > 0) {
            for (j=radius; j < image->height-radius; j++) {
                for (i=radius; i < image->width-radius; i++) {
                    float norm=(values[j*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;
                }
            }
        }
    }
}