分形压缩编码和解码程序

clear;
outfp = fopen（‘自己分形图像压缩.txt‘‘w‘）;


[DATAmap]=imread（‘E:\下载程序\2010年9月1日\分形图像压缩。\lena.bmp‘‘bmp‘）;
DATA=double（DATA）;
size1=size（DATA）;
tu=DATA（1:size1（1）1:size1（2）1）;

xx=size1（2）;
yy=size1（1）;
nrx=xx/8;
nry=yy/8;
ndx=xx/16;
ndy=yy/16;

DD=zeros（ndyndx864）;
cund1=zeros（1616）;
cund2=zeros（88）;
cund3=zeros（164）;
cund33=zeros（164）;
cund4=zeros（88）;

RRR=zeros（nrynrx6）;
% soD（y）D（x）U
cunr=zeros（164）;

for i=1:ndy
    for j=1:ndx
%         for k=1:16
            cund1=tu（1+16*（i-1）:16+16*（i-1）1+16*（j-1）:16+16*（j-1））;
            for l=1:8
                for m=1:8
                    cund2（lm）=（cund1（1+2*（l-1）1+2*（m-1））+cund1（2+2*（l-1）2+2*（m-1））+cund1（2+2*（l-1）1+2*（m-1））+cund1（1+2*（l-1）2+2*（m-1）））/4;
                end;
            end;
%         end;

        DD（ij11:64）=reshape（cund2[164]）;
        cund4=fliplr（cund2）;
        DD（ij21:64）=reshape（cund4[164]）;
        cund4=flipud（cund2）;
        DD（ij31:64）=reshape（cund4[164]）;
        cund4=flipud（fliplr（cund2））;
        DD（ij41:64）=reshape（cund4[164]）;
        cund4=rot90（flipud（cund2））;
        DD（ij51:64）=reshape（cund4[164]）;
        cund4=rot90（cund2）;
        DD（ij61:64）=reshape（cund4[164]）;
        cund4=rot90（rot90（rot90（cund2）））;
        DD（ij71:64）=reshape（cund4[164]）;
        cund4=cund2‘;
        DD（ij81:64）=reshape（cund4[164]）;
%         DD（ij11:81:8）=cund2;
%         DD（ij21:81:8）=fliplr（cund2）;
%         DD（ij31:81:8）=flipud（cund2）;
%         DD（ij41:81:8）=flipud（fliplr（cund2））;
%         DD（ij51:81:8）=Rot90（flipud（cund2））;
%         DD（ij61:81:8）=Rot90（cund2）;
%         DD（ij71:81:8）=Rot90（Rot90（Rot90（cund2）））;
%         DD（ij81:81:8）=cund2‘;
        
    end;
end;

for i=1:nry
    for j=1:nrx
        cunr=reshape（tu（1+8*（i-1）:8+8*（i-1）1+8*（j-1）:8+8*（j-1））[164]）;
        sumalpha=sum（cunr）;   %cunr  is ri
        sumalpha2=norm（cunr）^2;
        dx=1;
        dy=1;
        ut=1;
        minH=10^20;
        minot=0;
        minst=0;
        for k=1:ndy
            for l=1:ndx
                for m=1:8
                    cund3（1:64）=DD（klm1:64）;
                    sumbeta=sum（cund3）;  % cund3 is di
                    sumbeta2=norm（cund3）^2;
                    alphabeta=cunr*cund3‘;
                    if （64*sumbeta2-sumbeta^2）~=0
                    st=（64*alphabeta-sumalpha*sumbeta）/（64*sumbeta2-sumbeta^2）;
                    elseif （64*alphabeta-sumalpha*sumbeta）==0||st > 1 || st < -1
                        st=0;
                    else
                        st=10^20;
                    end;
                    ot=（sumalpha-st*sumbeta）/64;
                    H=（sumalpha2+st*（st*sumbeta2-2*alphabeta+2*