Centre of mass of an N-D image (1 < N < 3)

Function test()
	make/n=(100,100, 100)/free w
	w = gauss(x, 10, 3, y, 20, 3, z, 30, 3)
	wave temp = com3D(w)
	print temp
End
 
Function/wave com3D(w)
	Wave w
	variable ii, t1, t2, t3
	make/n=(wavedims(w))/free/d COM
	switch(wavedims(w))
		case 1:
			COM[0] = Pla_peakcentroid(w)
			break
		case 2:
			imagetransform sumallrows w
			COM[0] = Pla_peakcentroid(W_sumrows)
			imagetransform sumallcols w
			COM[1] = Pla_peakcentroid(W_sumcols)
			break
		case 3:
			imagetransform sumplanes w
			Wave temp = com3D(M_SumPlanes)
			COM[0] = temp[0]
			COM[1] = temp[1]
			imagetransform/g = 2 transposeVol w
			imagetransform sumplanes M_volumetranspose
                        imagetransform sumallrows M_sumplanes
			COM[2] = Pla_peakcentroid(W_sumrows)
			break
		default:
			break
	endswitch
	return COM
End
 
Threadsafe Function Pla_peakCentroid(ywave)
	Wave ywave
	variable retval
 
	duplicate/free ywave, W_integrate
	W_integrate = 0
 
	multithread W_integrate[] = Area(tempy, 0, p)
 
	W_integrate /= W_integrate[numpnts(W_tempy)-1]
	variable centroid = binarysearchinterp(W_integrate,0.5)
	if(numtype(centroid))
		retval = NaN
	else
		retval = centroid
	endif
 
	return retval
End