MATLAB程序大全

精品文檔 1歡迎下載 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é)果 原圖 精品文檔 2歡迎下載 結(jié)果 matlab 代碼如下 clear all close all clc img imread pan jpg 精品文檔 3歡迎下載 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 x 2 dis y 2 if l r1 theta 0 if y re m 2 theta atan2 dis y dis x end if y 1 角點(diǎn)區(qū)域 H 0 harris 角點(diǎn)檢測(cè)就用到了第三類判斷 當(dāng)然 在實(shí)際應(yīng)用的時(shí)候 H 和 K 的值肯定都不會(huì)是理想 所以我用的都是近似判斷 處理結(jié)果如下 原圖 平坦區(qū)域 精品文檔 9歡迎下載 邊緣區(qū)域 角點(diǎn)區(qū)域 好像也不全角點(diǎn) 求角點(diǎn)還是 harris 好了 結(jié)構(gòu)張量行列式與跡的關(guān)系 精品文檔 10歡迎下載 其中紅框?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 end lambda k K H 精品文檔 11歡迎下載 k k 1 end end figure 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)后 精品文檔 12歡迎下載 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 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 精品文檔 17歡迎下載 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è)上的 精品文檔 18歡迎下載 matlab 代碼如下 netplot 和 compresstable2matrix 和上一篇使用的一樣 main m clear all close all clc 初始化鄰接壓縮表 1 2 10 表示從節(jié)點(diǎn) 1 到節(jié)點(diǎn) 2 邊的權(quán)重為 10 b 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)到自己的距離為 0 pa zeros 1 m 存儲(chǔ)每個(gè)節(jié)點(diǎn)的前驅(qū) 在松弛過(guò)程中賦值 精品文檔 19歡迎下載 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 end A 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 pa j i end end 精品文檔 20歡迎下載 end end extract 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 end 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 end A compresstable2matrix re 從鄰接壓縮表構(gòu)造圖的矩陣表示 figure netplot A 1 形象表示 compresstable2matrix m function 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 end end 精品文檔 23歡迎下載 11 如此經(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)先搜索 精品文檔 24歡迎下載 深度優(yōu)先搜索 matlab 代碼如下 其中的畫(huà)圖函數(shù) netplot m BFS m clear 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 head queue 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 新節(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 精品文檔 25歡迎下載 tail tail 1 end A compresstable2matrix re figure netplot A 1 DFS m clear 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 擴(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 end A compresstable2matrix re 精品文檔 26歡迎下載 figure netplot A 1 compresstable2matrix m function 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 end end 12 模擬退火首先從某個(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é)果如下 初始旅行商路線 精品文檔 27歡迎下載 最終求得的旅行商路線 每次迭代求得的旅行距離 matlab 代碼如下 main m 精品文檔 28歡迎下載 clear all close all clc n 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 end l 1 統(tǒng)計(jì)迭代次數(shù) len l computer tour city n 每次迭代后的路線長(zhǎng)度 netplot city n 初始旅行路線 while temperature 0 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 最終旅行路線 精品文檔 29歡迎下載 figure plot len computer tour m function 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 end perturb tour m function 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 end netplot m function netplot city n 連線各城市 將路線畫(huà)出來(lái) hold on for i 1 n 1 plot city i x city i y r 精品文檔 30歡迎下載 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 end 13 還是這本書(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 m 精品文檔 31歡迎下載 clear all close all clc x y meshgrid 100 100 100 100 sigma 50 img 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 目標(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)度 end par 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ù) end select and recombin m function 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 精品文檔 32歡迎下載 變異 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 fit next par i fit next par i par i end 如果變異后適應(yīng)度高于種群最高適應(yīng)個(gè)體 則更新種群適應(yīng)度最高個(gè) 體 if next par i fit par best fit par best next par i end end end compute fit m function re compute fit par x par x y par y sigma 50 if x100 y100 re 0 超出范圍適應(yīng)度為 0 else 否則適應(yīng)度按目標(biāo)函數(shù)求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 14 這兩天在看 M Tim Jones 的 人工智能 書(shū)中不只介紹原理 而且都有相應(yīng)的 c 代碼實(shí)現(xiàn) 雖然代碼不完全 不過(guò)缺少的部分完全可以自己補(bǔ)完 差異演化和昨天實(shí)現(xiàn)的 PSO 很類似 都屬于優(yōu)化算法 算法步驟 1 設(shè)定種群個(gè)體個(gè)數(shù)和需要迭代的次數(shù) 當(dāng)然也可以設(shè)定條件 然后判斷是否停止迭代 2 定義交叉概率 CR 個(gè)體有一定概率進(jìn)行變異 如果變異則進(jìn)行第 3 步 如果不變異則下一代 個(gè)體和當(dāng)前個(gè)體一樣 精品文檔 33歡迎下載 3 在種群中隨機(jī)選出三個(gè)互不相同的個(gè)體進(jìn)行變異 變異公式如下 有博士論文總結(jié)了一大堆 變異公式 Xi g 1 Xr1 g F Xr2 g Xr3 g 其中 X 是種群中個(gè)體特征向量 這里就是 x 和 y 坐標(biāo) g 是種群的代數(shù) i 代表當(dāng)前個(gè)體 r1 r2 r3 是和 i 不同 并且也互不相同的個(gè)體 4 計(jì)算變異后個(gè)體的適應(yīng)度 如果變異后適應(yīng)度不如變異前 那么將變異后個(gè)體重新恢復(fù)為變 異前個(gè)體 5 比較變異后個(gè)體適應(yīng)度和種群最優(yōu)個(gè)體適應(yīng)度 將適應(yīng)度高的個(gè)體賦給種群最優(yōu)個(gè)體 最后當(dāng)然還是需要目標(biāo)函數(shù)才能計(jì)算適應(yīng)度 算法結(jié)果如下圖 小綠點(diǎn)代表種群最優(yōu)個(gè)體適應(yīng)度 matlab 代碼如下 main m clear all close all clc x y meshgrid 100 100 100 100 sigma 50 img 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 目標(biāo)函數(shù) 高 斯函數(shù) mesh img hold on n 20 種群個(gè)體的數(shù)量 初始化種群 定義結(jié)構(gòu)體 par struct for i 1 n 精品文檔 34歡迎下載 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 0 個(gè)體適應(yīng)度為 0 初始化 end par best par 1 初始化種群中最佳個(gè)體 for k 1 100 迭代次數(shù) plot3 par best x 100 par best y 100 par best fit g 畫(huà)出最佳 個(gè)體的位置 100 為相對(duì)偏移 par par best select and recombine par par best n 差異演 化函數(shù) end select and recombine m function next par par best select and recombine par par best n F 0 5 加速因子 CR 0 8 變異率 next par par 新種群 for i 1 n while 1 在原種群中任選三個(gè)互不相同的個(gè)體進(jìn)行交叉變異 r1 floor 1 20 rand r2 floor 1 20 rand r3 floor 1 20 rand if i r1 end end if rand next par i fit next par i par i 精品文檔 35歡迎下載 end 如果變異后適應(yīng)度高于種群最高適應(yīng)個(gè)體 則更新種群適應(yīng)度最高個(gè) 體 if next par i fit par best fit par best next par i end end end compute fit m function re compute fit par x par x y par y sigma 50 if x100 y100 re 0 超出范圍適應(yīng)度為 0 else 否則適應(yīng)度按目標(biāo)函數(shù)求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 15 算法沒(méi)有和圖像處理直接相關(guān) 不過(guò)對(duì)于圖像分類中的模式識(shí)別相關(guān)算法 也許會(huì)用到這個(gè)優(yōu) 化算法 算法步驟 1 首先確定粒子個(gè)數(shù)與迭代次數(shù) 2 對(duì)每個(gè)粒子隨機(jī)初始化位置與速度 3 采用如下公式更新每個(gè)粒子的位置與速度 Px Px Pv t 位置更新公式 Pv Pv c1 rand Gx Px c2 rand PBx Px 速度更新公式 這里 c1 和 c2 是加速因子 和梯度下降算法那里的加速因子我感覺(jué)很類似 Gx 是粒子群中最佳粒子的位置 PBx 為當(dāng)前粒子最佳位置 4 每次迭代 首先檢查新粒子適應(yīng)度是否高于原最優(yōu)適應(yīng)度 如果高于則對(duì)自己的位置和適應(yīng) 度進(jìn)行更新 然后再判斷此粒子適應(yīng)度是否高于全局最優(yōu)粒子 如果高于則更新全局最優(yōu)粒子 適應(yīng)度和位置 因?yàn)樽约翰皇侵饕芯窟@方面算法的 所以還有一些疑問(wèn) 自問(wèn)自答 精品文檔 36歡迎下載 1 算法需要目標(biāo)函數(shù) 如果沒(méi)有目標(biāo)函數(shù)怎么辦 也許就不用這個(gè)算法了 或者其他什么算法 先求出了目標(biāo)函數(shù)了 2 既然給了目標(biāo)函數(shù) 那么直接遍歷所有值再 max 應(yīng)該就能求得最佳位置 而 PSO 算法是不 是只是為了減少運(yùn)算量 比如我這里 200 200 的矩陣 本來(lái)需要計(jì)算 40000 次函數(shù) 而 PSO 只 計(jì)算了 100 次函數(shù)就得到近似最優(yōu)解了 難怪叫優(yōu)化算法 反正我暫時(shí)只能這樣理解了 其他細(xì)節(jié)代碼注釋的很清楚了 下圖展示了一個(gè) PSO 的運(yùn)行結(jié)果 目標(biāo)函數(shù)是高斯函數(shù) 綠點(diǎn)代表最佳粒子的位置 matlab 代碼如下 main m clear all close all clc x y meshgrid 100 100 100 100 sigma 50 img 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 目標(biāo)函數(shù) 高 斯函數(shù) mesh img hold on n 10 粒子群粒子個(gè)數(shù) 初始化粒子群 定義結(jié)構(gòu)體 結(jié)構(gòu)體中八個(gè)元素 分別是粒子坐標(biāo) 粒子速度 粒子適應(yīng)度 粒子最佳適應(yīng) 度 粒子最佳坐標(biāo) par struct for i 1 n par i x 100 200 rand 100 100 對(duì) x 位置隨機(jī)初始化 par i y 100 200 rand 100 100 對(duì) y 位置隨機(jī)初始化 par i vx 1 2 rand 1 1 對(duì) vx 速度隨機(jī)初始化 精品文檔 37歡迎下載 par i vy 1 2 rand 1 1 對(duì) vy 速度隨機(jī)初始化 par i fit 0 粒子適應(yīng)度為 0 初始化 par i bestfit 0 粒子最佳適應(yīng)度為 0 初始化 par i bestx par i x 粒子 x 最佳位置初始化 par i besty par i y 粒子 y 最佳位置初始化 end par best par 1 初始化粒子群中最佳粒子 for k 1 10 plot3 par best x 100 par best y 100 par best fit g 畫(huà)出最佳 粒子的位置 100 為相對(duì)偏移 for p 1 n par p par best update par par p par best 更新每個(gè)粒子 信息 end end update par m function par par best update par par par best Px Px Pv t 這里 t 1 Px 為當(dāng)前粒子的位置 Pv 為當(dāng)前粒子的速度 par x par x par vx par y par x par vy par fit compute fit par 計(jì)算當(dāng)前粒子適應(yīng)度 Pv Pv c1 rand Gx Px c2 rand PBx Px 這里 c1 c2 為加速因子 Gx 為具有最佳適應(yīng)度粒子的位置 PBx 為當(dāng)前粒子的最佳位置 c1 1 c2 1 par vx par vx c1 rand par best x par x c2 rand par bestx par x par vy par vy c1 rand par best y par y c2 rand par besty par y if par fit par bestfit 如果當(dāng)前粒子適應(yīng)度要好于當(dāng)前粒子最佳 適應(yīng)度 par bestfit par fit 則更新當(dāng)前粒子最佳適應(yīng)度 par bestx par x 更新當(dāng)前粒子最佳位置 par besty par y 精品文檔 38歡迎下載 if par bestfit par best fit 如果當(dāng)前粒子最佳適應(yīng)度好于最 佳粒子適應(yīng)度 par best fit par bestfit 則更新最佳粒子適應(yīng)度 par best x par x 更新最佳粒子位置 par best y par y end end end compute fit m function re compute fit par x par x y par y sigma 50 if x100 y100 re 0 超出范圍適應(yīng)度為 0 else 否則適應(yīng)度按目標(biāo)函數(shù)求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 16 Log Polar 就是所謂的極坐標(biāo)變換 通過(guò)變換能將笛卡爾坐標(biāo)系中的圖像映射到極坐標(biāo)系中 確切的來(lái)說(shuō)我這里不算是 Log Polar 因?yàn)?Log Polar 是將圖像從 x y 映射到 log r theta 而我是將圖像從 x y 映射到 r theta 原理是一樣的 變換公式是 x r cos theta y r sin theta 其中 r 代表極徑 theta 代表極角 這些應(yīng)該學(xué)過(guò)高中的都會(huì)的 下圖能形象的表示變換 精品文檔 39歡迎下載 似乎有通過(guò)此變換求圖像特征的算法 不過(guò)具體操作我也不清楚 我這里只是實(shí)現(xiàn)了變換的操作 原圖如下 以圖像中心作為極坐標(biāo)原點(diǎn)進(jìn)行極坐標(biāo)變換 matlab 代碼如下 精品文檔 40歡迎下載 clear all close all clc img double imread lena jpg imshow img m n size img ox oy ginput 獲得極坐標(biāo)變換的原點(diǎn) oy round oy ox round ox 求中心點(diǎn)到圖像四個(gè)角的距離 up left sqrt oy 0 2 ox 0 2 up right sqrt oy 0 2 ox n 2 down left sqrt oy m 2 ox 0 2 down right sqrt oy m 2 ox n 2 求中心點(diǎn)距離四角距離的最大值 作為變換后圖像的高 這個(gè)最大值也是極坐標(biāo)變換的極徑 radius round max up left up right down left down right angle 360 變換后圖像的寬 imgn zeros radius angle for i 1 radius 縱坐標(biāo)代表極徑 不同情況不一樣 for j 1 angle 橫坐標(biāo)代表極角 恒為 360 oy ox 作為極坐標(biāo)變換中心坐標(biāo) 需要作為偏移量相加 h oy round i sin j pi 180 w ox round i cos j pi 180 if h 0 最鄰近插值 因?yàn)榍懊媲?h w 用了 round end end end figure imshow imgn 17 通過(guò)本篇和上一篇的結(jié)合 應(yīng)該就能做出拉普拉斯圖像融合了 這里用的方法很簡(jiǎn)單 就是用模板和兩個(gè)圖像相乘 然后對(duì)處理后的兩個(gè)圖像再相加就可以了 精品文檔 41歡迎下載 拉普拉斯融合就是對(duì)金字塔的每一層圖像做這樣的操作 然后再重構(gòu)就行了 先看看這里的效果 原圖 apple 原圖 orange 模板圖像 實(shí)際處理時(shí)需要對(duì)模板進(jìn)行模糊 精品文檔 42歡迎下載 融合后圖像 matlab 代碼如下 clear all close all clc apple double imread apple jpg orange double imread orange jpg m n size apple maska ones m 220 zeros m n 220 masko 1 maska 精品文檔 43歡迎下載 maska imfilter maska fspecial gaussian 20 10 replicate same masko imfilter masko fspecial gaussian 20 10 replicate same re maska apple masko orange imshow re 18 過(guò)去研究過(guò)高斯金字塔 不過(guò)今天看來(lái) 當(dāng)時(shí)似乎搞錯(cuò)了 金字塔上下層應(yīng)該是采樣的關(guān)系 而不是縮放的關(guān)系 而且不同層的 sigma 應(yīng)該也有所不同 拉普拉斯金字塔在圖像融合中有所應(yīng)用 方法是首先對(duì)兩個(gè)待融合圖像求拉普拉斯殘差金字塔 然后用模板對(duì)每一級(jí)殘差圖像進(jìn)行融合得到融合后圖像的殘差金字塔 然后對(duì)這個(gè)金字塔進(jìn)行 重構(gòu)就能得到最終的融合圖像 圖像各尺度細(xì)節(jié)得到保留 注 融合時(shí)模板一般會(huì)先用高斯 函數(shù)模糊一下 不過(guò)這里不實(shí)現(xiàn)融合 只實(shí)現(xiàn)拉普拉斯金字塔的建立 建立金字塔的方框圖如下 摘自 數(shù)字圖像處理 第二版 通過(guò) j 1 級(jí)近似圖像和 j 級(jí)殘差圖像就能恢復(fù) j 級(jí)輸入圖像 近似濾波器就是拉普拉斯濾波模板 插入濾波器是對(duì)拉普拉斯濾波模板采樣形成的新濾波模板 使用的拉普拉斯濾波模板如下 w 1 256 1 4 6 4 1 4 16 24 16 4 6 24 36 24 6 4 16 24 16 4 1 4 6 4 1 產(chǎn)生的殘差金字塔 精品文檔 44歡迎下載 前面四張圖像都是殘差 最后一張圖像是原圖的 4 級(jí)近似 通過(guò)最后一張圖和殘差圖 就能重 構(gòu)回原圖 matlab 代碼如下 main m clear all close all clc img double imread lena jpg m n size img w 1 256 1 4 6 4 1 拉普拉斯濾波器 4 16 24 16 4 6 24 36 24 6 4 16 24 16 4 1 4 6 4 1 imgn 1 img for i 2 5 濾波 下采樣 imgn i imfilter imgn i 1 w replicate imgn i imgn i 1 2 size imgn i 1 1 1 2 size imgn i 2 1 i 1 級(jí)近似 end for i 5 1 2 調(diào)整圖像大小 imgn i 1 imgn i 1 1 2 size imgn i 1 1 2 size imgn i 2 end for i 1 4 獲得殘差圖像 i 級(jí)預(yù)測(cè)殘差 精品文檔 45歡迎下載 imgn i imgn i expand imgn i 1 w end for i 4 1 1 殘差圖像重構(gòu)原圖像 imgn i imgn i expand imgn i 1 w end imshow uint8 imgn 1 expand m function re expand img w img double img w w 4 m n size img M N size w 插入濾波器 w up left w 1 2 M 1 2 N w up right w 1 2 M 2 2 N w down left w 2 2 M 1 2 N w down right w 2 2 M 2 2 N img up left imfilter img w up left replicate same img up right imfilter img w up right replicate same img down left imfilter img w down left replicate same img down right imfilter img w down right replicate same re zeros m 2 n 2 上采樣 re 1 2 m 2 1 2 n 2 img up left re 2 2 m 2 1 2 n 2 img up right re 1 2 m 2 2 2 n 2 img down left re 2 2 m 2 2 2 n 2 img down right end 19 RANSAC 在圖像拼接中有所使用 有時(shí)候也在圖像理解的相關(guān)算法中有所使用 算法簡(jiǎn)介如下 摘自 圖像處理 分析與機(jī)器視覺(jué) 第 3 版 精品文檔 46歡迎下載 1 假設(shè)我們要將 n 個(gè)數(shù)據(jù)點(diǎn) X x1 x1 xn 擬合為一個(gè)由至少 m 個(gè)點(diǎn)決定的模型 m n 對(duì) 于直線 m 2 我這里實(shí)際是兩個(gè)不同均值 協(xié)方差高斯分布產(chǎn)生的數(shù)據(jù) 2 設(shè)迭代計(jì)數(shù) k 1 3 從 X 中隨機(jī)選取 m 個(gè)項(xiàng)并擬合一個(gè)模型 我這里直線擬合 選了 2 個(gè)項(xiàng) 4 給定偏差 計(jì)算 X 中相對(duì)于模型的殘差在偏差 的個(gè)數(shù) 如果元素個(gè)數(shù)大于一個(gè)閾值 t 根 據(jù)一致點(diǎn)集重新擬合模型 可以利用最小二乘或其變種 算法終止 我這里的偏差為 1 閾值為數(shù)據(jù)個(gè)數(shù)的 2 3 5 設(shè) k k 1 如果 k 小于一個(gè)事先給定的 K 跳至第 3 步 否則采用具有迄今最大的一致點(diǎn)集模 型 或算法失敗 運(yùn)行效果如下 紅圈是所有的數(shù)據(jù) 藍(lán)叉是符合擬合模型的數(shù)據(jù) matlab 代碼如下 main m close all clear all clc 真實(shí)的數(shù)據(jù) mu 0 0 均值 S 1 2 5 2 5 8 協(xié)方差 data1 mvnrnd mu S 200 產(chǎn)生 200 個(gè)高斯分布數(shù)據(jù) 噪聲數(shù)據(jù) mu 2 2 S 8 0 0 8 data2 mvnrnd mu S 100 產(chǎn)生 100 個(gè)噪聲數(shù)據(jù) 精品文檔 47歡迎下載 data data1 data2 plot data 1 data 2 ro 顯示全部數(shù)據(jù) K 100 設(shè)置最大迭代次數(shù) sigma 1 設(shè)置擬合直線與數(shù)據(jù)距離的偏差 pretotal 0 符合擬合模型的數(shù)據(jù)的個(gè)數(shù) k 1 while pretotal size data 1 2 3 產(chǎn)生兩個(gè)隨機(jī)索引 找樣本用 floor 向下取整 samp1