
%Adam Wang, 17 Mar 2008.
% Create 200 realizations of CT noise, low band, mid band, and high band
% noise.
% Find average wavelet power of each.
% Analyze local 

randn('state',2);

ph = phantom(256);
sino = radon(ph, 0:.5:359.5)';
recon = iradon(sino', 0:.5:359.5);

LPF = [zeros(241,1); ones(239,1); zeros(240,1)]; % low band (freq space)
MPF = [zeros(121,1); ones(120,1); zeros(239,1); ones(120,1); zeros(120,1)]; % mid band (freq space)
HPF = [ones(121,1); zeros(479,1); ones(120,1)]; % high band (freq space)

lpf = ifftc(LPF);
mpf = ifftc(MPF);
hpf = ifftc(HPF);

wc = [128.5 128.5;
      70.5  128.5;
      20.5  128.5;
      54.5  202.5;];
ws = 32;
num_realizations = 100;

num_ROI = size(wc,1); len = log(ws)/log(2)*2+1;
NPS_W = zeros(num_ROI,len); NPS_L = zeros(num_ROI,len);
NPS_M = zeros(num_ROI,len); NPS_H = zeros(num_ROI,len);

imageWhiteAll = zeros(256,256,200); imageLowAll = zeros(256,256,200);
imageHiAll = zeros(256,256,200); imageMidAll = zeros(256,256,200);
for i = 1:200
    tic
    whiteSino = randn(size(sino));
    hybrid = fftshift(fft(fftshift(whiteSino,1),[],1),1);

    filtSinoH = fftshift(ifft(fftshift(hybrid.*(HPF*ones(1,size(whiteSino,2))),1),[],1),1); % high temporal freq sinogram
    filtSinoL = fftshift(ifft(fftshift(hybrid.*(LPF*ones(1,size(whiteSino,2))),1),[],1),1); % low temporal freq sinogram
    filtSinoM = fftshift(ifft(fftshift(hybrid.*(MPF*ones(1,size(whiteSino,2))),1),[],1),1); % mid temporal freq sinogram

    imageWhiteAll(:,:,i) = iradon2(whiteSino', 0:.5:359.5, 'sinc', 256); % recon and store CT noise
    imageLowAll(:,:,i) = iradon2(filtSinoL', 0:.5:359.5, 'sinc', 256);  % recon and store low band noise
    imageHiAll(:,:,i) = iradon2(filtSinoH', 0:.5:359.5, 'sinc', 256);   % recon and store high band noise
    imageMidAll(:,:,i) = iradon2(filtSinoM', 0:.5:359.5, 'sinc', 256);  % recon and store mid band noise

    if(mod(i,20)==1) disp([i toc]); end
end
disp([i toc]);

folder = 'C:\Documents and Settings\Adam Wang\My Documents\_Second Year\PSYCH221\Project\data\';
imwrite(uint8((whiteSino+4)/8*256), [folder 'whiteSino.png'], 'png');
imwrite(uint8((filtSinoH+4)/8*256), [folder 'filtSinoH.png'], 'png');
imwrite(uint8((filtSinoM+4)/8*256), [folder 'filtSinoM.png'], 'png');
imwrite(uint8((filtSinoL+4)/8*256), [folder 'filtSinoL.png'], 'png');

imwrite(uint8((squeeze(imageWhiteAll(:,:,end))+.1)/.2*256), [folder 'imageWhite.png'], 'png');
imwrite(uint8((squeeze(imageHiAll(:,:,end))+.1)/.2*256), [folder 'imageHi.png'], 'png');
imwrite(uint8((squeeze(imageMidAll(:,:,end))+.1)/.2*256), [folder 'imageMid.png'], 'png');
imwrite(uint8((squeeze(imageLowAll(:,:,end))+.1)/.2*256), [folder 'imageLow.png'], 'png');


% whiteF1 = zeros(16,16); hiF1 = zeros(16,16); lowF1 = zeros(16,16); midF1 = zeros(16,16);
% for i = 1:200
%     whiteF1 = whiteF1 + abs(fft2c(squeeze(imageWhiteAll(121:136,121:136,i)))).^2;
%     hiF1 = hiF1 + abs(fft2c(squeeze(imageHiAll(121:136,121:136,i)))).^2;
%     midF1 = midF1 + abs(fft2c(squeeze(imageMidAll(121:136,121:136,i)))).^2;
%     lowF1 = lowF1 + abs(fft2c(squeeze(imageLowAll(121:136,121:136,i)))).^2;
% end
% figure, subplot(221); imshow(whiteF1,[]); subplot(222); imshow(hiF1,[]); subplot(223); imshow(midF1,[]); subplot(224); imshow(lowF1,[]);
% whiteF2 = zeros(16,16); hiF2 = zeros(16,16); lowF2 = zeros(16,16); midF2 = zeros(16,16);
% for i = 1:200
%     whiteF2 = whiteF2 + abs(fft2c(squeeze(imageWhiteAll(20:35,121:136,i)))).^2;
%     hiF2 = hiF2 + abs(fft2c(squeeze(imageHiAll(20:35,121:136,i)))).^2;
%     midF2 = midF2 + abs(fft2c(squeeze(imageMidAll(20:35,121:136,i)))).^2;
%     lowF2 = lowF2 + abs(fft2c(squeeze(imageLowAll(20:35,121:136,i)))).^2;
% end
% figure, subplot(221); imshow(whiteF2,[]); subplot(222); imshow(hiF2,[]); subplot(223); imshow(midF2,[]); subplot(224); imshow(lowF2,[]);


% Compute average wavelet power of each.
whiteW = zeros(256,256);
for i = 1:200
    whiteW = (W*squeeze(imageWhiteAll(:,:,i))).^2 + whiteW;
end
figure, imshow(whiteW,[])
lowW = zeros(256,256);
for i = 1:200
    lowW = (W*squeeze(imageLowAll(:,:,i))).^2 + lowW;
end
hiW = zeros(256,256);
for i = 1:200
    hiW = (W*squeeze(imageHiAll(:,:,i))).^2 + hiW;
end
midW = zeros(256,256);
for i = 1:200
    midW = (W*squeeze(imageMidAll(:,:,i))).^2 + midW;
end


% Compute local NPS for center and edge ROIs from average wavelet power
[NPS rho] = NPS_ROI_(midW,[128.5 128.5],32)
[NPS rho] = NPS_ROI_(highW,[128.5 128.5],32)
[NPS rho] = NPS_ROI_(hiW,[128.5 128.5],32)
[NPS rho] = NPS_ROI_(lowW,[128.5 128.5],32)
[NPS rho] = NPS_ROI_(whiteW,[220 200],32); NPS*1000
[NPS rho] = NPS_ROI_(lowW,[200 200],32); NPS*1000
[NPS rho] = NPS_ROI_(midW,[200 200],32); NPS*1e6
[NPS rho] = NPS_ROI_(hiW,[200 200],32); NPS*1e6

