"Heat map" advice/opinion

sjr51 wrote: Hello all,

I have a set of 2D coords and I'd like to plot the density of these points as a "heat map" (colour-coded map of probability at every location within a boundary). In the past, I have done one of two things:
1. 2D histogram - sample number of points in a grid and use that for representation, but this can look blocky. I have around 200 points so reducing bins can look weird with smaller bin sizes.
2. ImageInterpolate (voronoi) - for each point, calculate the number of neighbouring points within some arbitrary distance and use these values as z to make a surface with

ImageInterpolate

with voronoi option.

I'm wondering if a better approach would be to make a 2D gauss at each location and sum these together (perhaps with smoothing) to make the representation.
I'm not sure how to go about doing this other approach. I've tried to investigate the new

StatsKDE

option in IP7 without much success - does it work on 2D waves?

If anyone has:
- an opinion of what is the best way to make a heat map in these circumstances
- advice on how to do the third option.
I'd love to hear it. Otherwise, feel free to tell me to take this to stackexchange instead!

//******************************************************************************************************
// Makes a 2D histogram with Gaussian smoothing
Function TwoDHist_Gaussian (datasetnamestr, OutputFolderPath, Xmin, Xmax, Ymin, Ymax, pixelsize, GausRadius, DisplayNow)
	String datasetnamestr, OutputFolderPath
	variable Xmin, Xmax, Ymin, Ymax,  pixelsize, GausRadius, DisplayNow
	
	// make duplicates of the X and Y waves and trim them to Xmin, Xmax, Ymin, Ymax, plus a border equivalent to 3X the GausRadius
	// because cells that far away could affect pixels in the desired image space. Aterwards, we will crop the image.
	wave Xwave = $datasetnamestr + "x"
	wave Ywave = $datasetnamestr + "y"
	if (!(datafolderexists ("root:packages:twoDHist")))
		if (!(DataFolderExists ("root:packages")))
			newDataFolder root:packages
		endif
		newdatafolder root:packages:twoDHist
	endif
	duplicate/o Xwave root:packages:twoDHist:XWave
	duplicate/o Ywave root:packages:twoDHist:YWave
	Wave Xwave = root:packages:twoDHist:XWave
	Wave Ywave = root:packages:twoDHist:YWave
	
	sort  Xwave Xwave, Ywave
	// have to get rid of NaNs, they will be at the end after a sort
	variable ii
	For (ii = numpnts (Xwave) -1; numtype (Xwave [ii]) > 0; ii -= 1)
		Deletepoints ii, 1, Xwave, Ywave
	EndFor
	
	variable bs = BinarySearch(Xwave, xmin - (GausRadius * pixelsize * 3 ))
	if (bs > -1)
		DeletePoints 0, bs, Xwave, Ywave
	endif
	bs = BinarySearch(Xwave, Xmax + (GausRadius  * pixelsize * 3 ))
	if (bs > -1)
		DeletePoints bs, numpnts (Xwave), Xwave, Ywave
	endif

	sort Ywave Ywave, Xwave
	bs= BinarySearch(Ywave, ymin - (GausRadius  * pixelsize * 3))
	if (bs > -1)
		DeletePoints 0, bs, Xwave, Ywave
	endif
	bs= BinarySearch(Ywave, Ymax + (GausRadius * pixelsize  * 3))
	if (bs > -1)
		DeletePoints bs, numpnts (Ywave), Xwave, Ywave
	endif
	
	// Make a Gausian Blob of the required size.  After 3X the radius, the Gaussian function has fallen to about 3%, so we'll cut it off there
	// and make the blob size 6 times the radius, plus 1 for the center pixel. We'll also blank out the corners of the square image.
	variable GausPixWid = GausRadius * 6 + 1
	make/o/n = ((GausPixWid), (GausPixWid)) root:packages:twoDHist:GausBlob
	WAVE GausBlob = root:packages:twoDHist:GausBlob
	GausBlob = ((p - 3 * GausRadius)^2 + (q - 3 * GausRadius) ^2) <  (3 * GausRadius)^2 ? e^-(((p - 3 * GausRadius)^2 + (q - 3 * GausRadius)^2)/(2 * GausRadius^2)): 0
	
	// Make an  image bigger than we need, because of the margin of interaction issue
	variable xpixels = ceil((Xmax-Xmin)/pixelSize) + GausPixWid -1
	variable ypixels = ceil((Ymax-Ymin)/pixelSize) + GausPixWid -1
	make/o/n=((xpixels), (ypixels)), $OutputFolderPath
	Wave celim = $OutputFolderPath
	celim = 0
	
	//Now fill the image with the data
	// first transform the X,Y coordinates from mm into pixel space
	Xwave -= (Xmin  - (GausRadius  * pixelsize * 3)+ pixelsize/2)
	Ywave -= (Ymin  - (GausRadius  * pixelsize * 3) + pixelsize/2)
	Xwave /=PixelSize
	Ywave/=PixelSize
	variable numdata = numpnts (Xwave)
	For (ii = 0; ii < numdata; ii += 1)
		celim [max (0, (Xwave [ii] - (3 * GausRadius))), min (Xpixels-1, (Xwave [ii] +(3 *  GausRadius)))] [max (0, (Ywave [ii] - (3 * GausRadius))), min (Ypixels-1, (Ywave [ii] + (3 *GausRadius)))] +=  GausBlob [p - Xwave [ii] + (3 *GausRadius)] [q -Ywave [ii] + (3 *GausRadius)]
	endfor
	killwaves/z Xwave, Ywave, GausBlob		//  kill the temporary waves
	
	//delete the extra margins of the image
	DeletePoints /M=0  (xpixels - (3 *GausRadius)), (3 *GausRadius), celim
	DeletePoints /M=0  0, (3 * GausRadius), celim
	DeletePoints /M=1  (ypixels - (3 * GausRadius)), (3 *GausRadius), celim
	DeletePoints /M=1  0, (3 *GausRadius), celim
	SetScale/P x Xmin + PixelSize/2, PixelSize,"mm", celim
	SetScale/P y Ymin + PixelSize/2, PixelSize,"mm", celim
	
	// Display the image if so requested
	if (DisplayNow)
		display;appendimage celim
		ModifyGraph width={Plan,1,bottom,left}
		ModifyImage $nameofwave (celim) ctab= {*,*,Rainbow,1}
	endif
end

•make/n = 250 testx, testy
•testx = enoise (5); testy = enoise (5) + 10
•TwoDHist_Gaussian ("test", "root:OutPut", -6, 6, 3, 16, .05, 9, 1)
•append testy vs testx
•ModifyGraph mode=3,rgb=(0,0,0)

GausBlob = ((p - 3 * GausRadius)^2 + (q - 3 * GausRadius) ^2) <  (3 * GausRadius)^2 ? e^-(((p - 3 * GausRadius)^2 + (q - 3 * GausRadius)^2)/(2 * GausRadius^2)): 0 // after this line
variable blobsum = sum(GausBlob)  //add this line
GausBlob /= blobsum  // and this line