clear
% Number of Blocks in Grid
N=5;
M=5;
%Number of Levels - Slows after 3!
itermax=3;
% Initializations
dc=1/N;
Data1=zeros(N*M,6);
dx=1/N;
dy=1/M;
% Initial Blocks
figure(1)
hold on
cnt=0;
for i=1:N
    for j=1:M
        cnt=cnt+1;
        Data1(cnt,:)=[dx 0 0 dy (i-1)*dx (j-1)*dy];
        v=[(i-1)*dx (j-1)*dy; i*dx (j-1)*dy; i*dx j*dy; (i-1)*dx j*dy;];
        patch(v(:,1), v(:,2),[dc*i dc*j 0]);
        flg(i,j)=1;
    end
end
axis([0 1 0 1])
axis square
hold off
% Edit Blocks
[x0 y0 bb]=ginput(1); 
hold on
while bb<2,
    if bb==1
        for i=1:N
            if (x0>(i-1)*dx)&(x0<i*dx) 
                i0=i;
            end      
        end
        for j=1:M
            if (y0>(j-1)*dy)&(y0<j*dy) 
                j0=j;
            end      
        end   
        flg(i0,j0)=1-flg(i0,j0);
        fl=flg(i0,j0);
        cnt=j0+(i0-1)*M;
        v=[(i0-1)*dx (j0-1)*dy; i0*dx (j0-1)*dy; i0*dx j0*dy; (i0-1)*dx j0*dy;];
        patch(v(:,1), v(:,2),[(1-fl)+fl*dc*i0 (1-fl)+fl*dc*j0 1-fl]);
        Data1(cnt,:)=[fl*dx 0 0 fl*dy fl*(i0-1)*dx fl*(j0-1)*dy];
    end 
 [x0 y0,bb]=ginput(1);
end

tmp=[];
tmp2=[];
for i=1:N*M
    if Data1(i,:)==[0 0 0 0 0 0]
           tmp=[tmp i]; 
    else
        tmp2=[tmp2 i];
    end
end

% New Transformations 
Data=Data1(tmp2,:);

figure(2)
% Initial Rectangle
%v=[0 0; dx 0; dx dy; 0 dy];
v=[0 0; 1 0; 1 1 ;0 1];
Npts=size(v,1);
z=v;
patch(v(:,1), v(:,2),[1 0 0]);
axis([0 1 0 1])
hold

% Linear Transformations
A=[dx 0; 0 dy];

% Translations
B=Data(:,5:6);

% Cycle Through Levels
disp('Computing ...')
h = waitbar(0,'Please wait... ','Name','Fractal Blocks'); 
t0=cputime;

NTrans=length(Data);

for iter=1:itermax
    Mtri=size(z,1)/Npts;   % Number of blocks
    for j=1:Mtri
        waitbar((iter-1+j/Mtri)/itermax ,h ,['Iteration = ' num2str(iter) ' Rectangle = ' num2str(j)])
        v=z(4*(j-1)+1:4*j,:);
        for k=1:NTrans
            w=A*v';
            w=w'+repmat(B(k,:),4,1);
            if (j==1)&(k==1)
                z2=w;
            else
                z2=[z2;w];
            end
            patch(w(:,1), w(:,2),[1-0.2*iter 0.2*iter 0.2*iter]);
        end
    end
    z=z2;
end
disp('CPU Time =')
cputime-t0
close(h)
disp('Finishing up ...')
hold