"Heat map" advice/opinion
sjr51
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!
You might also want to check out the ImageFromXYZ operation.
FWIW: I think there is value in displaying the raw joint histogram results. The interpolation could lead you to misinterpret the results.
October 28, 2016 at 10:14 am - Permalink
October 28, 2016 at 10:33 pm - Permalink
This is my code from years ago to do the third option you suggest. See here, for published example involving 2D maps of labeled neurons in brain sections
It assumes that x and y wave names are the same, except for ending with "x" and "y" respectively, and are generated from dataSetNameStr
Xmin, Xmax, Ymin, Ymax are boundaries of the output image
Pixelsize for output image is in whatever units your data are in, and Gaussian radius is the radius of the Gaussian kernel in pixels
// 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
I have attached a picture of some random data ovulated on its 2D gaussian blobbed equivalent
•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)
Jamie Boyd, Ph.D.
October 31, 2016 at 05:28 pm - Permalink
November 1, 2016 at 01:04 am - Permalink
variable blobsum = sum(GausBlob) //add this line
GausBlob /= blobsum // and this line
so that the sum of the image will be equal to the number of points contributing to the image, for whatever size of kernel you use.
Jamie Boyd, Ph.D.
November 1, 2016 at 10:03 am - Permalink