Home › Code Snippets

2D (Joint) Histogram

Posted February 18th, 2010 by RGerkin
in
Average rating
(2 votes)

This function creates a joint histogram using two input waves, which must be of equal size. The result will be placed in the wave 'hist', which must be made beforehand and scaled to determine the histogram binning and boundaries. For example, 

make /o/n=1000 data1=gnoise(1), data2=gnoise(1)
make /o/n=(25,25) myHist
setscale x,-3,3,myHist
setscale y,-3,3,myHist
JointHistogram(data1,data2,myHist)
NewImage myHist

This is the function, which runs slightly differently depending on how many cores your processor has. In my experience, the MatrixOp version is faster unless you have 4 or more cores, in which case the Multithreaded Make version is faster. You don't have to worry about this, as it is taken care of automatically.

Function JointHistogram(w0,w1,hist)
	wave w0,w1,hist
 
	variable bins0=dimsize(hist,0)
	variable bins1=dimsize(hist,1)
	variable n=numpnts(w0)
	variable left0=dimoffset(hist,0)
	variable left1=dimoffset(hist,1)
	variable right0=left0+bins0*dimdelta(hist,0)
	variable right1=left1+bins1*dimdelta(hist,1)
 
	// Scale between 0 and the number of bins to create an index wave.  
	if(ThreadProcessorCount<4) // For older machines, matrixop is faster.  
		matrixop /free idx=round(bins0*(w0-left0)/(right0-left0))+bins0*round(bins1*(w1-left1)/(right1-left1))
	else // For newer machines with many cores, multithreading with make is faster.  
		make /free/n=(n) idx
		multithread idx=round(bins0*(w0-left0)/(right0-left0))+bins0*round(bins1*(w1-left1)/(right1-left1))
	endif
 
	// Compute the histogram and redimension it.  
	histogram /b={0,1,bins0*bins1} idx,hist
	redimension /n=(bins0,bins1) hist // Redimension to 2D.  
	setscale x,left0,right0,hist // Fix the histogram scaling in the x-dimension.  
	setscale y,left1,right1,hist // Fix the histogram scaling in the y-dimension.  
End

Back to top