Rails Against The Machine

In this entry I'm going to show you how to easily create beautiful heatmaps for google maps using the technique of The definitive heat map and borrowing code from Johan Liesén The definitive heat map – with Java

So without further ado here is the sourcecode:

I'm using Bear Bibeault's Front Man™ as a lightweight server for my map tiles.

The command:

package commands;

import java.io.IOException;
import java.io.OutputStream;
import javax.imageio.ImageIO;
import javax.servlet.ServletException;
import model.DensityMap;
import org.bibeault.frontman.Command;
import org.bibeault.frontman.CommandContext;
public class GetTileCommand implements Command {

    public void execute(CommandContext context) throws ServletException, IOException {
        int xCoord=Integer.parseInt(context.getRequest().getParameter("x"));
        int yCoord=Integer.parseInt(context.getRequest().getParameter("y"));
        int zoom=Integer.parseInt(context.getRequest().getParameter("zoom"));
        
        context.getResponse().setContentType("image/png");
        OutputStream output = context.getResponse().getOutputStream();
        DensityMap h=DensityMap.getInstance();
        ImageIO.write(h.getTile(xCoord,yCoord,zoom), "png", output);
    }

}

Loading the data and plotting it

package model;
import DensityHeatMap.HeatmapImage;
import Jama.Matrix;
import java.awt.image.BufferedImage;
import database.FileHandling;
import java.awt.Point;
import services.sharedfunctions.Conversions;

public class DensityMap {

    private static DensityMap instance = null;
    private Matrix coordinates; //Training data points
    private Matrix values; //Training data values

    protected DensityMap() {
        Matrix[] test_data = FileHandling.readTrainingData("data.csv");
        coordinates = test_data[0];
        values = test_data[1];
    }

    public static DensityMap getInstance() {
        if (instance == null) {
            instance = new DensityMap();
        }
        return instance;
    }
    
    public BufferedImage getTile(int x_tile, int y_tile, int zoom) {
        HeatmapImage h=new HeatmapImage(255,255,zoom);
        for (int i = 0; i < coordinates.getColumnDimension(); i++) {
            double lng=coordinates.get(0,i);
            double lat=coordinates.get(1,i);
            int[] x=Conversions.lng_to_pixel(lng, zoom);
            int[] y=Conversions.lat_to_pixel(lat, zoom);
            
            //relative to coordinate frame of tile i.e. 0,0 is top left hand corner.
            int x_pixels=x[1]+(x[0]-x_tile)*255;
            int y_pixels=y[1]+(y[0]-y_tile)*255;

            //we need to plot points which are just off the tile for continuity.
            //idealy the cut off should depend on zoom level and basis function size
            if((Math.abs(x[0]-x_tile)<=1)&&(Math.abs(y[0]-y_tile)<=1))
            {
                //x_Pixel               
                Point p=new Point(x_pixels, y_pixels);
                //CAUTION
                //use for densities only dont try to use for interpolation of valued functions
                h.addDotImage(p, 0.25f);
            }
            
        }
        return h.getImage();
    }
}

LatLong mercator conversions

package services.sharedfunctions;

public class Conversions {

    static int width = 256;
    static int height = 256;
    private static double earthRadius = 6367500.0;

    public static double pixel_to_lat(int pixel, int tile, int zoom) {
        double percent_down = (pixel + (height * tile)) / (height * Math.pow(2, zoom));//work out what % of the way down we are

        double from_equator = 1f - 2f * percent_down;//but in standard mercator projection the equator is zero 

        double y_mercator = from_equator * Math.PI;//and units run from + PI to -PI       

        double latitude = 2 * Math.atan(Math.exp(y_mercator)) - Math.PI / 2;
        return 180 * latitude / Math.PI;//convert to degrees
    }

    //returns tile, pixels
    public static int[] lat_to_pixel(double lat, int zoom) {
        //convert from degrees to radians
        double latitude = lat * Math.PI / 180;
        //convert to percator units that run from + PI to -PI   
        double y_mercator = Math.log(Math.tan(latitude / 2 + Math.PI / 4));
        //Now to percent up or down from equator
        double percent_down = (1f - (y_mercator / Math.PI)) / 2;
        //pixels
        double pixels = percent_down * (height * Math.pow(2, zoom));
        int tile = (int) Math.floor(pixels / height);
        int pixel = (int) (pixels - tile * height);
        int[] tmp = {tile, pixel};
        return tmp;
    }

    public static double pixel_to_lng(int pixel, int tile, int zoom) {
        double circumference = width * Math.pow(2, zoom); //Number of tiles forming x axis

        double radius = circumference / (2 * Math.PI);
        double longitude = (pixel + width * tile) / radius;//in radians its just the ratio

        longitude = 180 * longitude / Math.PI;//convert to degrees

        return longitude - 180;//change center from date line to grenich

    }

    public static int[] lng_to_pixel(double lng, int zoom) {
        //convert from degrees to radians
        double longitude = Math.PI * (lng + 180) / 180;
        double circumference = width * Math.pow(2, zoom); //Number of tiles forming x axis
        double radius = circumference / (2 * Math.PI);
        double pixels = longitude * radius;
        int tile = (int) Math.floor(pixels / width);
        int pixel = (int) (pixels - tile * width);
        int[] tmp = {tile, pixel};
        return tmp;
    }

    public static double lat_lon_to_x_meters(double lat, double lon, double lon_origin) {
        return (Math.PI * earthRadius / 180.0) * (lon - lon_origin) * Math.cos(lat * (Math.PI / 180.0));
    }

    public static double lat_lon_to_y_meters(double lat, double lat_origin) {
        return (Math.PI * earthRadius / 180.0) * (lat - lat_origin);
    }
}

Generating the heat map
just taken from Johan Liesén and refactored except that we now scale the size of the dotImage with zoom level.

package DensityHeatMap;

import java.awt.*;
import java.awt.image.*;
import services.sharedfunctions.Conversions;


public class HeatmapImage {

    private BufferedImage dotImage;
    private BufferedImage monochromeImage;
    private BufferedImage heatmapImage;
    private BufferedImage colorImage;
    private LookupTable colorTable;
    private LookupOp colorOp;

    public HeatmapImage(int width, int height,int zoom) {
        //create the image used to represent the basis functions
        double delta_lng=Conversions.pixel_to_lng(1,0,zoom)-Conversions.pixel_to_lng(0,0,zoom);
        int size=(int) (0.01 / delta_lng);
        dotImage = createFadedCircleImage(size);
        
        //Create a 1 pixel by 64 reference image to define the colours
        //we could load this instead of generating it.
        colorImage = createEvenlyDistributedGradientImage(new Dimension(
                64, 1), Color.WHITE, Color.RED, Color.YELLOW, Color.GREEN.darker(), Color.CYAN, Color.BLUE, new Color(0, 0, 0x33));
        //create a colour lookup table from the reference image
        colorTable = createColorLookupTable(colorImage, .5f);
        colorOp = new LookupOp(colorTable, null);

        //Create a monochrome image to hold our density plot
        monochromeImage = createCompatibleTranslucentImage(width, height);

        Graphics g = monochromeImage.getGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, width, height);
    }

    public void addDotImage(Point p, float alpha) {
        int circleRadius = dotImage.getWidth() / 2;

        Graphics2D g = (Graphics2D) monochromeImage.getGraphics();

        g.setComposite(BlendComposite.Multiply.derive(alpha));
        g.drawImage(dotImage, null, p.x - circleRadius, p.y - circleRadius);
    }

    public BufferedImage getImage() {
        return heatmapImage = colorize(colorOp);
    }
    //--------------------------------------------------------------------------

    private BufferedImage colorize(LookupOp colorOp) {
        return colorOp.filter(monochromeImage, null);
    }

    private BufferedImage colorize(LookupTable colorTable) {
        return colorize(new LookupOp(colorTable, null));
    }

    private static BufferedImage createFadedCircleImage(int size) {
        BufferedImage im = createCompatibleTranslucentImage(size, size);
        float radius = size / 2f;

        RadialGradientPaint gradient = new RadialGradientPaint(
                radius, radius, radius, new float[]{0f, 1f}, new Color[]{
                    Color.BLACK, new Color(0xffffffff, true)
                });

        Graphics2D g = (Graphics2D) im.getGraphics();

        g.setPaint(gradient);
        g.fillRect(0, 0, size, size);

        g.dispose();

        return im;
    }

    private static BufferedImage createCompatibleTranslucentImage(int width,
            int height) {
        GraphicsEnvironment ge = GraphicsEnvironment.getLocalGraphicsEnvironment();
        GraphicsDevice gd = ge.getDefaultScreenDevice();
        GraphicsConfiguration gc = gd.getDefaultConfiguration();

        return gc.createCompatibleImage(
                width, height, Transparency.TRANSLUCENT);
    }

    /**
     * Creates the color lookup table from an image
     * 
     * @param im
     * @return
     */
    public static LookupTable createColorLookupTable(BufferedImage im,
            float alpha) {
        int tableSize = 256;
        Raster imageRaster = im.getData();
        double sampleStep = 1D * im.getWidth() / tableSize; // Sample pixels
        // evenly
        byte[][] colorTable = new byte[4][tableSize];
        int[] pixel = new int[1]; // Sample pixel
        Color c;

        for (int i = 0; i < tableSize; ++i) {
            imageRaster.getDataElements((int) (i * sampleStep), 0, pixel);

            c = new Color(pixel[0]);

            colorTable[0][i] = (byte) c.getRed();
            colorTable[1][i] = (byte) c.getGreen();
            colorTable[2][i] = (byte) c.getBlue();
            colorTable[3][i] = (byte) (alpha * 0xff);
        }

        LookupTable lookupTable = new ByteLookupTable(0, colorTable);

        return lookupTable;
    }

    public static BufferedImage createEvenlyDistributedGradientImage(
            Dimension size, Color... colors) {
        BufferedImage im = createCompatibleTranslucentImage(
                size.width, size.height);
        Graphics2D g = im.createGraphics();

        float[] fractions = new float[colors.length];
        float step = 1f / colors.length;

        for (int i = 0; i < colors.length; ++i) {
            fractions[i] = i * step;
        }

        LinearGradientPaint gradient = new LinearGradientPaint(
                0, 0, size.width, 1, fractions, colors,
                MultipleGradientPaint.CycleMethod.REPEAT);

        g.setPaint(gradient);
        g.fillRect(0, 0, size.width, size.height);

        g.dispose();

        return im;
    }
}

BlendComposite by Romain Guy

package DensityHeatMap;



import java.awt.Composite;
import java.awt.CompositeContext;
import java.awt.RenderingHints;
import java.awt.image.ColorModel;
import java.awt.image.DataBuffer;
import java.awt.image.Raster;
import java.awt.image.WritableRaster;

public final class BlendComposite implements Composite {
    public enum BlendingMode {
        NORMAL,
        AVERAGE,
        MULTIPLY,
        SCREEN,
        DARKEN,
        LIGHTEN,
        OVERLAY,
        HARD_LIGHT,
        SOFT_LIGHT,
        DIFFERENCE,
        NEGATION,
        EXCLUSION,
        COLOR_DODGE,
        INVERSE_COLOR_DODGE,
        SOFT_DODGE,
        COLOR_BURN,
        INVERSE_COLOR_BURN,
        SOFT_BURN,
        REFLECT,
        GLOW,
        FREEZE,
        HEAT,
        ADD,
        SUBTRACT,
        STAMP,
        RED,
        GREEN,
        BLUE,
        HUE,
        SATURATION,
        COLOR,
        LUMINOSITY
    }

    public static final BlendComposite Normal = new BlendComposite(BlendingMode.NORMAL);
    public static final BlendComposite Average = new BlendComposite(BlendingMode.AVERAGE);
    public static final BlendComposite Multiply = new BlendComposite(BlendingMode.MULTIPLY);
    public static final BlendComposite Screen = new BlendComposite(BlendingMode.SCREEN);
    public static final BlendComposite Darken = new BlendComposite(BlendingMode.DARKEN);
    public static final BlendComposite Lighten = new BlendComposite(BlendingMode.LIGHTEN);
    public static final BlendComposite Overlay = new BlendComposite(BlendingMode.OVERLAY);
    public static final BlendComposite HardLight = new BlendComposite(BlendingMode.HARD_LIGHT);
    public static final BlendComposite SoftLight = new BlendComposite(BlendingMode.SOFT_LIGHT);
    public static final BlendComposite Difference = new BlendComposite(BlendingMode.DIFFERENCE);
    public static final BlendComposite Negation = new BlendComposite(BlendingMode.NEGATION);
    public static final BlendComposite Exclusion = new BlendComposite(BlendingMode.EXCLUSION);
    public static final BlendComposite ColorDodge = new BlendComposite(BlendingMode.COLOR_DODGE);
    public static final BlendComposite InverseColorDodge = new BlendComposite(BlendingMode.INVERSE_COLOR_DODGE);
    public static final BlendComposite SoftDodge = new BlendComposite(BlendingMode.SOFT_DODGE);
    public static final BlendComposite ColorBurn = new BlendComposite(BlendingMode.COLOR_BURN);
    public static final BlendComposite InverseColorBurn = new BlendComposite(BlendingMode.INVERSE_COLOR_BURN);
    public static final BlendComposite SoftBurn = new BlendComposite(BlendingMode.SOFT_BURN);
    public static final BlendComposite Reflect = new BlendComposite(BlendingMode.REFLECT);
    public static final BlendComposite Glow = new BlendComposite(BlendingMode.GLOW);
    public static final BlendComposite Freeze = new BlendComposite(BlendingMode.FREEZE);
    public static final BlendComposite Heat = new BlendComposite(BlendingMode.HEAT);
    public static final BlendComposite Add = new BlendComposite(BlendingMode.ADD);
    public static final BlendComposite Subtract = new BlendComposite(BlendingMode.SUBTRACT);
    public static final BlendComposite Stamp = new BlendComposite(BlendingMode.STAMP);
    public static final BlendComposite Red = new BlendComposite(BlendingMode.RED);
    public static final BlendComposite Green = new BlendComposite(BlendingMode.GREEN);
    public static final BlendComposite Blue = new BlendComposite(BlendingMode.BLUE);
    public static final BlendComposite Hue = new BlendComposite(BlendingMode.HUE);
    public static final BlendComposite Saturation = new BlendComposite(BlendingMode.SATURATION);
    public static final BlendComposite Color = new BlendComposite(BlendingMode.COLOR);
    public static final BlendComposite Luminosity = new BlendComposite(BlendingMode.LUMINOSITY);

    private float alpha;
    private BlendingMode mode;

    private BlendComposite(BlendingMode mode) {
        this(mode, 1.0f);
    }

    private BlendComposite(BlendingMode mode, float alpha) {
        this.mode = mode;
        setAlpha(alpha);
    }

    public static BlendComposite getInstance(BlendingMode mode) {
        return new BlendComposite(mode);
    }

    public static BlendComposite getInstance(BlendingMode mode, float alpha) {
        return new BlendComposite(mode, alpha);
    }

    public BlendComposite derive(BlendingMode mode) {
        return this.mode == mode ? this : new BlendComposite(mode, getAlpha());
    }

    public BlendComposite derive(float alpha) {
        return this.alpha == alpha ? this : new BlendComposite(getMode(), alpha);
    }

    public float getAlpha() {
        return alpha;
    }

    public BlendingMode getMode() {
        return mode;
    }

    private void setAlpha(float alpha) {
        if (alpha < 0.0f || alpha > 1.0f) {
            throw new IllegalArgumentException(
                    "alpha must be comprised between 0.0f and 1.0f");
        }

        this.alpha = alpha;
    }

    @Override
    public int hashCode() {
        return Float.floatToIntBits(alpha) * 31 + mode.ordinal();
    }

    @Override
    public boolean equals(Object obj) {
        if (!(obj instanceof BlendComposite)) {
            return false;
        }

        BlendComposite bc = (BlendComposite) obj;

        if (mode != bc.mode) {
            return false;
        }

        return alpha == bc.alpha;
    }

    public CompositeContext createContext(ColorModel srcColorModel,
                                          ColorModel dstColorModel,
                                          RenderingHints hints) {
        return new BlendingContext(this);
    }

    private static final class BlendingContext implements CompositeContext {
        private final Blender blender;
        private final BlendComposite composite;

        private BlendingContext(BlendComposite composite) {
            this.composite = composite;
            this.blender = Blender.getBlenderFor(composite);
        }

        public void dispose() {
        }

        public void compose(Raster src, Raster dstIn, WritableRaster dstOut) {
            if (src.getSampleModel().getDataType() != DataBuffer.TYPE_INT ||
                dstIn.getSampleModel().getDataType() != DataBuffer.TYPE_INT ||
                dstOut.getSampleModel().getDataType() != DataBuffer.TYPE_INT) {
                throw new IllegalStateException(
                        "Source and destination must store pixels as INT.");
            }

            int width = Math.min(src.getWidth(), dstIn.getWidth());
            int height = Math.min(src.getHeight(), dstIn.getHeight());

            float alpha = composite.getAlpha();

            int[] srcPixel = new int[4];
            int[] dstPixel = new int[4];
            int[] srcPixels = new int[width];
            int[] dstPixels = new int[width];

            for (int y = 0; y < height; y++) {
                src.getDataElements(0, y, width, 1, srcPixels);
                dstIn.getDataElements(0, y, width, 1, dstPixels);
                for (int x = 0; x < width; x++) {
                    // pixels are stored as INT_ARGB
                    // our arrays are [R, G, B, A]
                    int pixel = srcPixels[x];
                    srcPixel[0] = (pixel >> 16) & 0xFF;
                    srcPixel[1] = (pixel >>  8) & 0xFF;
                    srcPixel[2] = (pixel      ) & 0xFF;
                    srcPixel[3] = (pixel >> 24) & 0xFF;

                    pixel = dstPixels[x];
                    dstPixel[0] = (pixel >> 16) & 0xFF;
                    dstPixel[1] = (pixel >>  8) & 0xFF;
                    dstPixel[2] = (pixel      ) & 0xFF;
                    dstPixel[3] = (pixel >> 24) & 0xFF;

                    int[] result = blender.blend(srcPixel, dstPixel);

                    // mixes the result with the opacity
                    dstPixels[x] = ((int) (dstPixel[3] + (result[3] - dstPixel[3]) * alpha) & 0xFF) << 24 |
                                   ((int) (dstPixel[0] + (result[0] - dstPixel[0]) * alpha) & 0xFF) << 16 |
                                   ((int) (dstPixel[1] + (result[1] - dstPixel[1]) * alpha) & 0xFF) <<  8 |
                                    (int) (dstPixel[2] + (result[2] - dstPixel[2]) * alpha) & 0xFF;
                }
                dstOut.setDataElements(0, y, width, 1, dstPixels);
            }
        }
    }

    private static abstract class Blender {
        public abstract int[] blend(int[] src, int[] dst);

        private static void RGBtoHSL(int r, int g, int b, float[] hsl) {
            float var_R = (r / 255f);
            float var_G = (g / 255f);
            float var_B = (b / 255f);

            float var_Min;
            float var_Max;
            float del_Max;

            if (var_R > var_G) {
                var_Min = var_G;
                var_Max = var_R;
            } else {
                var_Min = var_R;
                var_Max = var_G;
            }
            if (var_B > var_Max) {
                var_Max = var_B;
            }
            if (var_B < var_Min) {
                var_Min = var_B;
            }

            del_Max = var_Max - var_Min;

            float H, S, L;
            L = (var_Max + var_Min) / 2f;

            if (del_Max - 0.01f <= 0.0f) {
                H = 0;
                S = 0;
            } else {
                if (L < 0.5f) {
                    S = del_Max / (var_Max + var_Min);
                } else {
                    S = del_Max / (2 - var_Max - var_Min);
                }

                float del_R = (((var_Max - var_R) / 6f) + (del_Max / 2f)) / del_Max;
                float del_G = (((var_Max - var_G) / 6f) + (del_Max / 2f)) / del_Max;
                float del_B = (((var_Max - var_B) / 6f) + (del_Max / 2f)) / del_Max;

                if (var_R == var_Max) {
                    H = del_B - del_G;
                } else if (var_G == var_Max) {
                    H = (1 / 3f) + del_R - del_B;
                } else {
                    H = (2 / 3f) + del_G - del_R;
                }
                if (H < 0) {
                    H += 1;
                }
                if (H > 1) {
                    H -= 1;
                }
            }

            hsl[0] = H;
            hsl[1] = S;
            hsl[2] = L;
        }

        private static void HSLtoRGB(float h, float s, float l, int[] rgb) {
            int R, G, B;

            if (s - 0.01f <= 0.0f) {
                R = (int) (l * 255.0f);
                G = (int) (l * 255.0f);
                B = (int) (l * 255.0f);
            } else {
                float var_1, var_2;
                if (l < 0.5f) {
                    var_2 = l * (1 + s);
                } else {
                    var_2 = (l + s) - (s * l);
                }
                var_1 = 2 * l - var_2;

                R = (int) (255.0f * hue2RGB(var_1, var_2, h + (1.0f / 3.0f)));
                G = (int) (255.0f * hue2RGB(var_1, var_2, h));
                B = (int) (255.0f * hue2RGB(var_1, var_2, h - (1.0f / 3.0f)));
            }

            rgb[0] = R;
            rgb[1] = G;
            rgb[2] = B;
        }

        private static float hue2RGB(float v1, float v2, float vH) {
            if (vH < 0.0f) {
                vH += 1.0f;
            }
            if (vH > 1.0f) {
                vH -= 1.0f;
            }
            if ((6.0f * vH) < 1.0f) {
                return (v1 + (v2 - v1) * 6.0f * vH);
            }
            if ((2.0f * vH) < 1.0f) {
                return (v2);
            }
            if ((3.0f * vH) < 2.0f) {
                return (v1 + (v2 - v1) * ((2.0f / 3.0f) - vH) * 6.0f);
            }
            return (v1);
        }

        public static Blender getBlenderFor(BlendComposite composite) {
            switch (composite.getMode()) {
                case NORMAL:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return src;
                        }
                    };
                case ADD:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.min(255, src[0] + dst[0]),
                                Math.min(255, src[1] + dst[1]),
                                Math.min(255, src[2] + dst[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case AVERAGE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                (src[0] + dst[0]) >> 1,
                                (src[1] + dst[1]) >> 1,
                                (src[2] + dst[2]) >> 1,
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case BLUE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0],
                                src[1],
                                dst[2],
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case COLOR:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            float[] srcHSL = new float[3];
                            RGBtoHSL(src[0], src[1], src[2], srcHSL);
                            float[] dstHSL = new float[3];
                            RGBtoHSL(dst[0], dst[1], dst[2], dstHSL);

                            int[] result = new int[4];
                            HSLtoRGB(srcHSL[0], srcHSL[1], dstHSL[2], result);
                            result[3] = Math.min(255, src[3] + dst[3]);

                            return result;
                        }
                    };
                case COLOR_BURN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - dst[0]) << 8) / src[0])),
                                src[1] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - dst[1]) << 8) / src[1])),
                                src[2] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - dst[2]) << 8) / src[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case COLOR_DODGE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0] == 255 ? 255 :
                                    Math.min((dst[0] << 8) / (255 - src[0]), 255),
                                src[1] == 255 ? 255 :
                                    Math.min((dst[1] << 8) / (255 - src[1]), 255),
                                src[2] == 255 ? 255 :
                                    Math.min((dst[2] << 8) / (255 - src[2]), 255),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case DARKEN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.min(src[0], dst[0]),
                                Math.min(src[1], dst[1]),
                                Math.min(src[2], dst[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case DIFFERENCE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.abs(dst[0] - src[0]),
                                Math.abs(dst[1] - src[1]),
                                Math.abs(dst[2] - src[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case EXCLUSION:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] + src[0] - (dst[0] * src[0] >> 7),
                                dst[1] + src[1] - (dst[1] * src[1] >> 7),
                                dst[2] + src[2] - (dst[2] * src[2] >> 7),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case FREEZE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0] == 0 ? 0 : Math.max(0, 255 - (255 - dst[0]) * (255 - dst[0]) / src[0]),
                                src[1] == 0 ? 0 : Math.max(0, 255 - (255 - dst[1]) * (255 - dst[1]) / src[1]),
                                src[2] == 0 ? 0 : Math.max(0, 255 - (255 - dst[2]) * (255 - dst[2]) / src[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case GLOW:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] == 255 ? 255 : Math.min(255, src[0] * src[0] / (255 - dst[0])),
                                dst[1] == 255 ? 255 : Math.min(255, src[1] * src[1] / (255 - dst[1])),
                                dst[2] == 255 ? 255 : Math.min(255, src[2] * src[2] / (255 - dst[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case GREEN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0],
                                dst[1],
                                src[2],
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case HARD_LIGHT:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0] < 128 ? dst[0] * src[0] >> 7 :
                                    255 - ((255 - src[0]) * (255 - dst[0]) >> 7),
                                src[1] < 128 ? dst[1] * src[1] >> 7 :
                                    255 - ((255 - src[1]) * (255 - dst[1]) >> 7),
                                src[2] < 128 ? dst[2] * src[2] >> 7 :
                                    255 - ((255 - src[2]) * (255 - dst[2]) >> 7),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case HEAT:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] == 0 ? 0 : Math.max(0, 255 - (255 - src[0]) * (255 - src[0]) / dst[0]),
                                dst[1] == 0 ? 0 : Math.max(0, 255 - (255 - src[1]) * (255 - src[1]) / dst[1]),
                                dst[2] == 0 ? 0 : Math.max(0, 255 - (255 - src[2]) * (255 - src[2]) / dst[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case HUE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            float[] srcHSL = new float[3];
                            RGBtoHSL(src[0], src[1], src[2], srcHSL);
                            float[] dstHSL = new float[3];
                            RGBtoHSL(dst[0], dst[1], dst[2], dstHSL);

                            int[] result = new int[4];
                            HSLtoRGB(srcHSL[0], dstHSL[1], dstHSL[2], result);
                            result[3] = Math.min(255, src[3] + dst[3]);

                            return result;
                        }
                    };
                case INVERSE_COLOR_BURN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - src[0]) << 8) / dst[0])),
                                dst[1] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - src[1]) << 8) / dst[1])),
                                dst[2] == 0 ? 0 :
                                    Math.max(0, 255 - (((255 - src[2]) << 8) / dst[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case INVERSE_COLOR_DODGE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] == 255 ? 255 :
                                    Math.min((src[0] << 8) / (255 - dst[0]), 255),
                                dst[1] == 255 ? 255 :
                                    Math.min((src[1] << 8) / (255 - dst[1]), 255),
                                dst[2] == 255 ? 255 :
                                    Math.min((src[2] << 8) / (255 - dst[2]), 255),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case LIGHTEN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.max(src[0], dst[0]),
                                Math.max(src[1], dst[1]),
                                Math.max(src[2], dst[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case LUMINOSITY:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            float[] srcHSL = new float[3];
                            RGBtoHSL(src[0], src[1], src[2], srcHSL);
                            float[] dstHSL = new float[3];
                            RGBtoHSL(dst[0], dst[1], dst[2], dstHSL);

                            int[] result = new int[4];
                            HSLtoRGB(dstHSL[0], dstHSL[1], srcHSL[2], result);
                            result[3] = Math.min(255, src[3] + dst[3]);

                            return result;
                        }
                    };
                case MULTIPLY:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                (src[0] * dst[0]) >> 8,
                                (src[1] * dst[1]) >> 8,
                                (src[2] * dst[2]) >> 8,
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case NEGATION:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                255 - Math.abs(255 - dst[0] - src[0]),
                                255 - Math.abs(255 - dst[1] - src[1]),
                                255 - Math.abs(255 - dst[2] - src[2]),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case OVERLAY:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] < 128 ? dst[0] * src[0] >> 7 :
                                    255 - ((255 - dst[0]) * (255 - src[0]) >> 7),
                                dst[1] < 128 ? dst[1] * src[1] >> 7 :
                                    255 - ((255 - dst[1]) * (255 - src[1]) >> 7),
                                dst[2] < 128 ? dst[2] * src[2] >> 7 :
                                    255 - ((255 - dst[2]) * (255 - src[2]) >> 7),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case RED:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0],
                                dst[1],
                                dst[2],
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case REFLECT:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                src[0] == 255 ? 255 : Math.min(255, dst[0] * dst[0] / (255 - src[0])),
                                src[1] == 255 ? 255 : Math.min(255, dst[1] * dst[1] / (255 - src[1])),
                                src[2] == 255 ? 255 : Math.min(255, dst[2] * dst[2] / (255 - src[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case SATURATION:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            float[] srcHSL = new float[3];
                            RGBtoHSL(src[0], src[1], src[2], srcHSL);
                            float[] dstHSL = new float[3];
                            RGBtoHSL(dst[0], dst[1], dst[2], dstHSL);

                            int[] result = new int[4];
                            HSLtoRGB(dstHSL[0], srcHSL[1], dstHSL[2], result);
                            result[3] = Math.min(255, src[3] + dst[3]);

                            return result;
                        }
                    };
                case SCREEN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                255 - ((255 - src[0]) * (255 - dst[0]) >> 8),
                                255 - ((255 - src[1]) * (255 - dst[1]) >> 8),
                                255 - ((255 - src[2]) * (255 - dst[2]) >> 8),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case SOFT_BURN:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] + src[0] < 256 ?
                                        (dst[0] == 255 ? 255 :
                                        Math.min(255, (src[0] << 7) / (255 - dst[0]))) :
                                            Math.max(0, 255 - (((255 - dst[0]) << 7) / src[0])),
                                dst[1] + src[1] < 256 ?
                                        (dst[1] == 255 ? 255 :
                                        Math.min(255, (src[1] << 7) / (255 - dst[1]))) :
                                            Math.max(0, 255 - (((255 - dst[1]) << 7) / src[1])),
                                dst[2] + src[2] < 256 ?
                                        (dst[2] == 255 ? 255 :
                                        Math.min(255, (src[2] << 7) / (255 - dst[2]))) :
                                            Math.max(0, 255 - (((255 - dst[2]) << 7) / src[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case SOFT_DODGE:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                dst[0] + src[0] < 256 ?
                                    (src[0] == 255 ? 255 :
                                        Math.min(255, (dst[0] << 7) / (255 - src[0]))) :
                                            Math.max(0, 255 - (((255 - src[0]) << 7) / dst[0])),
                                dst[1] + src[1] < 256 ?
                                    (src[1] == 255 ? 255 :
                                        Math.min(255, (dst[1] << 7) / (255 - src[1]))) :
                                            Math.max(0, 255 - (((255 - src[1]) << 7) / dst[1])),
                                dst[2] + src[2] < 256 ?
                                    (src[2] == 255 ? 255 :
                                        Math.min(255, (dst[2] << 7) / (255 - src[2]))) :
                                            Math.max(0, 255 - (((255 - src[2]) << 7) / dst[2])),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case SOFT_LIGHT:
                    break;
                case STAMP:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.max(0, Math.min(255, dst[0] + 2 * src[0] - 256)),
                                Math.max(0, Math.min(255, dst[1] + 2 * src[1] - 256)),
                                Math.max(0, Math.min(255, dst[2] + 2 * src[2] - 256)),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
                case SUBTRACT:
                    return new Blender() {
                        @Override
                        public int[] blend(int[] src, int[] dst) {
                            return new int[] {
                                Math.max(0, src[0] + dst[0] - 256),
                                Math.max(0, src[1] + dst[1] - 256),
                                Math.max(0, src[2] + dst[2] - 256),
                                Math.min(255, src[3] + dst[3])
                            };
                        }
                    };
            }
            throw new IllegalArgumentException("Blender not implement for " +
                                               composite.getMode().name());
        }
    }
}

But before I go I have a couple of words of caution about this technique:

1)Intensity values are 'added' they are not averaged or interpolated so this is NOT an alternative to krigging or gausian processes

Only use this technique to generate heat maps of densities!!
e.g. population,clicks,customers,suicides
Do not use this technique to generate heatmaps of values.
e.g. house prices, travel times

try putting a low value next to a high value to see what I mean!

2)There are no error bars or probability estimates.
Furthermore the higher the values of a datapoint or (equivelently the more datapoints) at a specific coordinate the more the blob spreads. This may wrongly imply that nearby areas have a higher density than they do. e.g. in a map of population density a sky scraper in the middle of a park would not appear as a bright white but small dot surrounded by blue but rather the park would appear red.

try putting loads of points in the same place to see what I mean.