MATLAB程序大全

1、.1全景图到穹景图这个程序我最初是用FreeImage写的，这两天改成了matlab，再不贴上来，我就要忘了。看到一篇文章有这样的变换，挺有意思的，就拿来试了一下，文章点此。全景图到穹顶图变换，通俗的说就是将全景图首尾相接做成一个圆环的样子。先看下面这张图：下面的矩形就是我们要处理的全景图，上面的矩形是变换后的图像。下面图像的底边对应穹顶图的内圆，顶边对应穹顶图的外圆，当然，反过来也是可以的。程序流程：1.定义穹顶图内圆和外圆的半径，变换后的像素就填充在这个内外半径的圆环中。2.遍历穹顶图，当所处理当前像素位于圆环内，则通过极坐标反变换去全景图中寻找相应位置的像素进行填充。3.遍历完图像就行了

2、。用的技巧和图像旋转或放大缩小都是类似的。处理结果：原图：结果：matlab代码如下：clear all;close all;clc;img=imread(pan.jpg);imshow(img);m,n=size(img);r1=100; %内环半径r2=r1+m; %外环半径imgn=zeros(2*r2,2*r2);re_m,re_n=size(imgn);for y=1:re_m for x=1:re_n dis_x=x-re_n/2; dis_y=y-re_m/2; l=sqrt(dis_x2+dis_y2); if l=r1 theta=0; if yre_m/2 theta=at

3、an2(dis_y,dis_x); end if y=1 & yy=1 & xx0 & K=0;角点区域：H0 & K0;harris角点检测就用到了第三类判断。当然，在实际应用的时候H和K的值肯定都不会是理想，所以我用的都是近似判断。处理结果如下：原图：平坦区域：边缘区域：角点区域（好像也不全角点，求角点还是harris好了）：结构张量行列式与迹的关系：其中红框为平坦区域，黄框为边缘区域，铝框为角点区域。matlab代码如下：clear all; close all; clc;img=double(imread(lena.jpg);m n=size(img);imshow(img,)Ix I

4、y=gradient(img);Ix2=Ix.2;Iy2=Iy.2;Ixy=Ix.*Iy;k=1;lambda=zeros(m*n,2);for i=1:m for j=1:n st=Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j); %结构张量 K=det(st); %求行列式 H=trace(st); %求迹 %所有的判断都是近似的 % if H50 & abs(K)50 & abs(K)0.01*10(-9) %认为是角点区域 img(i,j)=255; end lambda(k,:)=K H; k=k+1; endendfigure;plot(lambda(:,1

5、),lambda(:,2),.);ylabel(trace);xlabel(det);figure;imshow(img,)6.模糊集图像增强算法有很多变种。不过主要就是以下三步。1.设计隶属度函数将图像从空间域变换到模糊集域。2.设计模糊增强算子，在模糊集域对图像进行处理。3.根据第1步的隶属度函数重新将图像从模糊集域变换到空间域。这和频域处理中的变换反变换不是很像么。我使用的隶属度函数和模糊增强算子在这篇论文里，也算相关算法的经典论文了。处理结果如下：原图：模糊集增强后：matlab代码如下：clear all; close all; clc;img=double(imread(lena.

6、jpg);imshow(img,)m n=size(img);Fe=1;%控制参数Fd=128;xmax=max(max(img);u=(1+(xmax-img)/Fd).(-Fe); %空间域变换到模糊域%也可以多次迭代for i=1:m %模糊域增强算子 for j=1:n if u(i,j)0 & isempty(find(flag=j,1) queue(head)=j; head=head+1; flag=flag j; pa(j)=i; end end tail=tail+1; end if pa(m)=0 %如果搜索不到汇节点，退出循环 break; end path=; i=m;

7、 %从汇节点开始 k=0; %路径包含的边的个数 while i=1 %使用前趋构造从源节点到汇节点的路径 path=path;pa(i) i A(pa(i),i); %存入路径 i=pa(i); %使用前趋表反向搜寻，借鉴Dijsktra中的松弛方法 k=k+1; end Mi=min(path(:,3); %寻找增广路径中最小的那条边 for i=1:k A(path(i,1),path(i,2)=A(path(i,1),path(i,2)-Mi; %增广路径中每条路径减去最小的边 maxflow(path(i,1),path(i,2)=maxflow(path(i,1),path(i,2

8、)+Mi; %最大流，原网络包含这个网络，我只能这样表示了 end %使用新的图A进入下一循环，从新开始找增广路径end figure;netplot(maxflow,1)9.单元最短路径图的相关算法也算是自己的一个软肋了，当年没选修图论也是一大遗憾。图像处理中，也有使用图论算法作为基础的相关算法，比如图割，这个算法就需要求最大流、最小割。所以熟悉一下图论算法对于图像处理还是很有帮助的。Dijkstra和Bellman-Ford类似，都是解决单源最短路径问题，不同的是这个方法只能解决边为非负的问题，实现的好的Dijkstra算法运行时间要快于Bellman-ford。算法步骤如下：1.首先设置

9、队列，所有节点入列，源节点值为0，其他节点值为无穷。2.然后在队列中找值最小的节点并出列。3.计算出列的节点所有后继节点的距离。4.松弛方法，如果新计算的距离小于上次计算的距离，那么更新距离，即将后继节点值设为较小的距离，并将后继节点的前趋设为当前的出列节点。5.对剩余的节点队列继续找最小值并出列，不断循环2、3、4步直到队列中没有节点了。步骤是上面没错，不过我程序中没有完全按照上述的步骤实现。不同的地方在于我没有做出列操作，而是通过标记节点的形式实现的。运行结果如下，图（是图不是图片）是算法导论367页上的：matlab代码如下，netplot和compresstable2matrix和上一

10、篇使用的一样：main.mclear all;close all;clc%初始化邻接压缩表，1 2 10 表示从节点1到节点2，边的权重为10b=1 2 10;1 4 5;2 3 1; 2 4 2; 3 5 4;4 2 3; 4 3 9; 4 5 2;5 1 7; 5 3 6;m=max(max(b(:,1:2); %压缩表中最大值就是邻接矩阵的宽与高A=compresstable2matrix(b); %从邻接压缩表构造图的矩阵表示netplot(A,1) %形象表示S=inf(1,m); %从开始的源点到每一个节点的距离S(1)=0; %源点到自己的距离为0pa=zeros(1,m); %

11、存储每个节点的前驱，在松弛过程中赋值pa(1)=1;%源点的前趋是自己visit=zeros(1,m); %标记某个节点是否访问过了index=1; %从index节点开始搜索 %判断是否对所有节点都找的最短路径了。可能会有源点没有路径到目标节点的情况，那就无限循环了while sum(visit)=m %没有出队列操作，不过通过visit来等价的表示了 visit(index)=1; %标记第index节点为已入列的节点 S pa=relax(S,pa,A,visit,index,m); %松弛，如果两个节点间有更短的距离，则用更短的距离 index=extract_min(S,visit,

12、index,m); %使用已访问的最小的节点作为下一次搜索的开始节点 end%最终我们需要的就是这两个值S %源点到其他每一点的距离pa %其他每一节点的前趋%算法到此结束，下面只是为了形象的表示而写的。re=;for i=2:m re=re;pa(i) i A(pa(i),i);endA=compresstable2matrix(re); %从邻接压缩表构造图的矩阵表示figure;netplot(A,1) %形象表示relax.m%边缘松弛，使用更短的距离作为节点的值function S pa=relax(S,pa,A,visit,index,m) i=index; for j=1:m i

13、f A(i,j)=inf & visit(j)=1 %搜索没有标记过的节点 if S(j)S(i)+A(i,j) %将较小的值赋给正在搜寻的节点 S(j)=S(i)+A(i,j); pa(j)=i; end end end endextract_min.m%提取队列中尚未标记的最小的值的序号function index=extract_min(S,visit,index,m) Mi=inf; for j=1:m if visit(j)=1 if S(j)S(i)+A(i,j) S(j)=S(i)+A(i,j); %边缘松弛，取两节点间最小权值作为实际权值 pa(j)=i; %寻找前趋 end

14、end end end endend%最终我们需要的就是这两个值S %源点到其他每一点的距离pa %其他每一节点的前趋%算法到此结束，下面只是为了形象的表示而写的。re=;for i=2:m re=re;pa(i) i A(pa(i),i);endA=compresstable2matrix(re); %从邻接压缩表构造图的矩阵表示figure;netplot(A,1) %形象表示compresstable2matrix.mfunction A=compresstable2matrix(b) n =size(b); m=max(max(b(:,1:2); A=inf(m,m); for i=1

15、:n A(b(i,1),b(i,2)=b(i,3); endend11.如此经典的算法竟一直没有单独的实现过，真是遗憾啊。广度优先搜索在过去实现的二值图像连通区域标记和prim最小生成树算法时已经无意识的用到了，深度优先搜索倒是没用过。这次单独的将两个算法实现出来，因为算法本身和图像没什么关系，所以更纯粹些。广度优先搜索是从某一节点开始，搜索与其线连接的所有节点，按照广度方向像外扩展，直到不重复遍历所有节点。深度优先搜索是从某一节点开始，沿着其搜索到的第一个节点不断深入下去，当无法再深入的时候，回溯节点，然后再在回溯中的某一节点开始沿另一个方向深度搜索，直到不重复的遍历所有节点。广度优先搜索用

16、的是队列作为临时节点存放处；深度优先搜索可以递归实现（算法导论就是用递归实现的伪代码），不过我这里是用栈作为临时节点存放处。感觉也没什么好介绍的了，抄算法导论上的介绍也没什么意思，所有的内容都是书上的，真正学东西还是要看书。下面是运行结果：原连通图：广度优先搜索：深度优先搜索：matlab代码如下，其中的画图函数netplot.m。BFS.mclear all;close all;clc%初始化邻接压缩表b=1 2;1 3;1 4;2 4; 2 5;3 6;4 6;4 7;m=max(b(:); %压缩表中最大值就是邻接矩阵的宽与高A=compresstable2matrix(b); %从邻接

17、压缩表构造图的矩阵表示netplot(A,1) %形象表示head=1; %队列头tail=1; %队列尾，开始队列为空，tail=headqueue(head)=1; %向头中加入图第一个节点head=head+1; %队列扩展flag=1; %标记某个节点是否访问过了re=; %最终结果while tail=head %判断队列是否为空 i=queue(tail); %取队尾节点 for j=1:m if A(i,j)=1 & isempty(find(flag=j,1) %如果节点相连并且没有访问过 queue(head)=j; %新节点入列 head=head+1; %扩展队列 fla

18、g=flag j; %对新节点进行标记 re=re;i j; %将边存入结果 end end tail=tail+1; endA=compresstable2matrix(re);figure;netplot(A,1)DFS.mclear all;close all;clc%初始化邻接压缩表b=1 2;1 3;1 4;2 4; 2 5;3 6;4 6;4 7;m=max(b(:); %压缩表中最大值就是邻接矩阵的宽与高A=compresstable2matrix(b); %从邻接压缩表构造图的矩阵表示netplot(A,1) %形象表示top=1; %堆栈顶stack(top)=1; %将第一

19、个节点入栈flag=1; %标记某个节点是否访问过了re=; %最终结果while top=0 %判断堆栈是否为空 pre_len=length(stack); %搜寻下一个节点前的堆栈长度 i=stack(top); %取堆栈顶节点 for j=1:m if A(i,j)=1 & isempty(find(flag=j,1) %如果节点相连并且没有访问过 top=top+1; %扩展堆栈 stack(top)=j; %新节点入栈 flag=flag j; %对新节点进行标记 re=re;i j; %将边存入结果 break; end end if length(stack)=pre_len

20、%如果堆栈长度没有增加，则节点开始出栈 stack(top)=; top=top-1; end endA=compresstable2matrix(re);figure;netplot(A,1)compresstable2matrix.mfunction A=compresstable2matrix(b) n =size(b); m=max(b(:); A=zeros(m,m); for i=1:n A(b(i,1),b(i,2)=1; A(b(i,2),b(i,1)=1; endend12.模拟退火首先从某个初始候选解开始，当温度大于0时执行循环。在循环中，通过随机扰动产生一个新的解，然后求

21、得新解和原解之间的能量差，如果差小于0，则采用新解作为当前解。如果差大于0，则采用一个当前温度与能量差成比例的概率来选择是否接受新解。温度越低，接受的概率越小，差值越大，同样接受概率越小。是否接受的概率用此公式计算：p=exp(-E/T)。这里E为新解与原解的差，T为当前的温度。由于温度随迭代次数逐渐降低，因此获得一个较差的解的概率较小。典型的模拟退火算法还使用了蒙特卡洛循环，在温度降低之前，通过多次迭代来找到当前温度下比较好的解。这里使用模拟退火解旅行商问题，因为这个问题本身是一个NP难问题，所以也就求不到最优解，不过应该可以求得一个比较好的解，然后再手工优化。具体到程序中，这里的随机扰动就

22、是随机置换两个城市的位置，能量就是旅行商路线的总长度。算法结果如下：初始旅行商路线：最终求得的旅行商路线：每次迭代求得的旅行距离：matlab代码如下：main.mclear all;close all;clcn=20; %城市个数temperature=100*n; %初始温度iter=100; %内部蒙特卡洛循环迭代次数%随机初始化城市坐标city=struct();for i=1:n city(i).x=floor(1+100*rand(); city(i).y=floor(1+100*rand();endl=1; %统计迭代次数len(l)=computer_tour(city,n);

23、 %每次迭代后的路线长度 netplot(city,n); %初始旅行路线while temperature0.001 %停止迭代温度 for i=1:iter %多次迭代扰动，一种蒙特卡洛方法，温度降低之前多次实验 len1=computer_tour(city,n); %计算原路线总距离 tmp_city=perturb_tour(city,n); %产生随机扰动 len2=computer_tour(tmp_city,n); %计算新路线总距离 delta_e=len2-len1; %新老距离的差值，相当于能量 if delta_erand() %以概率选择是否接受新解 city=tmp

24、_city; %可能得到较差的解 end end end l=l+1; len(l)=computer_tour(city,n); %计算新路线距离 temperature=temperature*0.99; %温度不断下降 end figure;netplot(city,n); %最终旅行路线figure;plot(len) computer_tour.mfunction len=computer_tour(city,n) %计算路线总长度，每个城市只计算和下家城市之间的距离。 len=0; for i=1:n-1 len=len+sqrt(city(i).x-city(i+1).x)2+(

25、city(i).y-city(i+1).y)2); end len=len+sqrt(city(n).x-city(1).x)2+(city(n).y-city(1).y)2);endperturb_tour.mfunction city=perturb_tour(city,n) %随机置换两个不同的城市的坐标 %产生随机扰动 p1=floor(1+n*rand(); p2=floor(1+n*rand(); while p1=p2 p1=floor(1+n*rand(); p2=floor(1+n*rand(); end tmp=city(p1); city(p1)=city(p2); ci

26、ty(p2)=tmp;endnetplot.mfunction netplot(city,n) %连线各城市，将路线画出来 hold on; for i=1:n-1 plot(city(i).x,city(i).y,r*); line(city(i).x city(i+1).x,city(i).y city(i+1).y); %只连线当前城市和下家城市 end plot(city(n).x,city(n).y,r*); line(city(n).x city(1).x,city(n).y city(1).y); %最后一家城市连线第一家城市 hold off;end13.还是这本书上的内容，不

27、过我看演化计算这一章是倒着看的，这里练习的算法正好和书中介绍的顺序是相反的。演化策略是最古老的的演化算法之一，和上一篇DE算法类似，都是基于种群的随机演化产生最优解的算法。算法步骤如下：1.设定种群个体数和需要迭代的次数。2.选择父代中的个体按照公式z1=sqrt(-2*ln(u1)*sin(2*pi*u2)*m，z2=sqrt(-2*ln(u1)*cos(2*pi*u2)*m进行演化。这里u1,u2都是随机值，m是控制因子，演化次数越多m，m越小，父代通过与z1，z2相加得到后代。3.计算后代的适应性。4.选择后代中最优的适应性作为全局最优适应性。其实整个过程和DE非常类似。过程都是随机变异

28、，求适应性，再找最优。我还试着将z1和z2横设为1，竟也能得到非常好的解。算法结果如下：matlab代码如下：main.mclear all;close all;clc;x y=meshgrid(-100:100,-100:100);sigma=50;img = (1/(2*pi*sigma2)*exp(-(x.2+y.2)/(2*sigma2); %目标函数，高斯函数mesh(img);hold on;n=50; %种群个体的数量iter=100; %迭代次数%初始化种群，定义结构体par=struct();for i=1:n par(i).x=-100+200*rand(); %个体的x特

29、征在-100 100随机初始化 par(i).y=-100+200*rand(); %个体的y特征在-100 100随机初始化 par(i).fit=compute_fit(par(i); %个体在x,y处的适应度endpar_best=par(1); %初始化种群中最佳个体for k=1:iter %迭代次数 plot3(par_best.x+100,par_best.y+100,par_best.fit,g*); %画出最佳个体的位置，+100为相对偏移 par par_best=select_and_recombin(par,par_best,n,k,iter); %差异演化函数ends

30、elect_and_recombin.mfunction next_par par_best=select_and_recombin(par,par_best,n,k,iter) mul=(iter-k)/iter; %限制进化因子，代数越高变异越小 next_par=par; %新种群 for i=1:n %产生变异随机数 u1=rand(); u2=rand(); z1=sqrt(-2*log(u1)*sin(2*pi*u2)*mul; z2=sqrt(-2*log(u1)*cos(2*pi*u2)*mul; %变异 next_par(i).x=par(i).x+z1; next_par(i).y=par(i).y+z2; %计算变异后个体的适应度 next_par(i).fit=compute_fit(next_par(i); %如果新个体没有变异前个体适应度高，新个体还原为旧个体 if par(i).fitnext_par(i).fit next_par(i)=par(i); end %如果变异后适应度高于种群最高适应个体，则更新种群适应度最高个体 if next_par(i).fitpar_best.fit par_best=next_par(i); end end endcompute_fit.mfunction re=compute_fit(p