function A=fasthaar2(show,A);
%
% first attempt of a FAST 2dim Haar wavelet transform
% ovrwriting the original matrix
% working from the finest level up, replace each
% 2x2 block by average in top left, and differences
% along horizontal, vertical, diagonal in other 3
% this is implemented by reshape and matrixmult
% then recursively consider only top left entries
% of each block, and do THE SAME operation to each
%
% eventually want to code the inverse similarily
% and to display basis elements and lin comb's
% as rastered images
% quality measure is energy in each frequency
% and in vertical / horizontal variation and
% diagonal dominance
%
% first attempt August 2008
n=5;
N=2^n;
if nargin < 1; show=1; else show=0; end
if nargin < 2
   q=double(imread('einstein.gif'));
   q=64-63*q/max(max(q));
   colormap(gray);
   image(q);
   A=q(46-N/2:45+N/2,36-N/2:35+N/2)/32-1;
end
image(32*(A+1))
pause
T=[1 1 1 1;1 1 -1 -1;1 -1 1 -1;1 -1 -1 1]/2;
for k=0:n-1
    for i=0:2^(k+1):N-2
    for j=0:2^(k+1):N-2
        M=[i+1:i+2,j+1:j+2];
        A(i+1:2^k:i+2^k+1,j+1:2^k:j+2^k+1)=...
            reshape(T*reshape(A(i+1:2^k:i+2^k+1,j+1:2^k:j+2^k+1),4,1),2,2);
        if show, image(32*(A+1)); pause; end;
    end
    end
end
end