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