
% Adam Wang, 17 Mar 2008
% Performs localized NPS analysis based on wavelet decomposition

% Inputs: im - image, wc - (x,y) coord of center of ROI window, ws - window
% size (should be power of 2).

% Outputs: NPS - power spectrum at freqs in vector rho, rho - frequencies

function [NPS rho] = NPS_ROI(im, wc, ws)

levels = log(ws)/log(2);

imSize = size(im,1);

W = Wavelet; % need @Wavelet folder!
imW = W*im;

imW(1:imSize/ws, 1:imSize/ws) = W'*imW(1:imSize/ws, 1:imSize/ws);

rho = zeros(1,1+2*levels);
NPS = zeros(1,1+2*levels);

for i = 1:levels

    rows = round((wc(1)-1)/2^i+1 - ws/2^(i+1) + 1/2):round((wc(1)-1)/2^i+1 + ws/2^(i+1) - 1/2); % finds center of window within subblock
    cols = round((wc(2)-1)/2^i+1 - ws/2^(i+1) + 1/2):round((wc(2)-1)/2^i+1 + ws/2^(i+1) - 1/2);

    NPS(i*2-1) = mean(mean(imW(imSize/2^i+rows, imSize/2^i+cols).^2)); % D coeffs
    rho(i*2-1) = (imSize/2^i+1 + imSize/2^(i-1))/2 * sqrt(2);
    
    NPS(i*2) = mean(mean([imW(rows, imSize/2^i+cols) imW(imSize/2^i+rows, cols)].^2)); % H and V coeffs
    rho(i*2) = norm([(imSize/2^i+1 + imSize/2^(i-1))/2 (1 + imSize/2^i)/2]); 
end

rho(end) = (1 + imSize/2^i)/2 * sqrt(2); % highest-level A coeffs
NPS(end) = imW(rows,cols)^2;

rho = rho/2;

rho = fliplr(rho);
NPS = fliplr(NPS);

    