;compute w(theta) from an input image and random catalog
function get_wtheta_img, img, random, non_masked_img, rmin, rmax, nbin, spix

nxpix = n_elements(img[*,0])
nypix = n_elements(img[0,*])


;Get mean
mean_img = mean(img[non_masked_img])
;remove mean and normalize by mean (= compute the overdensity)
img[non_masked_img] -= mean_img
img[non_masked_img] /= mean_img

;Add zero values around
img_zero    = zero_around(img)
random_zero = zero_around(random)


ft_img    = fft(img_zero, /double)
ft_random = fft(random_zero, /double)
img_zero    = 0
random_zero = 0

square_mod_img    = (abs(ft_img))^2
ft_img = 0
square_mod_random = (abs(ft_random))^2
ft_random = 0


xi_img    = real_part(fft(square_mod_img, /inverse, /double))
xi_random = real_part(fft(square_mod_random, /inverse, /double))

;get xi(r)
nxpix = n_elements(xi_img[*,0])
nypix = n_elements(xi_img[0,*])
xi_r  = dblarr(nxpix, nypix)

use       = where(xi_random ne 0, count)
xi_r[use] = xi_img[use]/xi_random[use]

;shift subscripts
xi_r      = shift(xi_r,nxpix/2 -1L, nypix/2 -1L)
xi_img    = shift(xi_img,nxpix/2 -1L, nypix/2 -1L)
xi_random = shift(xi_random,nxpix/2 -1L, nypix/2 -1L)

;get position array
mx  = indgen(nxpix) - (nxpix/2 -1)
rx = mx*spix

my  = indgen(nypix) - (nypix/2 -1)
ry  = my*spix


;Compute distances (in order to integrate xi(r) over r)
dist = dblarr(nxpix, nypix)
for iy = 0L, nypix -1L do begin
    dist[*, iy] = ((rx)^2 + (ry[iy])^2)
endfor
dist = sqrt(dist)

;Set up log binning
rmin_log  = alog10(rmin)
rmax_log  = alog10(rmax)
bin_log   = (rmax_log-rmin_log)/nbin
r_arr     = 10^(rmin_log +bin_log*dindgen(nbin +1L))


num        = dblarr(nbin)
den        = dblarr(nbin)
;Loop over bins in r and integrate 'xi(r)'
if (keyword_set(verbose)) then print, 'Integrate to get w(theta)'
for ibin = 0, nbin-1 do begin
    search = where(dist GT r_arr[ibin] and dist LE r_arr[ibin +1])
    dr               = r_arr[ibin +1] - r_arr[ibin]
    dr2              = dist[search]*dr*2.*!dpi
    num[ibin]        = total(xi_img[search]*dr2)
    den[ibin]        = total(xi_random[search]*dr2)
endfor

wtheta    = num/den

;stop
return, wtheta
end