clear
tic
Image1=imread('lena.bmp');

[imagem imagen]=size(Image1);
Sr=4;Sd=8;
Rnum=(imagem/Sr)*(imagen/Sr);
Dnum=(imagem/Sd)*(imagen/Sd);
Image2=zeros(Dnum,Sr,Sr);
Image2=blkproc(Image1,[Sd/Sr,Sd/Sr],'mean(mean(x))');
% there are no  eight tranformation for simpleness
RBlocks=zeros(Rnum,Sr,Sr);
DBlocks=zeros(Dnum,Sd,Sd);
DBlocksReduce=zeros(Dnum*8,Sr,Sr);
for i=1:imagem/Sr
    for j=1:imagen/Sr
        k=(i-1)*imagen/Sr+j;
        RBlocks(k,:,:)=Image1((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr);
       
    end
end
for i=1:imagem/Sd
    for j=1:imagen/Sd
        k=(i-1)*imagen/Sd+j;
        m=Sr;n=Sr;
        DBlocksReduce(k,:,:)=Image2((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr);
        DBlocksReduce(k+Dnum,:,:)=DBlocksReduce(k,m:-1:1,:);              % 行上下翻轉(zhuǎn)  y軸方向旋轉(zhuǎn)180度
        DBlocksReduce(k+2*Dnum,:,:)=DBlocksReduce(k,:,n:-1:1);              % 列左右翻轉(zhuǎn)  x軸方向旋轉(zhuǎn)180度
        DBlocksReduce(k+3*Dnum,:,:)=DBlocksReduce(k,m:-1:1,n:-1:1);         % 先行翻，再列翻  旋轉(zhuǎn)180度
        DBlocksReduce(k+4*Dnum,:,:)=reshape(DBlocksReduce(k,:,:),Sr,Sr)';                       % 關(guān)于y=-x對(duì)稱(chēng)  
        A=reshape( DBlocksReduce(k+3*Dnum,:,:),Sr,Sr)';
        DBlocksReduce(k+5*Dnum,:,:)=A(:,n:-1:1);     % 關(guān)于y=x對(duì)稱(chēng)
        DBlocksReduce(k+6*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),90);           % 逆時(shí)針旋轉(zhuǎn)90度   
        DBlocksReduce(k+7*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),270);          % 逆時(shí)針旋轉(zhuǎn)270度 
        
        
        DBlocks(k,:,:)=Image1((i-1)*Sd+1:i*Sd,(j-1)*Sd+1:j*Sd);
    end
end
RandDbest=zeros(Rnum,1)+256^3;
RandDbesta=zeros(Rnum,1);
RandDbestb=zeros(Rnum,1);
RandDbestj=zeros(Rnum,1);
for i=1:Rnum
    x=reshape(RBlocks(i,:,:),Sr*Sr,1);
    meanx=mean(x);
    for j=1:Dnum*8
        y=reshape(DBlocksReduce(j,:,:),Sr*Sr,1);
        meany=mean(y);
        a=(x-meanx)'*(y-meany)/((y-meany)'*(y-meany));
        b=(meanx-a*meany);
        c=(x-a*y-b)'*(x-a*y-b);
        if (RandDbest(i)>c)&(abs(a)<1)
            RandDbest(i)=c;
            RandDbesta(i)=a;
            RandDbestb(i)=b;
            RandDbestj(i)=j;
        end
    end
end
%iteration limit
toc
tic
m=8;
e=mean(mean(Image1));
Image3=e*ones(imagem,imagen);
for L=1:m
    Image4=blkproc(Image3,[Sd/Sr,Sd/Sr],'mean(mean(x))');
% 
%     for i=1:imagem/Sd
%         for j=1:imagen/Sd
%             k=(i-1)*imagen/Sd+j;
%             DBlocksReduce(k,:,:)=Image4((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr);
%             m=Sr;n=Sr;
%             DBlocksReduce(k,:,:)=Image4((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr);
%             DBlocksReduce(k+Dnum,:,:)=DBlocksReduce(k,m:-1:1,:);              % 行上下翻轉(zhuǎn)  y軸方向旋轉(zhuǎn)180度
%             DBlocksReduce(k+2*Dnum,:,:)=DBlocksReduce(k,:,n:-1:1);              % 列左右翻轉(zhuǎn)  x軸方向旋轉(zhuǎn)180度
%             DBlocksReduce(k+3*Dnum,:,:)=DBlocksReduce(k,m:-1:1,n:-1:1);         % 先行翻，再列翻  旋轉(zhuǎn)180度
%             DBlocksReduce(k+4*Dnum,:,:)=reshape(DBlocksReduce(k,:,:),Sr,Sr)';                       % 關(guān)于y=-x對(duì)稱(chēng)  
%             A=reshape( DBlocksReduce(k+3*Dnum,:,:),Sr,Sr)';
%             DBlocksReduce(k+5*Dnum,:,:)=A(:,n:-1:1);     % 關(guān)于y=x對(duì)稱(chēng)
%             DBlocksReduce(k+6*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),90);           % 逆時(shí)針旋轉(zhuǎn)90度   
%             DBlocksReduce(k+7*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),270);          % 逆時(shí)針旋轉(zhuǎn)270度 
%         end
%     end
    for i=1:imagem/Sr
        for j=1:imagen/Sr
            m=Sr;n=Sr;
            k=(i-1)*imagen/Sr+j;
            l=RandDbestj(k);
            k1=mod(l-1,Dnum)+1;
            l1=(l-k1)/Dnum+1;
%             if (k1==0)&(l1>0)
%                 k1=Dnum;
%                 l1=l1-1;
%             end
            %
%             k1=k1+1;
            %
            j1=mod(k1-1,imagen/Sd)+1;
            i1=(k1-j1)/(imagen/Sd)+1;
%             if (j1==0)&(i1>1)
%                 j1=imagen/Sd;
%                 i1=i1-1;
%             end
            DBlocksReduce(k1,:,:)=Image4((i1-1)*Sr+1:i1*Sr,(j1-1)*Sr+1:j1*Sr);
            switch l1-1
                case 0
                    DBlocksReduce(l,:,:)=Image4((i1-1)*Sr+1:i1*Sr,(j1-1)*Sr+1:j1*Sr);
                case 1
                    DBlocksReduce(l,:,:)=DBlocksReduce(k1,m:-1:1,:);  
                case 2
                    DBlocksReduce(l,:,:)=DBlocksReduce(k1,:,n:-1:1);  
                case 3
                    DBlocksReduce(l,:,:)=DBlocksReduce(k1,m:-1:1,n:-1:1);  
                case 4
                    DBlocksReduce(l,:,:)=reshape(DBlocksReduce(k1,:,:),Sr,Sr)'; 
                case 5
                    DBlocksReduce(k1+3*Dnum,:,:)=DBlocksReduce(k1,m:-1:1,n:-1:1); 
                    A=reshape( DBlocksReduce(k1+3*Dnum,:,:),Sr,Sr)';
                    DBlocksReduce(l,:,:)=A(:,n:-1:1); 
                case 6
                    DBlocksReduce(l,:,:)=imrotate(reshape(DBlocksReduce(k1,:,:),Sr,Sr),90); 
                case 7
                     DBlocksReduce(l,:,:)=imrotate(reshape(DBlocksReduce(k1,:,:),Sr,Sr),270); 
            end
            RBlocks(k,:,:)=RandDbesta(k)*DBlocksReduce(l,:,:)+RandDbestb(k);
            Image3((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr)=reshape(RBlocks(k,:,:),Sr,Sr);
        end
    end
    error=double(Image1)-Image3;
     Ps1(L)=20*log10(255/(sqrt(mean(mean(error.^2)))));
end
imshow(uint8(Image3))

toc
save('sa.mat')