MATLAB程序大全.doc

1全景圖到穹景圖這個(gè)程序我最初是用FreeImage寫(xiě)的，這兩天改成了matlab，再不貼上來(lái)，我就要忘了?？吹揭黄恼掠羞@樣的變換，挺有意思的，就拿來(lái)試了一下，文章點(diǎn)此。全景圖到穹頂圖變換，通俗的說(shuō)就是將全景圖首尾相接做成一個(gè)圓環(huán)的樣子。先看下面這張圖：下面的矩形就是我們要處理的全景圖，上面的矩形是變換后的圖像。下面圖像的底邊對(duì)應(yīng)穹頂圖的內(nèi)圓，頂邊對(duì)應(yīng)穹頂圖的外圓，當(dāng)然，反過(guò)來(lái)也是可以的。程序流程：1.定義穹頂圖內(nèi)圓和外圓的半徑，變換后的像素就填充在這個(gè)內(nèi)外半徑的圓環(huán)中。2.遍歷穹頂圖，當(dāng)所處理當(dāng)前像素位于圓環(huán)內(nèi)，則通過(guò)極坐標(biāo)反變換去全景圖中尋找相應(yīng)位置的像素進(jìn)行填充。3.遍歷完圖像就行了。用的技巧和圖像旋轉(zhuǎn)或放大縮小都是類似的。處理結(jié)果：原圖：結(jié)果：matlab代碼如下：clear all;close all;clc;img=imread(pan.jpg);imshow(img);m,n=size(img);r1=100; %內(nèi)環(huán)半徑r2=r1+m; %外環(huán)半徑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=atan2(dis_y,dis_x); end if y=1 & yy=1 & xx0 & K=0;角點(diǎn)區(qū)域：H0 & K0;harris角點(diǎn)檢測(cè)就用到了第三類判斷。當(dāng)然，在實(shí)際應(yīng)用的時(shí)候H和K的值肯定都不會(huì)是理想，所以我用的都是近似判斷。處理結(jié)果如下：原圖：平坦區(qū)域：邊緣區(qū)域：角點(diǎn)區(qū)域（好像也不全角點(diǎn)，求角點(diǎn)還是harris好了）：結(jié)構(gòu)張量行列式與跡的關(guān)系：其中紅框?yàn)槠教箙^(qū)域，黃框?yàn)檫吘墔^(qū)域，鋁框?yàn)榻屈c(diǎn)區(qū)域。matlab代碼如下：clear all; close all; clc;img=double(imread(lena.jpg);m n=size(img);imshow(img,)Ix Iy=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); %結(jié)構(gòu)張量 K=det(st); %求行列式 H=trace(st); %求跡 %所有的判斷都是近似的 % if H50 & abs(K)50 & abs(K)0.01*10(-9) %認(rèn)為是角點(diǎn)區(qū)域 img(i,j)=255; end lambda(k,:)=K H; k=k+1; endendfigure;plot(lambda(:,1),lambda(:,2),.);ylabel(trace);xlabel(det);figure;imshow(img,)6.模糊集圖像增強(qiáng)算法有很多變種。不過(guò)主要就是以下三步。1.設(shè)計(jì)隸屬度函數(shù)將圖像從空間域變換到模糊集域。2.設(shè)計(jì)模糊增強(qiáng)算子，在模糊集域?qū)D像進(jìn)行處理。3.根據(jù)第1步的隸屬度函數(shù)重新將圖像從模糊集域變換到空間域。這和頻域處理中的變換反變換不是很像么。我使用的隸屬度函數(shù)和模糊增強(qiáng)算子在這篇論文里，也算相關(guān)算法的經(jīng)典論文了。處理結(jié)果如下：原圖：模糊集增強(qiáng)后：matlab代碼如下：clear all; close all; clc;img=double(imread(lena.jpg);imshow(img,)m n=size(img);Fe=1;%控制參數(shù)Fd=128;xmax=max(max(img);u=(1+(xmax-img)/Fd).(-Fe); %空間域變換到模糊域%也可以多次迭代for i=1:m %模糊域增強(qiáng)算子 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 %如果搜索不到匯節(jié)點(diǎn)，退出循環(huán) break; end path=; i=m; %從匯節(jié)點(diǎn)開(kāi)始 k=0; %路徑包含的邊的個(gè)數(shù) while i=1 %使用前趨構(gòu)造從源節(jié)點(diǎn)到匯節(jié)點(diǎn)的路徑 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)+Mi; %最大流，原網(wǎng)絡(luò)包含這個(gè)網(wǎng)絡(luò)，我只能這樣表示了 end %使用新的圖A進(jìn)入下一循環(huán)，從新開(kāi)始找增廣路徑end figure;netplot(maxflow,1)9.單元最短路徑圖的相關(guān)算法也算是自己的一個(gè)軟肋了，當(dāng)年沒(méi)選修圖論也是一大遺憾。圖像處理中，也有使用圖論算法作為基礎(chǔ)的相關(guān)算法，比如圖割，這個(gè)算法就需要求最大流、最小割。所以熟悉一下圖論算法對(duì)于圖像處理還是很有幫助的。Dijkstra和Bellman-Ford類似，都是解決單源最短路徑問(wèn)題，不同的是這個(gè)方法只能解決邊為非負(fù)的問(wèn)題，實(shí)現(xiàn)的好的Dijkstra算法運(yùn)行時(shí)間要快于Bellman-ford。算法步驟如下：1.首先設(shè)置隊(duì)列，所有節(jié)點(diǎn)入列，源節(jié)點(diǎn)值為0，其他節(jié)點(diǎn)值為無(wú)窮。2.然后在隊(duì)列中找值最小的節(jié)點(diǎn)并出列。3.計(jì)算出列的節(jié)點(diǎn)所有后繼節(jié)點(diǎn)的距離。4.松弛方法，如果新計(jì)算的距離小于上次計(jì)算的距離，那么更新距離，即將后繼節(jié)點(diǎn)值設(shè)為較小的距離，并將后繼節(jié)點(diǎn)的前趨設(shè)為當(dāng)前的出列節(jié)點(diǎn)。5.對(duì)剩余的節(jié)點(diǎn)隊(duì)列繼續(xù)找最小值并出列，不斷循環(huán)2、3、4步直到隊(duì)列中沒(méi)有節(jié)點(diǎn)了。步驟是上面沒(méi)錯(cuò)，不過(guò)我程序中沒(méi)有完全按照上述的步驟實(shí)現(xiàn)。不同的地方在于我沒(méi)有做出列操作，而是通過(guò)標(biāo)記節(jié)點(diǎn)的形式實(shí)現(xiàn)的。運(yùn)行結(jié)果如下，圖（是圖不是圖片）是算法導(dǎo)論367頁(yè)上的：matlab代碼如下，netplot和compresstable2matrix和上一篇使用的一樣：main.mclear all;close all;clc%初始化鄰接壓縮表，1 2 10 表示從節(jié)點(diǎn)1到節(jié)點(diǎn)2，邊的權(quán)重為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); %從鄰接壓縮表構(gòu)造圖的矩陣表示netplot(A,1) %形象表示S=inf(1,m); %從開(kāi)始的源點(diǎn)到每一個(gè)節(jié)點(diǎn)的距離S(1)=0; %源點(diǎn)到自己的距離為0pa=zeros(1,m); %存儲(chǔ)每個(gè)節(jié)點(diǎn)的前驅(qū)，在松弛過(guò)程中賦值pa(1)=1;%源點(diǎn)的前趨是自己visit=zeros(1,m); %標(biāo)記某個(gè)節(jié)點(diǎn)是否訪問(wèn)過(guò)了index=1; %從index節(jié)點(diǎn)開(kāi)始搜索 %判斷是否對(duì)所有節(jié)點(diǎn)都找的最短路徑了。可能會(huì)有源點(diǎn)沒(méi)有路徑到目標(biāo)節(jié)點(diǎn)的情況，那就無(wú)限循環(huán)了while sum(visit)=m %沒(méi)有出隊(duì)列操作，不過(guò)通過(guò)visit來(lái)等價(jià)的表示了 visit(index)=1; %標(biāo)記第index節(jié)點(diǎn)為已入列的節(jié)點(diǎn) S pa=relax(S,pa,A,visit,index,m); %松弛，如果兩個(gè)節(jié)點(diǎn)間有更短的距離，則用更短的距離 index=extract_min(S,visit,index,m); %使用已訪問(wèn)的最小的節(jié)點(diǎn)作為下一次搜索的開(kāi)始節(jié)點(diǎn) end%最終我們需要的就是這兩個(gè)值S %源點(diǎn)到其他每一點(diǎn)的距離pa %其他每一節(jié)點(diǎn)的前趨%算法到此結(jié)束，下面只是為了形象的表示而寫(xiě)的。re=;for i=2:m re=re;pa(i) i A(pa(i),i);endA=compresstable2matrix(re); %從鄰接壓縮表構(gòu)造圖的矩陣表示figure;netplot(A,1) %形象表示relax.m%邊緣松弛，使用更短的距離作為節(jié)點(diǎn)的值function S pa=relax(S,pa,A,visit,index,m) i=index; for j=1:m if A(i,j)=inf & visit(j)=1 %搜索沒(méi)有標(biāo)記過(guò)的節(jié)點(diǎn) if S(j)S(i)+A(i,j) %將較小的值賦給正在搜尋的節(jié)點(diǎn) S(j)=S(i)+A(i,j); pa(j)=i; end end end endextract_min.m%提取隊(duì)列中尚未標(biāo)記的最小的值的序號(hào)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); %邊緣松弛，取兩節(jié)點(diǎn)間最小權(quán)值作為實(shí)際權(quán)值 pa(j)=i; %尋找前趨 end end end end endend%最終我們需要的就是這兩個(gè)值S %源點(diǎn)到其他每一點(diǎn)的距離pa %其他每一節(jié)點(diǎn)的前趨%算法到此結(jié)束，下面只是為了形象的表示而寫(xiě)的。re=;for i=2:m re=re;pa(i) i A(pa(i),i);endA=compresstable2matrix(re); %從鄰接壓縮表構(gòu)造圖的矩陣表示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:n A(b(i,1),b(i,2)=b(i,3); endend11.如此經(jīng)典的算法竟一直沒(méi)有單獨(dú)的實(shí)現(xiàn)過(guò)，真是遺憾啊。廣度優(yōu)先搜索在過(guò)去實(shí)現(xiàn)的二值圖像連通區(qū)域標(biāo)記和prim最小生成樹(shù)算法時(shí)已經(jīng)無(wú)意識(shí)的用到了，深度優(yōu)先搜索倒是沒(méi)用過(guò)。這次單獨(dú)的將兩個(gè)算法實(shí)現(xiàn)出來(lái)，因?yàn)樗惴ū旧砗蛨D像沒(méi)什么關(guān)系，所以更純粹些。廣度優(yōu)先搜索是從某一節(jié)點(diǎn)開(kāi)始，搜索與其線連接的所有節(jié)點(diǎn)，按照廣度方向像外擴(kuò)展，直到不重復(fù)遍歷所有節(jié)點(diǎn)。深度優(yōu)先搜索是從某一節(jié)點(diǎn)開(kāi)始，沿著其搜索到的第一個(gè)節(jié)點(diǎn)不斷深入下去，當(dāng)無(wú)法再深入的時(shí)候，回溯節(jié)點(diǎn)，然后再在回溯中的某一節(jié)點(diǎn)開(kāi)始沿另一個(gè)方向深度搜索，直到不重復(fù)的遍歷所有節(jié)點(diǎn)。廣度優(yōu)先搜索用的是隊(duì)列作為臨時(shí)節(jié)點(diǎn)存放處；深度優(yōu)先搜索可以遞歸實(shí)現(xiàn)（算法導(dǎo)論就是用遞歸實(shí)現(xiàn)的偽代碼），不過(guò)我這里是用棧作為臨時(shí)節(jié)點(diǎn)存放處。感覺(jué)也沒(méi)什么好介紹的了，抄算法導(dǎo)論上的介紹也沒(méi)什么意思，所有的內(nèi)容都是書(shū)上的，真正學(xué)東西還是要看書(shū)。下面是運(yùn)行結(jié)果：原連通圖：廣度優(yōu)先搜索：深度優(yōu)先搜索：matlab代碼如下，其中的畫(huà)圖函數(shù)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); %從鄰接壓縮表構(gòu)造圖的矩陣表示netplot(A,1) %形象表示head=1; %隊(duì)列頭tail=1; %隊(duì)列尾，開(kāi)始隊(duì)列為空，tail=headqueue(head)=1; %向頭中加入圖第一個(gè)節(jié)點(diǎn)head=head+1; %隊(duì)列擴(kuò)展flag=1; %標(biāo)記某個(gè)節(jié)點(diǎn)是否訪問(wèn)過(guò)了re=; %最終結(jié)果while tail=head %判斷隊(duì)列是否為空 i=queue(tail); %取隊(duì)尾節(jié)點(diǎn) for j=1:m if A(i,j)=1 & isempty(find(flag=j,1) %如果節(jié)點(diǎn)相連并且沒(méi)有訪問(wèn)過(guò) queue(head)=j; %新節(jié)點(diǎn)入列 head=head+1; %擴(kuò)展隊(duì)列 flag=flag j; %對(duì)新節(jié)點(diǎn)進(jìn)行標(biāo)記 re=re;i j; %將邊存入結(jié)果 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); %從鄰接壓縮表構(gòu)造圖的矩陣表示netplot(A,1) %形象表示top=1; %堆棧頂stack(top)=1; %將第一個(gè)節(jié)點(diǎn)入棧flag=1; %標(biāo)記某個(gè)節(jié)點(diǎn)是否訪問(wèn)過(guò)了re=; %最終結(jié)果while top=0 %判斷堆棧是否為空 pre_len=length(stack); %搜尋下一個(gè)節(jié)點(diǎn)前的堆棧長(zhǎng)度 i=stack(top); %取堆棧頂節(jié)點(diǎn) for j=1:m if A(i,j)=1 & isempty(find(flag=j,1) %如果節(jié)點(diǎn)相連并且沒(méi)有訪問(wèn)過(guò) top=top+1; %擴(kuò)展堆棧 stack(top)=j; %新節(jié)點(diǎn)入棧 flag=flag j; %對(duì)新節(jié)點(diǎn)進(jìn)行標(biāo)記 re=re;i j; %將邊存入結(jié)果 break; end end if length(stack)=pre_len %如果堆棧長(zhǎng)度沒(méi)有增加，則節(jié)點(diǎn)開(kāi)始出棧 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.模擬退火首先從某個(gè)初始候選解開(kāi)始，當(dāng)溫度大于0時(shí)執(zhí)行循環(huán)。在循環(huán)中，通過(guò)隨機(jī)擾動(dòng)產(chǎn)生一個(gè)新的解，然后求得新解和原解之間的能量差，如果差小于0，則采用新解作為當(dāng)前解。如果差大于0，則采用一個(gè)當(dāng)前溫度與能量差成比例的概率來(lái)選擇是否接受新解。溫度越低，接受的概率越小，差值越大，同樣接受概率越小。是否接受的概率用此公式計(jì)算：p=exp(-E/T)。這里E為新解與原解的差，T為當(dāng)前的溫度。由于溫度隨迭代次數(shù)逐漸降低，因此獲得一個(gè)較差的解的概率較小。典型的模擬退火算法還使用了蒙特卡洛循環(huán)，在溫度降低之前，通過(guò)多次迭代來(lái)找到當(dāng)前溫度下比較好的解。這里使用模擬退火解旅行商問(wèn)題，因?yàn)檫@個(gè)問(wèn)題本身是一個(gè)NP難問(wèn)題，所以也就求不到最優(yōu)解，不過(guò)應(yīng)該可以求得一個(gè)比較好的解，然后再手工優(yōu)化。具體到程序中，這里的隨機(jī)擾動(dòng)就是隨機(jī)置換兩個(gè)城市的位置，能量就是旅行商路線的總長(zhǎng)度。算法結(jié)果如下：初始旅行商路線：最終求得的旅行商路線：每次迭代求得的旅行距離：matlab代碼如下：main.mclear all;close all;clcn=20; %城市個(gè)數(shù)temperature=100*n; %初始溫度iter=100; %內(nèi)部蒙特卡洛循環(huán)迭代次數(shù)%隨機(jī)初始化城市坐標(biāo)city=struct();for i=1:n city(i).x=floor(1+100*rand(); city(i).y=floor(1+100*rand();endl=1; %統(tǒng)計(jì)迭代次數(shù)len(l)=computer_tour(city,n); %每次迭代后的路線長(zhǎng)度 netplot(city,n); %初始旅行路線while temperature0.001 %停止迭代溫度 for i=1:iter %多次迭代擾動(dòng)，一種蒙特卡洛方法，溫度降低之前多次實(shí)驗(yàn) len1=computer_tour(city,n); %計(jì)算原路線總距離 tmp_city=perturb_tour(city,n); %產(chǎn)生隨機(jī)擾動(dòng) len2=computer_tour(tmp_city,n); %計(jì)算新路線總距離 delta_e=len2-len1; %新老距離的差值，相當(dāng)于能量 if delta_erand() %以概率選擇是否接受新解 city=tmp_city; %可能得到較差的解 end end end l=l+1; len(l)=computer_tour(city,n); %計(jì)算新路線距離 temperature=temperature*0.99; %溫度不斷下降 end figure;netplot(city,n); %最終旅行路線figure;plot(len) computer_tour.mfunction len=computer_tour(city,n) %計(jì)算路線總長(zhǎng)度，每個(gè)城市只計(jì)算和下家城市之間的距離。 len=0; for i=1:n-1 len=len+sqrt(city(i).x-city(i+1).x)2+(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) %隨機(jī)置換兩個(gè)不同的城市的坐標(biāo) %產(chǎn)生隨機(jī)擾動(dòng) 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); city(p2)=tmp;endnetplot.mfunction netplot(city,n) %連線各城市，將路線畫(huà)出來(lái) 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); %只連線當(dāng)前城市和下家城市 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.還是這本書(shū)上的內(nèi)容，不過(guò)我看演化計(jì)算這一章是倒著看的，這里練習(xí)的算法正好和書(shū)中介紹的順序是相反的。演化策略是最古老的的演化算法之一，和上一篇DE算法類似，都是基于種群的隨機(jī)演化產(chǎn)生最優(yōu)解的算法。算法步驟如下：1.設(shè)定種群個(gè)體數(shù)和需要迭代的次數(shù)。2.選擇父代中的個(gè)體按照公式z1=sqrt(-2*ln(u1)*sin(2*pi*u2)*m，z2=sqrt(-2*ln(u1)*cos(2*pi*u2)*m進(jìn)行演化。這里u1,u2都是隨機(jī)值，m是控制因子，演化次數(shù)越多m，m越小，父代通過(guò)與z1，z2相加得到后代。3.計(jì)算后代的適應(yīng)性。4.選擇后代中最優(yōu)的適應(yīng)性作為全局最優(yōu)適應(yīng)性。其實(shí)整個(gè)過(guò)程和DE非常類似。過(guò)程都是隨機(jī)變異，求適應(yīng)性，再找最優(yōu)。我還試著將z1和z2橫設(shè)為1，竟也能得到非常好的解。算法結(jié)果如下：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); %目標(biāo)函數(shù)，高斯函數(shù)mesh(img);hold on;n=50; %種群個(gè)體的數(shù)量iter=100; %迭代次數(shù)%初始化種群，定義結(jié)構(gòu)體par=struct();for i=1:n par(i).x=-100+200*rand(); %個(gè)體的x特征在-100 100隨機(jī)初始化 par(i).y=-100+200*rand(); %個(gè)體的y特征在-100 100隨機(jī)初始化 par(i).fit=compute_fit(par(i); %個(gè)體在x,y處的適應(yīng)度endpar_best=par(1); %初始化種群中最佳個(gè)體for k=1:iter %迭代次數(shù) plot3(par_best.x+100,par_best.y+100,par_best.fit,g*); %畫(huà)出最佳個(gè)體的位置，+100為相對(duì)偏移 par par_best=select_and_recombin(par,par_best,n,k,iter); %差異演化函數(shù)endselect_and_recombin.mfunction next_par par_best=select_and_recombin(par,par_best,n,k,iter) mul=(iter-k)/iter; %限制進(jìn)化因子，代數(shù)越高變異越小 next_par=par; %新種群 for i=1:n %產(chǎn)生變異隨機(jī)數(shù) 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; %計(jì)算變異后個(gè)體的適應(yīng)度 next_par(i).fit=compute_fit(next_par(i); %如果新個(gè)體沒(méi)有變異前個(gè)體適應(yīng)度高，新個(gè)體還原為舊個(gè)體 if par(i).fitnext_par(i).fit next_par(i)=par(i); end %如果變異后適應(yīng)度高于種群最高適應(yīng)個(gè)體，則更新種群適應(yīng)度最高個(gè)體 if next_par(i).fitpar_best.fit par_best=next_par(i); end end endcompute_fit.mfunction re=compute_fit(par)