function [power_mean,power_var]=imagett_power(Image,nbx,nby,nx,ny)
%
%   [power_mean,power_var]=imagett_power(Image,nbx,nby,nx,ny)
%
% functions used to compute the power spectrum density of image
% input:
%       Image: the input image
%       nbx: subimage number in range /column
%       nby: subimage number in azimuth /row
%       nx: boxsize in range
%       ny: boxsize in azimuth
% output:
%       power_mean: the mean of PSD
%       power_var: the variance of PSD
%
power_mean=zeros(nx,ny); 
power_mean_s=zeros(nx,ny);
for i=1:nbx  
    for j=1:nby
        % seperate the image to several subimages
        subimage=Image((i-1)*nx+1:i*nx,(j-1)*ny+1:j*ny);
        subimage=reshape(subimage,nx,ny); 
        mean1=nanmean(nanmean(subimage(subimage(:)~=0))); 
        if isnan(mean1)
            continue
        end
        % compute the fft of subimage
        subimage(isnan(subimage(:)))=0;
        sub_fft=ifft2(subimage-mean1);
        sub_fft_abs2=abs(sub_fft).*abs(sub_fft);
        [endx,endy]=size(sub_fft_abs2);
        sub_fft_abs2_shift=sub_fft_abs2([endx-nx/2+1:endx,1:endx-nx/2],[endy-ny/2+1:endy,1:endy-ny/2]); 
        
        % compute the power spectrum density
        power=sub_fft_abs2_shift;
        power_mean=power_mean+power;
        power_mean_s=power_mean_s+power.*power;
        clear power sub_fft sub_fft_abs2 sub_fft_abs2_shift
    end
end
power_var=(power_mean_s-(power_mean.*power_mean)/(nbx*nby))/(nbx*nby-1);
power_mean=power_mean/(nbx*nby);
end





