分形插值算法和MATLAB实验

一，分形插值算法

——分形图的递归算法1，分形的定义

分形（Fractal）一词，是法国人B.B.Mandelbrot 创造出来的，其原意包含了不规则、支离破碎等意思。Mandelbrot 基于对不规则的几何对象长期地、系统地研究，于1973 年提出了分维数和分形几何的设想。分形几何是一门以非规则几何形状为研究对象的几何学，用以描述自然界中普遍存在着的不规则对象。分形几何有其显明的特征，一是自相似性；分形作为一个数学集合, 其内部具有精细结构, 即在所有比例尺度上其组成部分应包含整体, 而且彼此是相似的。其定义有如下两种描述：

定义 1如果一个集合在欧式空间中的 Hausdorff 维数H D 恒大于其拓扑维数

r D ，则称该集合为分形集，简称分形。

定义 2组成部分以某种方式与整体相似的形体叫分形。

对于定义 1 的理解需要一定的数学基础，不仅要知道什么是Hausdorff 维数，而且要知道什么是拓扑维数，看起来很抽象，也不容易推广。定义 2 比较笼统的说明了自然界中的物质只要局部和局部或者局部和整体之间存在自相似性，那么这个物质就是分形。正是这一比较“模糊”的概念被人们普遍接受，同时也促进了分形的发展。

根据自相似性的程度，分形可分为有规分形和无规分形。有规分形是指具有严格的自相似的分形，比如，三分康托集，Koch 曲线。无规分形是指具有统计意义上的自相似性的分形，比如，曲折的海岸线，漂浮的云等。本文主要研究有规分形。

2. 分形图的递归算法

2.1 三分康托集

1883 年，德国数学家康托(G.Cantor)提出了如今广为人知的三分康托集。三分康托集是很容易构造的，然而，它却显示出许多最典型的分形特征。它是从单位区间出发，再由这个区间不断地去掉部分子区间的过程构造出来的(如图2.1)。

其详细构造过程是：第一步，把闭区间[0，1]平均分为三段，去掉中间的 1/3 部分段，则只剩下两个闭区间[0，1/3]和[2/3，1]。第二步，再将剩下的两个闭区间各自平均分为三段，同样去掉中间的区间段，这时剩下四段闭区间：[0，1/9]，[2/9，1/3]，[2/3，7/9]和[8/9，1]。第三步，重复删除每个小区间中间的 1/3 段。如此不断的分割下去，最后剩下的各个小区间段就构成了三分康托集。三分康托集的 Hausdorff 维数是0.6309。

图2.2 三分康托集的构造过程

2．2 Koch 曲线

1904 年，瑞典数学家柯赫构造了“Koch 曲线”几何图形。Koch 曲线大于一维，具有无限的长度，但是又小于二维，并且生成的图形的面积为零。它和三分康托集一样，是一个典型的分形。根据分形的次数不同，生成的Koch 曲线也有很多种，比如三次 Koch 曲线，四次 Koch 曲线等。下面以三次 Koch 曲线为例，介绍 Koch 曲线的构造方法，其它的可依此类推。

三次Koch 曲线的构造过程主要分为三大步骤：第一步，给定一个初始图形一条线段；第二步，将这条线段中间的 1/3 处向外折起；第三步，按照第二步的方法不断的把各段线段中间的 1/3 处向外折起。这样无限的进行下去，最后即可构造出Koch 曲线。其图例构造过程如下图2.2所示(迭代了 6 次的图形)。

图 2.2 Koch曲线生成过程

2.3 Sierpinski 垫片的递归算法

在 Cantor 集与Koch 集的基础上继续讨论Sierpinski 垫片的生成算法，同样先建立模型如下图2.3，同时展示生成后的图形。

图2.3 Sierpinski 垫片模型

其中先设定正三角形中心点坐标为(x, y)，边长为L，迭代深度为n，三角形顶点分别为( xi,yi ),i= 1,2,3。且下一层正三角形的中心坐标为（xi，yi），i=01,02,03。

2.4 分支结构分形递归算法

研究如下图（图 2.4）的分支结构图的递归算法。

图 2.4分支结构分形图

细分此分支结构，建立模型如下，其中取 A 为起点，且记 A 点坐标为(x，y)，B 点坐标为（x1，x2），线段AB = L, 2 .3BC = BD = L 且设定 AB 与水平面的夹角为alpha, 递归深度为n。

3 基于MATLAB的各种递归算法实例

3.1 三分康托集的MATLAB实例

Cantor 三分集的递归算法设计为构造一个函数为Cantor(ax,ay,bx,by)，

其中(ax,ay)和(bx,by)分别代表初始位置端点的坐标。取 err 为递归终止的极小量，取 h 为不同层次线段之间的距离，且设出(cx,cy)和(dx,dy)如图3.1

图 3.1 建立Cantor 三分集模型

（1），MATLAB程序

Cantor(ax,ay,bx,by)

{ plot([ax,bx],[ay,by]); %画出每个层次上的图形

if (bx-ax)>err%设定有限次递归

{ cx = ax + ( bx – ax ) / 3;

cy = ay – h;

dx = bx – ( bx – ax ) / 3;

dy = by – h;

ay = ay – h;

by = by – h;

Cantor(ax,ay,cx,cy); %内部调用Cantor函数Cantor(dx,dy,bx,by); %内部调用Cantor函数}

}

（2） ,生成结果图形

图 3.2 Cantor 三分集结果

3.2 Koch 曲线的MATLAB实例

在设计Koch曲线的递归算法时参考Cantor集的方法，先建立一个模型如下，

图 3.3 Koch曲线模型

（1），MATLAB程序

Koch(ax,ay,bx,by)

{ err=1; %递归终极小量

line([ax bx],[ay by])；hold on;%画出图形

if ((bx-ax)^2+(by-ay)^2)>err%递归终止条件

{ cx=ax+(bx-ax)/3;

cy=ay+(by-ay)/3;

ex=bx-(bx-ax)/3;

ey=by-(by-ay)/3;

L=sqrt((ex-cx)^2+(ey-cy)^2);%计算 L

alpha=atan((ey-cy)/(ex-cx)); %计算a

if((ex-cx)<0)

{

alpha=alpha+pi; %图形旋转

}

dy=cy+sin(alpha+pi/3)*L;

dx=cx+cos(alpha+pi/3)*L;

Koch(ax,ay,cx,cy); %递归调用

Koch(ex,ey,bx,by);

Koch(cx,cy,dx,dy);

Koch(dx,dy,ex,ey);

（2），生成结果图形

图 3.4 Koch曲线生成结果

3.3 Sierpinski 垫片的MATLAB实例

（1），MATLAB程序

Sierpinski(x, y, L, n)

{if (n==1)

x1=x-L/2;

y1=y-L/2*tan(pi/6);

x2=x+L/2;

y2=y-L/2*tan(pi/6);

x3=x;

y3=y+L/2*tan(pi/6);

line([x1 x2],[y1 y2],'Color','r','LineWidth',2); hold on;

line([x2 x3],[y2 y3],'Color','g','LineWidth',2); hold on;

line([x3 x1],[y3 y1],'Color','y','LineWidth',2); hold on;

else

x01=x-L/4;

y01=y-L/4*tan(pi/6);

x02=x+L/4;

y02=y-L/4*tan(pi/6);

x03=x;

y03=y+L/2*tan(pi/6);

Sierpinski(x01,y01,L/2,n-1); %递归调用Sierpinski(x02,y02,L/2,n-1);

Sierpinski(x03,y03,L/2,n-1);

end

}

（2），生成结果图形

图 3.5 n=6时Sierpinski垫片生成图

图3.6 n=10 时Sierpinski 垫片生成图

3.4 分支结构分形的MATLAB实例

（1），MATLAB程序

Ramus(x, y,alpha, L, n)

if n>0

x1=x+cos(alpha)*L;

y1=y+sin(alpha)*L;

line([x,x1],[y,y1],'Color','g','LineWidth',2);

alpha_L=alpha+pi/8;

alpha_R=alpha-pi/8;

L=2*L/3;

Ramus(x1,y1,alpha_L,L,n-1); %递归调用

Ramus(x1,y1,alpha_R,L,n-1);

End

（2），生成结果图形

图 3.7α=π/4时分支结构生成图

图 3.7α=π/3时分支结构生成图

二， MATLAB实验

1.图像的二值化

程序如下：

I=imread('C:\Users\Administrator.SD-20111228XIYE\Desktop\001.j pg');

I1=rgb2gray(I);

level=graythresh(I1);

I2=im2bw(I1,level);

I3=~I2;

I4=bwareaopen(I3,50);

I5=~I4;

figure,imshow(I5)

原图：

二值化后的图像：

2.三种插值算法

程序：

I=imread('C:\Users\Administrator.SD-20111228XIYE\Desktop\001.jp g ');

figure, imshow(I);

A=imresize(I,2,'nearest'); %最邻近插值

figure, imshow(A);

B=imresize(I,2,'bilinear');%双线性插值

figure, imshow(B);

C=imresize(I,2,'bicubic'); %三次插值

figure, imshow(C);

（1）最邻近插值

（2）双线性插值

（3）立方卷积插值

3.抠图实验

MATLAB程序：

Main

close all

clear all

clc

img_name = input('请输入图像名字(图像必须为RGB图像，输入0结束)：','s'); % 当输入0时结束

while ~strcmp(img_name,'0')

%进行人脸识别

facedetection(img_name);

img_name = input('请输入图像名(图像必须为RGB图像，输入0结束)：','s'); end

子程序1 RGB

function rgbPic = bw2rgb(bwPic)

bwPicSize = size(bwPic);

rgbPic = zeros(bwPicSize(1),bwPicSize(2),3); rgbPic(bwPic==255)=255;

rgbPic(:,:,2) = rgbPic(:,:,1); % nice job

rgbPic(:,:,3) = rgbPic(:,:,1);

rgbPic = im2uint8(rgbPic);

子程序2 人脸识别function facedetection(img_name)

I = imread(img_name);

I2 = imread('dog.png');

doge = imresize(I2,[size(I,1) size(I,2)]);

gray = rgb2gray(I);

YCbCr = rgb2ycbcr(I);

heigth = size(gray,1);

width = size(gray,2);

for i = 1:heigth

for j = 1:width

Y = YCbCr(i,j,1);

Cb = YCbCr(i,j,2);

Cr = YCbCr(i,j,3);

if(Y < 80)

gray(i,j) = 0;

else

if(skin(Y,Cb,Cr) == 1)

gray(i,j) = 255;

else

gray(i,j) = 0;

end

end

end

end

SE=strel('arbitrary',eye(5));

gray = imopen(gray,SE);

gray = imclose(gray,SE);

a1=(255-bw2rgb(gray))/255;

next1=immultiply(I,a1);

a2=bw2rgb(gray)/255;

next2=immultiply(doge,a2);

add=imadd(next1,next2);

[L,num] = bwlabel(gray,8);

STATS = regionprops(L,'BoundingBox');

n = 1;

result = zeros(n,4);

figure,

subplot(1,2,1),imshow(I),

subplot(1,2,2),imshow(doge);

figure,

subplot(1,2,1),imshow(next1),

subplot(1,2,2),imshow(next2);

figure,

imshow(add);

hold on;

for i = 1:num

box = STATS(i).BoundingBox;

x = box(1);

y = box(2);

w = box(3);

h = box(4);

ratio = h/w;

ux = uint8(x);

uy = uint8(y);

if ux > 1

ux = ux - 1;

end

if uy > 1

uy = uy - 1;

end

if w < 20 || h < 20 || w*h < 400

continue

elseif ratio < 2.0 && ratio > 0.6 && findeye(gray,ux,uy,w,h) == 1

result(n,:) = [ux uy w h];

n = n+1;

end

end

if size(result,1) == 1 && result(1,1) > 0

rectangle('Position',[result(1,1),result(1,2),result(1,3),result(1,4)],'EdgeColor','r');

for m = 1:size(result,1)

m1 = result(m,1);

m2 = result(m,2);

m3 = result(m,3);

m4 = result(m,4);

if m1 + m3 < width && m2 + m4 < heigth

rectangle('Position',[m1,m2,m3,m4],'EdgeColor','r');

end

end

end

子程序3 find eye

function eye = findeye(bImage,x,y,w,h)

part = zeros(h,w);

for i = y:(y+h)

for j = x:(x+w)

if bImage(i,j) == 0

part(i-y+1,j-x+1) = 255;

else

part(i-y+1,j-x+1) = 0;

end

end

end

[L,num] = bwlabel(part,8);

if num < 2

eye = 0;

else

eye = 1;

end

子程序4 skin

function result = skin(Y,Cb,Cr)

a = 25.39;

b = 14.03;

ecx = 1.60;

ecy = 2.41;

sita = 2.53;

cx = 109.38;

cy = 152.02;

xishu = [cos(sita) sin(sita);-sin(sita) cos(sita)];