求迭代分形的MATLAB程序

function

cycle6(px,py,r,n)

%px,py和r分别是大圆圆心坐标和半径,n是迭代次数

%例子:

%cycle6(0,0,1,4)axis

equal

t=linspace(-pi,pi)

plot(r*cos(t)+px,r*sin(t)+py)

theta=pi/3*(0:5)

newr=r/3

newpx=px+2*newr*cos(theta)

newpy=py+2*newr*sin(theta)

if

n>1

for

i=1:6

hold

oncycle6(newpx(i),newpy(i),newr,n-1)

end

end

代码有很多小错误，我帮你修改了下，

这是函数文件

function S1tree(n)

clc

S='F'a=pi/10A=pi/2z=0zA=[0,pi/2]

p1='FF+[+F+F]-[+F]'

p2='F[+F]F[-F[+F]]'

p3='FF-[-F+F+F]+[+F-F-F]'

for k=2:n

  c=rand(1)

  if c>=0.7      S=strrep(S,'F',p1)

  elseif c>=0.35 S=strrep(S,'F',p2)

  else           S=strrep(S,'F',p3)

  end

end

figure

for k=1:length(S)

  switch S(k)

      case 'F'

          plot(real(z+2*exp(i*A)),imag(z+2*exp(i*A)),'g','LineWidth',2)

          hold on

          z=z+2*exp(i*A)

      case '+'

          A=A+a

      case '-'

          A=A-a

      case '['

          zA=[zA[z,A]]

      case ']'

          z=zA(end,1)

          A=zA(end,2)

          zA(end,:)=[]

      otherwise

  end

end

在主窗口中输入

S1tree(7)

画出的图如下（由于每次运行S1tree(7)代码产生随机数不一样，得到的图不一样但是类似）

分形编码的基本思路是：先采用一种合适的初级压缩方法对图像进行压缩，得到一组压缩编码，然后解码这组编码，得到一幅解压缩图像．对解码图像与原始图像求差值，得到一差值图像，然后对该差值图像进行适当的编码．对差值图像的编码与初级编码共同构成对原始图像的编码．这种方法需要选择合适的初级编码方法与差值编码方法，使得这两者相结合，可以得到一种综合性能较好的编码方法。

其定义如图。

%%%%%%%%%%%%%%

clear

tic

%Image1=imread('pic\cameraman.tif')

xianshi

number=input(' input the number:')

Image1=suoxiao('pic\cameraman.tif',number)

[imagem imagen]=size(Image1)

Sr=4Sd=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))')

%压缩image1为原来1/2

% there are no  eight tranformation for simpleness

RBlocks=zeros(Rnum,Sr,Sr)

DBlocks=zeros(Dnum,Sd,Sd)

DBlocksReduce=zeros(Dnum*8,Sr,Sr)

%%取R块，K记标号----------------------------------

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

%取R块，K记标号----------------------------------

for i=1:imagem/Sd

for j=1:imagen/Sd

k=(i-1)*imagen/Sd+j

m=Srn=Sr

DBlocksReduce(k,:,:)=Image2((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr)

DBlocksReduce(k+Dnum,:,:)=DBlocksReduce(k,m:-1:1,:)              % 行上下翻转===(x轴对称)

DBlocksReduce(k+2*Dnum,:,:)=DBlocksReduce(k,:,n:-1:1)              % 列左右翻转 ==== y轴对称

DBlocksReduce(k+3*Dnum,:,:)=DBlocksReduce(k,m:-1:1,n:-1:1)         % 先行翻，再列翻  旋转180度

DBlocksReduce(k+4*Dnum,:,:)=reshape(DBlocksReduce(k,:,:),Sr,Sr)'                       % 关于y=-x对称

A=reshape( DBlocksReduce(k+3*Dnum,:,:),Sr,Sr)'

DBlocksReduce(k+5*Dnum,:,:)=A(:,n:-1:1)     % 关于y=x对称

DBlocksReduce(k+6*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),90)           % 逆时针旋转90度

DBlocksReduce(k+7*Dnum,:,:)=imrotate(reshape(DBlocksReduce(k,:,:),Sr,Sr),270)          % 逆时针旋转270度

DBlocks(k,:,:)=Image1((i-1)*Sd+1:i*Sd,(j-1)*Sd+1:j*Sd)

end

end

RandDbest=zeros(Rnum,1)+256^3

RandDbests=zeros(Rnum,1)

RandDbesto=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)

s=(x-meanx)'*(y-meany)/((y-meany)'*(y-meany))%计算s

o=(meanx-s*meany)%计算o

c=(x-s*y-o)'*(x-s*y-o)%距离

if (RandDbest(i)>c)&(abs(s)<1)

RandDbest(i)=c

RandDbests(i)=s

RandDbesto(i)=o

RandDbestj(i)=j%可以找到对应变换和D块

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/Sr

for j=1:imagen/Sr

m=Srn=Sr

k=(i-1)*imagen/Sr+j

l=RandDbestj(k)

k1=mod(l-1,Dnum)+1%第几个D

l1=(l-k1)/Dnum+1%变换号

%R对应D在Image4的起始点

j1=mod(k1-1,imagen/Sd)+1

i1=(k1-j1)/(imagen/Sd)+1

%变换------------------------------------------------------------------------

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,:,:)=RandDbests(k)*DBlocksReduce(l,:,:)+RandDbesto(k)

%生成R---------------------------

Image3((i-1)*Sr+1:i*Sr,(j-1)*Sr+1:j*Sr)=reshape(RBlocks(k,:,:),Sr,Sr)%更新迭代图象

end

end

wucha=double(Image1)-Image3%误差图

Ps1(L)=20*log10(255/(sqrt(mean(mean(wucha.^2)))))

PSNR=psnr(wucha)

figure

imshow(uint8(Image3))

end

toc

figure

wucha=uint8(wucha)

imshow(wucha)

figure

imshow(uint8(Image1)),title('原图')

save('sa.mat')

zhifangtu(wucha)%%%%分形主函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

子函数:1：

function b=suoxiao(filename,bili)

a=imread(filename)

a=double(a)

[m,n]=size(a)

i=1

while i=m/bili

j=1

while j=n/bili

k=mean(mean(a(bili*(i-1)+1:bili*(i-1)+bili,bili*(j-1)+1:bili*(j-1)+bili)))

b(i,j)=k

j=j+1

end

i=i+1

end

%b=uint8(b)

size(b)

%imshow(b)

子函数2：

%clc

function zhifangtu(a)

J=a

%计算灰度图象的直方图数据，a为如象数组

L=256 %灰度级

Ps = zeros(L,1) %统计直方图结果数据

nk=zeros(L,1)

[row,col]=size(a)

n=row*col %总像素个数

for i = 1:row

for j = 1:col

num = double(a(i,j))+1 %获取像素点灰度级

nk(num) = nk(num)+1 %统计nk

end

end

%计算直方图概率估计

for i=1:L

Ps(i)=nk(i)/n

end

figure

subplot(3,1,1)imshow(J),title('误差图')

subplot(3,1,2),plot(nk),title('直方图（nk）')

subplot(3,1,3),plot(Ps),title('直方图（Ps）')

子函数3：

function PSNR=psnr(a)

[m,n]=size(a)

a=uint8(a)

a=double(a)

imagesize=m*n

MSE=sum(dot(a,a))/ imagesize

PSNR=10*log10(255^2/MSE)

%%%%%%%%%%%%%%%%%%%%%%

说明：

1、因为本程序时间长，FX中先选择图片的大小

2、编码与解码

3、做误差图和只方图

4：画出每次迭代的解码图象

---