MATLAB程序大全

精品文档 1欢迎下载 1 全景图到穹景图 这个程序我最初是用 FreeImage 写的 这两天改成了 matlab 再不贴上来 我就要忘了 看到一篇文章有这样的变换 挺有意思的 就拿来试了一下 文章点此 全景图到穹顶图变换 通俗的说就是将全景图首尾相接做成一个圆环的样子 先看下面这张图 下面的矩形就是我们要处理的全景图 上面的矩形是变换后的图像 下面图像的底边对应穹顶 图的内圆 顶边对应穹顶图的外圆 当然 反过来也是可以的 程序流程 1 定义穹顶图内圆和外圆的半径 变换后的像素就填充在这个内外半径的圆环中 2 遍历穹顶图 当所处理当前像素位于圆环内 则通过极坐标反变换去全景图中寻找相应位置 的像素进行填充 3 遍历完图像就行了 用的技巧和图像旋转或放大缩小都是类似的 处理结果 原图 精品文档 2欢迎下载 结果 matlab 代码如下 clear all close all clc img imread pan jpg 精品文档 3欢迎下载 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 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 角点区域 H 0 harris 角点检测就用到了第三类判断 当然 在实际应用的时候 H 和 K 的值肯定都不会是理想 所以我用的都是近似判断 处理结果如下 原图 平坦区域 精品文档 9欢迎下载 边缘区域 角点区域 好像也不全角点 求角点还是 harris 好了 结构张量行列式与迹的关系 精品文档 10欢迎下载 其中红框为平坦区域 黄框为边缘区域 铝框为角点区域 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 结构张量 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 模糊集图像增强 算法有很多变种 不过主要就是以下三步 1 设计隶属度函数将图像从空间域变换到模糊集域 2 设计模糊增强算子 在模糊集域对图像进行处理 3 根据第 1 步的隶属度函数重新将图像从模糊集域变换到空间域 这和频域处理中的变换反变换不是很像么 我使用的隶属度函数和模糊增强算子在这篇论文里 也算相关算法的经典论文了 处理结果如下 原图 模糊集增强后 精品文档 12欢迎下载 matlab 代码如下 clear all close all clc img double imread lena 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 head head 1 flag flag j pa j i end end tail tail 1 end if pa m 0 如果搜索不到汇节点 退出循环 break end path i m 从汇节点开始 k 0 路径包含的边的个数 while i 1 使用前趋构造从源节点到汇节点的路径 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 最大流 原网络包含这个网络 我只能这样表示了 end 使用新的图 A 进入下一循环 从新开始找增广 路径 end figure netplot maxflow 1 9 单元最短路径 图的相关算法也算是自己的一个软肋了 当年没选修图论也是一大遗憾 图像处理中 也有使用图论算法作为基础的相关算法 比如图割 这个算法就需要求最大流 最小割 所以熟悉一下图论算法对于图像处理还是很有帮助的 Dijkstra 和 Bellman Ford 类似 都是解决单源最短路径问题 不同的是这个方法只能解决边 为非负的问题 实现的好的 Dijkstra 算法运行时间要快于 Bellman ford 算法步骤如下 1 首先设置队列 所有节点入列 源节点值为 0 其他节点值为无穷 2 然后在队列中找值最小的节点并出列 3 计算出列的节点所有后继节点的距离 4 松弛方法 如果新计算的距离小于上次计算的距离 那么更新距离 即将后继节点值设为较 小的距离 并将后继节点的前趋设为当前的出列节点 5 对剩余的节点队列继续找最小值并出列 不断循环 2 3 4 步直到队列中没有节点了 步骤是上面没错 不过我程序中没有完全按照上述的步骤实现 不同的地方在于我没有做出列 操作 而是通过标记节点的形式实现的 运行结果如下 图 是图不是图片 是算法导论 367 页上的 精品文档 18欢迎下载 matlab 代码如下 netplot 和 compresstable2matrix 和上一篇使用的一样 main m clear all close all clc 初始化邻接压缩表 1 2 10 表示从节点 1 到节点 2 边的权重为 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 从邻接压缩表构造图的矩阵表示 netplot A 1 形象表示 S inf 1 m 从开始的源点到每一个节点的距离 S 1 0 源点到自己的距离为 0 pa zeros 1 m 存储每个节点的前驱 在松弛过程中赋值 精品文档 19欢迎下载 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 index m 使用已访问的最小的节点作为下 一次搜索的开始节点 end 最终我们需要的就是这两个值 S 源点到其他每一点的距离 pa 其他每一节点的前趋 算法到此结束 下面只是为了形象的表示而写的 re for i 2 m re re pa i i A pa i i end A compresstable2matrix re 从邻接压缩表构造图的矩阵表示 figure netplot A 1 形象表示 relax m 边缘松弛 使用更短的距离作为节点的值 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 提取队列中尚未标记的最小的值的序号 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 end end end end end 最终我们需要的就是这两个值 S 源点到其他每一点的距离 pa 其他每一节点的前趋 算法到此结束 下面只是为了形象的表示而写的 re for i 2 m re re pa i i A pa i i end A compresstable2matrix re 从邻接压缩表构造图的矩阵表示 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 如此经典的算法竟一直没有单独的实现过 真是遗憾啊 广度优先搜索在过去实现的二值图像连通区域标记和 prim 最小生成树算法时已经无意识的用到 了 深度优先搜索倒是没用过 这次单独的将两个算法实现出来 因为算法本身和图像没什么关系 所以更纯粹些 广度优先搜索是从某一节点开始 搜索与其线连接的所有节点 按照广度方向像外扩展 直到 不重复遍历所有节点 深度优先搜索是从某一节点开始 沿着其搜索到的第一个节点不断深入下去 当无法再深入的 时候 回溯节点 然后再在回溯中的某一节点开始沿另一个方向深度搜索 直到不重复的遍历 所有节点 广度优先搜索用的是队列作为临时节点存放处 深度优先搜索可以递归实现 算法导论就是用 递归实现的伪代码 不过我这里是用栈作为临时节点存放处 感觉也没什么好介绍的了 抄算法导论上的介绍也没什么意思 所有的内容都是书上的 真正 学东西还是要看书 下面是运行结果 原连通图 广度优先搜索 精品文档 24欢迎下载 深度优先搜索 matlab 代码如下 其中的画图函数 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 从邻接压缩表构造图的矩阵表示 netplot A 1 形象表示 head 1 队列头 tail 1 队列尾 开始队列为空 tail head queue head 1 向头中加入图第一个节点 head head 1 队列扩展 flag 1 标记某个节点是否访问过了 re 最终结果 while tail head 判断队列是否为空 i queue tail 取队尾节点 for j 1 m if A i j 1 新节点入列 head head 1 扩展队列 flag flag j 对新节点进行标记 re re i j 将边存入结果 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 从邻接压缩表构造图的矩阵表示 netplot A 1 形象表示 top 1 堆栈顶 stack top 1 将第一个节点入栈 flag 1 标记某个节点是否访问过了 re 最终结果 while top 0 判断堆栈是否为空 pre len length stack 搜寻下一个节点前的堆栈长度 i stack top 取堆栈顶节点 for j 1 m if A i j 1 扩展堆栈 stack top j 新节点入栈 flag flag j 对新节点进行标记 re re i j 将边存入结果 break end end if length stack pre len 如果堆栈长度没有增加 则节点开始出栈 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 模拟退火首先从某个初始候选解开始 当温度大于 0 时执行循环 在循环中 通过随机扰动产生一个新的解 然后求得新解和原解之间的能量差 如果差小于 0 则采用新解作为当前解 如果差大于 0 则采用一个当前温度与能量差成比例的概率来选择是否接受新解 温度越低 接受的概率越小 差值越大 同样接受概率越小 是否接受的概率用此公式计算 p exp E T 这里 E 为新解与原解的差 T 为当前的温度 由于温度随迭代次数逐渐降低 因此获得一个较差的解的概率较小 典型的模拟退火算法还使用了蒙特卡洛循环 在温度降低之前 通过多次迭代来找到当前温度 下比较好的解 这里使用模拟退火解旅行商问题 因为这个问题本身是一个 NP 难问题 所以也就求不到最优解 不过应该可以求得一个比较好的解 然后再手工优化 具体到程序中 这里的随机扰动就是随机置换两个城市的位置 能量就是旅行商路线的总长度 算法结果如下 初始旅行商路线 精品文档 27欢迎下载 最终求得的旅行商路线 每次迭代求得的旅行距离 matlab 代码如下 main m 精品文档 28欢迎下载 clear all close all clc n 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 end l 1 统计迭代次数 len l computer tour city n 每次迭代后的路线长度 netplot city n 初始旅行路线 while temperature 0 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 city 可能得到较差的解 end end end l l 1 len l computer tour city n 计算新路线距离 temperature temperature 0 99 温度不断下降 end figure netplot city n 最终旅行路线 精品文档 29欢迎下载 figure plot len computer tour m function len computer tour city n 计算路线总长度 每个城市只计算 和下家城市之间的距离 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 随机置换两个不同的城市的坐标 产生随机扰动 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 连线各城市 将路线画出来 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 只连 线当前城市和下家城市 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 还是这本书上的内容 不过我看演化计算这一章是倒着看的 这里练习的算法正好和书中介绍 的顺序是相反的 演化策略是最古老的的演化算法之一 和上一篇 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 非常类似 过程都是随机变异 求适应性 再找最优 我还试着将 z1 和 z2 横设为 1 竟也能得到非常好的解 算法结果如下 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 目标函数 高 斯函数 mesh img hold on n 50 种群个体的数量 iter 100 迭代次数 初始化种群 定义结构体 par struct for i 1 n par i x 100 200 rand 个体的 x 特征在 100 100 随机初始 化 par i y 100 200 rand 个体的 y 特征在 100 100 随机初始 化 par i fit compute fit par i 个体在 x y 处的适应度 end par 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 差异演化函数 end select and recombin m function 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 精品文档 32欢迎下载 变异 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 fit next par i fit next par i par i end 如果变异后适应度高于种群最高适应个体 则更新种群适应度最高个 体 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 超出范围适应度为 0 else 否则适应度按目标函数求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 14 这两天在看 M Tim Jones 的 人工智能 书中不只介绍原理 而且都有相应的 c 代码实现 虽然代码不完全 不过缺少的部分完全可以自己补完 差异演化和昨天实现的 PSO 很类似 都属于优化算法 算法步骤 1 设定种群个体个数和需要迭代的次数 当然也可以设定条件 然后判断是否停止迭代 2 定义交叉概率 CR 个体有一定概率进行变异 如果变异则进行第 3 步 如果不变异则下一代 个体和当前个体一样 精品文档 33欢迎下载 3 在种群中随机选出三个互不相同的个体进行变异 变异公式如下 有博士论文总结了一大堆 变异公式 Xi g 1 Xr1 g F Xr2 g Xr3 g 其中 X 是种群中个体特征向量 这里就是 x 和 y 坐标 g 是种群的代数 i 代表当前个体 r1 r2 r3 是和 i 不同 并且也互不相同的个体 4 计算变异后个体的适应度 如果变异后适应度不如变异前 那么将变异后个体重新恢复为变 异前个体 5 比较变异后个体适应度和种群最优个体适应度 将适应度高的个体赋给种群最优个体 最后当然还是需要目标函数才能计算适应度 算法结果如下图 小绿点代表种群最优个体适应度 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 目标函数 高 斯函数 mesh img hold on n 20 种群个体的数量 初始化种群 定义结构体 par struct for i 1 n 精品文档 34欢迎下载 par i x 100 200 rand 个体的 x 特征在 100 100 随机初始化 par i y 100 200 rand 个体的 y 特征在 100 100 随机初始化 par i fit 0 个体适应度为 0 初始化 end par best par 1 初始化种群中最佳个体 for k 1 100 迭代次数 plot3 par best x 100 par best y 100 par best fit g 画出最佳 个体的位置 100 为相对偏移 par par best select and recombine par par best n 差异演 化函数 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 在原种群中任选三个互不相同的个体进行交叉变异 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 如果变异后适应度高于种群最高适应个体 则更新种群适应度最高个 体 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 超出范围适应度为 0 else 否则适应度按目标函数求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 15 算法没有和图像处理直接相关 不过对于图像分类中的模式识别相关算法 也许会用到这个优 化算法 算法步骤 1 首先确定粒子个数与迭代次数 2 对每个粒子随机初始化位置与速度 3 采用如下公式更新每个粒子的位置与速度 Px Px Pv t 位置更新公式 Pv Pv c1 rand Gx Px c2 rand PBx Px 速度更新公式 这里 c1 和 c2 是加速因子 和梯度下降算法那里的加速因子我感觉很类似 Gx 是粒子群中最佳粒子的位置 PBx 为当前粒子最佳位置 4 每次迭代 首先检查新粒子适应度是否高于原最优适应度 如果高于则对自己的位置和适应 度进行更新 然后再判断此粒子适应度是否高于全局最优粒子 如果高于则更新全局最优粒子 适应度和位置 因为自己不是主要研究这方面算法的 所以还有一些疑问 自问自答 精品文档 36欢迎下载 1 算法需要目标函数 如果没有目标函数怎么办 也许就不用这个算法了 或者其他什么算法 先求出了目标函数了 2 既然给了目标函数 那么直接遍历所有值再 max 应该就能求得最佳位置 而 PSO 算法是不 是只是为了减少运算量 比如我这里 200 200 的矩阵 本来需要计算 40000 次函数 而 PSO 只 计算了 100 次函数就得到近似最优解了 难怪叫优化算法 反正我暂时只能这样理解了 其他细节代码注释的很清楚了 下图展示了一个 PSO 的运行结果 目标函数是高斯函数 绿点代表最佳粒子的位置 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 目标函数 高 斯函数 mesh img hold on n 10 粒子群粒子个数 初始化粒子群 定义结构体 结构体中八个元素 分别是粒子坐标 粒子速度 粒子适应度 粒子最佳适应 度 粒子最佳坐标 par struct for i 1 n par i x 100 200 rand 100 100 对 x 位置随机初始化 par i y 100 200 rand 100 100 对 y 位置随机初始化 par i vx 1 2 rand 1 1 对 vx 速度随机初始化 精品文档 37欢迎下载 par i vy 1 2 rand 1 1 对 vy 速度随机初始化 par i fit 0 粒子适应度为 0 初始化 par i bestfit 0 粒子最佳适应度为 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 画出最佳 粒子的位置 100 为相对偏移 for p 1 n par p par best update par par p par best 更新每个粒子 信息 end end update par m function par par best update par par par best Px Px Pv t 这里 t 1 Px 为当前粒子的位置 Pv 为当前粒子的速度 par x par x par vx par y par x par vy par fit compute fit par 计算当前粒子适应度 Pv Pv c1 rand Gx Px c2 rand PBx Px 这里 c1 c2 为加速因子 Gx 为具有最佳适应度粒子的位置 PBx 为当前粒子的最佳位置 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 如果当前粒子适应度要好于当前粒子最佳 适应度 par bestfit par fit 则更新当前粒子最佳适应度 par bestx par x 更新当前粒子最佳位置 par besty par y 精品文档 38欢迎下载 if par bestfit par best fit 如果当前粒子最佳适应度好于最 佳粒子适应度 par best fit par bestfit 则更新最佳粒子适应度 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 超出范围适应度为 0 else 否则适应度按目标函数求解 re 1 2 pi sigma 2 exp x 2 y 2 2 sigma 2 end end 16 Log Polar 就是所谓的极坐标变换 通过变换能将笛卡尔坐标系中的图像映射到极坐标系中 确切的来说我这里不算是 Log Polar 因为 Log Polar 是将图像从 x y 映射到 log r theta 而我是将图像从 x y 映射到 r theta 原理是一样的 变换公式是 x r cos theta y r sin theta 其中 r 代表极径 theta 代表极角 这些应该学过高中的都会的 下图能形象的表示变换 精品文档 39欢迎下载 似乎有通过此变换求图像特征的算法 不过具体操作我也不清楚 我这里只是实现了变换的操作 原图如下 以图像中心作为极坐标原点进行极坐标变换 matlab 代码如下 精品文档 40欢迎下载 clear all close all clc img double imread lena jpg imshow img m n size img ox oy ginput 获得极坐标变换的原点 oy round oy ox round ox 求中心点到图像四个角的距离 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 求中心点距离四角距离的最大值 作为变换后图像的高 这个最大值也是极坐标变换的极径 radius round max up left up right down left down right angle 360 变换后图像的宽 imgn zeros radius angle for i 1 radius 纵坐标代表极径 不同情况不一样 for j 1 angle 横坐标代表极角 恒为 360 oy ox 作为极坐标变换中心坐标 需要作为偏移量相加 h oy round i sin j pi 180 w ox round i cos j pi 180 if h 0 最邻近插值 因为前面求 h w 用了 round end end end figure imshow imgn 17 通过本篇和上一篇的结合 应该就能做出拉普拉斯图像融合了 这里用的方法很简单 就是用模板和两个图像相乘 然后对处理后的两个图像再相加就可以了 精品文档 41欢迎下载 拉普拉斯融合就是对金字塔的每一层图像做这样的操作 然后再重构就行了 先看看这里的效果 原图 apple 原图 orange 模板图像 实际处理时需要对模板进行模糊 精品文档 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 过去研究过高斯金字塔 不过今天看来 当时似乎搞错了 金字塔上下层应该是采样的关系 而不是缩放的关系 而且不同层的 sigma 应该也有所不同 拉普拉斯金字塔在图像融合中有所应用 方法是首先对两个待融合图像求拉普拉斯残差金字塔 然后用模板对每一级残差图像进行融合得到融合后图像的残差金字塔 然后对这个金字塔进行 重构就能得到最终的融合图像 图像各尺度细节得到保留 注 融合时模板一般会先用高斯 函数模糊一下 不过这里不实现融合 只实现拉普拉斯金字塔的建立 建立金字塔的方框图如下 摘自 数字图像处理 第二版 通过 j 1 级近似图像和 j 级残差图像就能恢复 j 级输入图像 近似滤波器就是拉普拉斯滤波模板 插入滤波器是对拉普拉斯滤波模板采样形成的新滤波模板 使用的拉普拉斯滤波模板如下 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 产生的残差金字塔 精品文档 44欢迎下载 前面四张图像都是残差 最后一张图像是原图的 4 级近似 通过最后一张图和残差图 就能重 构回原图 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 级近似 end for i 5 1 2 调整图像大小 imgn i 1 imgn i 1 1 2 size imgn i 1 1 2 size imgn i 2 end for i 1 4 获得残差图像 i 级预测残差 精品文档 45欢迎下载 imgn i imgn i expand imgn i 1 w end for i 4 1 1 残差图像重构原图像 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 在图像拼接中有所使用 有时候也在图像理解的相关算法中有所使用 算法简介如下 摘自 图像处理 分析与机器视觉 第 3 版 精品文档 46欢迎下载 1 假设我们要将 n 个数据点 X x1 x1 xn 拟合为一个由至少 m 个点决定的模型 m n 对 于直线 m 2 我这里实际是两个不同均值 协方差高斯分布产生的数据 2 设迭代计数 k 1 3 从 X 中随机选取 m 个项并拟合一个模型 我这里直线拟合 选了 2 个项 4 给定偏差 计算 X 中相对于模型的残差在偏差 的个数 如果元素个数大于一个阈值 t 根 据一致点集重新拟合模型 可以利用最小二乘或其变种 算法终止 我这里的偏差为 1 阈值为数据个数的 2 3 5 设 k k 1 如果 k 小于一个事先给定的 K 跳至第 3 步 否则采用具有迄今最大的一致点集模 型 或算法失败 运行效果如下 红圈是所有的数据 蓝叉是符合拟合模型的数据 matlab 代码如下 main m close all clear all clc 真实的数据 mu 0 0 均值 S 1 2 5 2 5 8 协方差 data1 mvnrnd mu S 200 产生 200 个高斯分布数据 噪声数据 mu 2 2 S 8 0 0 8 data2 mvnrnd mu S 100 产生 100 个噪声数据 精品文档 47欢迎下载 data data1 data2 plot data 1 data 2 ro 显示全部数据 K 100 设置最大迭代次数 sigma 1 设置拟合直线与数据距离的偏差 pretotal 0 符合拟合模型的数据的个数 k 1 while pretotal size data 1 2 3 产生两个随机索引 找样本用 floor 向下取整 samp1