% this function gives the colored noise variances at the different scales in
% the form of a image matrix. 
%
% differs from colored_variance function by using cdwt2_b
%
% inputs: 
%
% runs:                    number of runs used to determine variance
% N:                       size of image (assumed square)
% L:                       levels of decomposition
% h:                       error diffusion filter
% K:                       Floyd linear model gain
%
% usage [varmat1,varmat2] = colored_variances_b(sigman,runs,N,L,h,K)
%

function [varmat1,varmat2] = colored_variances_b(runs,N,L,h,K)

% create the  noise variance image. difft variance at different scales.
varmat1 = 0;
varmat2 = 0;

% For Floyd

noise_filter = -h;
xsize = size(noise_filter,2);
ysize = size(noise_filter,1);
noise_filter(1,ceil(xsize/2)) = 1;
noise_filter = noise_filter/K;

for iter = 1:runs
  junk = conv2(noise_filter, randn(N+1,N+1));
  junk = junk(2:2+N-1,4:4+N-1);
  %junk = real(ifft2(fft2(randn(N,N)) .* fft2(noise_filter)));
  [convmat1, convmat2] = cdwt2_b(junk, L);

  test = convmat1;
  for levs = 1:L
    start = N*(2^(-levs))+1;
    finish = N*(2^(-levs+1));
    
    % HH mat
    tempmat = (abs(test(start:finish,start:finish))).^2;
    varmat(start:finish,start:finish) = mean(tempmat(:));
    
    % HL mat
    tempmat = (abs(test(1:start-1,start:finish))).^2;
    varmat(1:start-1,start:finish) =  mean(tempmat(:));
    
    % LH mat
    tempmat = (abs(test(start:finish,1:start-1))).^2;
    varmat(start:finish,1:start-1) =  mean(tempmat(:));
  end
  tempmat = (abs(test(1:start-1,1:start-1))).^2;  
  varmat(1:start-1,1:start-1) = mean(tempmat(:));;
  varmat1 = varmat+varmat1 ;
  
  varmat = 0;
  test = convmat2;
  for levs = 1:L
    start = N*(2^(-levs))+1;
    finish = N*(2^(-levs+1));
    
    % HH mat
    tempmat = (abs(test(start:finish,start:finish))).^2;
    varmat(start:finish,start:finish) = mean(tempmat(:));
    
    % HL mat
    tempmat = (abs(test(1:start-1,start:finish))).^2;
    varmat(1:start-1,start:finish) =  mean(tempmat(:));
    
    % LH mat
    tempmat = (abs(test(start:finish,1:start-1))).^2;
    varmat(start:finish,1:start-1) =  mean(tempmat(:));
  end
  tempmat = (abs(test(1:start-1,1:start-1))).^2;  
  varmat(1:start-1,1:start-1) = mean(tempmat(:));;
  varmat2 = varmat+varmat2;
end

varmat1 = varmat1/iter;
varmat2 = varmat2/iter;


