随笔分类 - Matlab练习程序
摘要：RANSAC在图像拼接中有所使用，有时候也在图像理解的相关算法中有所使用。算法简介如下(摘自《图像处理、分析与机器视觉（第3版）》)：1.假设我们要将n个数据点X={x1,x1,...,xn}拟合为一个由至少m个点决定的模型（m<=n，对于直线，m=2）。（我这里实际是两个不同均值、协方差高斯分布产生的数据）2.设迭代计数k=1。3.从X中随机选取m个项并拟合一个模型。（我这里直线拟合，选了2个项）4.给定偏差ε，计算X中相对于模型的残差在偏差ε的个数，如果元素个数大于一个阈值t，根据一致点集重新拟合模型（可以利用最小二乘或其变种），算法终止。（我这里的偏差为1，阈值为数据个数的2/3）
阅读全文
摘要：曾经想要实现过Bertalmio图像修复算法，无奈自身实力不够，耗费两天时间也没能实现。昨天博客上有人问到TV模型，这个模型我过去是没听说过的，于是就找来相关论文研究了一下，发现TV模型也可以用来修复图像，于是就有了想实现的想法。用到的偏微分方程技巧和各项异性扩散很像。先看看效果吧：原lena:随手截的噪声图：合成的需要修复的图：修复后的图（没有处理边界）：对于从来没有接触过图像修复的我来说，效果真是惊艳了。下面介绍运算步骤：和各项异性扩散类似，整个算法也是基于迭代的，迭代公式如下：其中Io代表当前处理的像素，Ip代表邻域像素，p就可以取news四邻域。H定义如下：这里的lambda为自定义的
阅读全文
摘要：HOG（Histogram of Oriented Gradient）方向梯度直方图，主要用来提取图像特征，最常用的是结合svm进行行人检测。算法流程图如下（这篇论文上的）：下面我再结合自己的程序，表述一遍吧：1.对原图像gamma校正，img=sqrt(img);2.求图像竖直边缘，水平边缘，边缘强度，边缘斜率。3.将图像每16*16（取其他也可以）个像素分到一个cell中。对于256*256的lena来说，就分成了16*16个cell了。4.对于每个cell求其梯度方向直方图。通常取9（取其他也可以）个方向（特征），也就是每360/9=40度分到一个方向，方向大小按像素边缘强度加权。最后归
阅读全文
摘要：radon变换就是图像在不同方向上的投影。下图f(x,y)可以代表图像，R(x')就是图像向右下方的投影。数学上是按投影方向进行线积分，在图像领域就是按照投影方向累加像素就行了。matlab中有radon函数，为了熟悉原理我又写了一个。通过旋转图像映射，效率不怎么样。canny后Lena:45度投影：代码如下：main.mclear all; close all; clc;img=imread('lena.jpg');ed=edge(img,'canny');r=Radon(ed,45); %检测直线什么的，可以投影到0-179度上%r=radon(ed
阅读全文
摘要：Renyi熵是对通常的香农熵的扩展，算是q阶广义熵。公式如下：其中P和香农熵公式中的P一样，是概率。当q=1时公式退化为香农熵公式。（如何证明？看wiki吧）有用此公式寻找图像最佳二值化阈值的。首先定义前景区域A，背景区域B。那么前景与背景区域像素相应的Renyi熵就如下定义：其中k是当前取的灰度级，P(A)是像素在A区域的概率，P(B)类似。当然，这里说区域，不是指空间区域，是像素灰度级区域。最后图像Renyi熵求最佳阈值定义如下：这里得到的K就是分割阈值。分割效果如下：原图：分割后：在编程时还需要确定阶数q，我取的是2。matlab代码如下：clear all;close all;clc;
阅读全文
摘要：普通的直方图就是统计图像中像素灰度出现的次数。二维直方图的其中一维就是普通的直方图，即统计图像中像素灰度出现的次数。另一维的概念就比较模糊了，我看很多论文中都是用的当前像素邻域均值灰度出现的次数。按照这个意思，其实用当前像素邻域中值灰度出现的次数也是可以的。或者用什么厉害的算法处理一下，在统计出现灰度的次数也是可以的。所以按照这样推广，其实构成N维直方图也是未尝不可的。至于有什么用处嘛，我还需要再研究研究。lena的二维直方图，1维是普通直方图，2维是半径为1的邻域均值直方图：matlab代码如下：clear all;close all;clc;img=imread('lena.jpg
阅读全文
摘要：中值滤波是很经典的算法了。今天看论文又知道还有一种叫自适应中值滤波的算法RAMF。原论文在这里。RAMF主要通过以下两步来处理图像。1.首先确定最大的滤波半径，然后用一个合适的半径r对图像进行滤波。计算当前滤波半径像素灰度的Imin,Imax,Imed，然后判断Imed是否在[Imin,Imax]中间，如果在则向下进行，否则扩大当前半径r继续滤波直到r等于最大滤波半径。2.如果当前处理的像素img(i,j)在[Imin,Imax]之间，则输出当前像素，否则输出当前滤波半径中值像素Imed。看下效果吧：噪声图像：RAMF算法：普通3*3中值滤波：matlab代码如下：clear all;clos
阅读全文
摘要：关于纹理合成最经典的论文应该就数Efros的Texture Synthesis by Non-parametric Sampling这篇论文了，引用量近2000。这里的合成是基于样例的，就是先有一个小的纹理图像，然后合成一个大的。我个人的理解是，先产生一个大的随机图像，然后对随机图像中每个像素及其邻域像素和小的纹理图像做对比，找出随机图像当前邻域和小的纹理图像最接近的邻域，然后将小的纹理图像最接近的像素赋值给大的随机图像。这样每处理一个随机图像中的像素，就需要遍历小的纹理图像的所有像素，所以运行时间是很慢的。先看下效果：小的纹理图像：合成的纹理图像：matlab代码如下：clear all;c
阅读全文
摘要：所谓的最终腐蚀并不是把图像不断腐蚀直到一片黑为止，那样还有什么意思呢。最终腐蚀的意思是连续腐蚀过程中连通子区域消失前所有剩余部分的并。看下图有一个连通区域：那么最终腐蚀就是把图像不断腐蚀到只剩圆环最中间的白色区域为止。如下：所以一般最终腐蚀之后，如果再腐蚀一下，那么图像就全没了。最终腐蚀的通常算法会先对原图像使用距离变换，然后求区域极大值。不过这种方法怎么说呢，速度倒是快，不过我实验了很多次，求得的结果很不稳定。所以我就自己开发了一个算法，虽然有些慢，不过很稳定。1.首先对原图像进行连通区域标记。2.对标记后图像不断进行腐蚀，每腐蚀一次，重新标记连通区域。3.每腐蚀并且标记一次之后，检查腐蚀后
阅读全文
摘要：顾名思义，条件膨胀是带有限定条件的膨胀。限定条件通常有两个。1.膨胀始终在原始集合内。如下图：原始集合图像：集合内部的子集图像：如果对后一个图进行规则1的条件膨胀，那么最后的结果就会得到第一幅图。所以这里第一幅图即是原始集合图像又是最终条件膨胀图像。2.膨胀之后各自粒子互不连接。如下图：原图：按规则2条件膨胀需首先标记连通区域：然后再进行互补侵犯的膨胀：根据具体应用环境，可以使用单一的条件或两个条件同时使用。分水岭分割的最后筑坝那一步就是两个条件同时使用到了。单独规则1的matlab代码：clear all;close all;clc;mask=imread('p1.jpg')
阅读全文
摘要：目标内边界的像素全都在目标里面，目标外边界的像素全都不在目标上，是包围着目标的。二值图像内外边界的计算都是有两种方法的，所以一共是4种算法，不过实际用到跟踪的只有一个而已。首先是内边界跟踪：第一种方法不是跟踪方法。步骤是先对原图像腐蚀，然后用原图像减去腐蚀后的图像就得到边界了。第二种方法是跟踪方法。步骤如下：1.遍历图像。2.标记第一个遇见像素块的前景像素（i,j）。3.对这个像素周围八邻域逆时针搜索，如果搜索到周围有前景像素，那么更新坐标（i,j）为（i',j'），并标记。4.不断执行第3步直到再次遇见此像素块第一次标记的像素。5.继续执行第1步。然后是外边界跟踪：第一种方法
阅读全文
摘要：聚类算法，不是分类算法。分类算法是给一个数据，然后判断这个数据属于已分好的类中的具体哪一类。聚类算法是给一大堆原始数据，然后通过算法将其中具有相似特征的数据聚为一类。这里的k-means聚类，是事先给出原始数据所含的类数，然后将含有相似特征的数据聚为一个类中。所有资料中还是Andrew Ng介绍的明白。首先给出原始数据{x1,x2,...,xn}，这些数据没有被标记的。初始化k个随机数据u1,u2,...,uk。这些xn和uk都是向量。根据下面两个公式迭代就能求出最终所有的u,这些u就是最终所有类的中心位置。公式一：意思就是求出所有数据和初始化的随机数据的距离，然后找出距离每个初始数据最近的数
阅读全文
摘要：图像区域分裂可以将图像中相似的区域归为一个子区域，否则这个区域分裂。分裂的步骤还是很好理解的。对一个图像而言，首先确定一个分裂准则，如果当前图像区域中有像素没有满足这个准则，那么本区域继续分裂，否则分裂停止。我这里用的准则就是当前像素灰度是否和本区域平均灰度差距过大，如果相差很大则继续分裂，否则分裂停止并将本区域赋值为本区域平均值。分裂具体步骤：1.输入图像作为初始区域。2.判断本区域是否有像素灰度与本区域平均灰度差距过大。若有，则本区域分解为四部分，对每部分继续执行步骤2；若无，结束。（显然的递归嘛，当然，递归和迭代一般是能互换的）先看看分裂效果：matlab代码如下：main.mclear
阅读全文
摘要：主要是用来平滑图像的，克服了高斯模糊的缺陷，各向异性扩散在平滑图像时是保留图像边缘的（和双边滤波很像）。通常我们有将图像看作矩阵的，看作图的，看作随机过程的，记得过去还有看作力场的。这次新鲜，将图像看作热量场了。每个像素看作热流，根据当前像素和周围像素的关系，来确定是否要向周围扩散。比如某个邻域像素和当前像素差别较大，则代表这个邻域像素很可能是个边界，那么当前像素就不向这个方向扩散了，这个边界也就得到保留了。先看下效果吧：具体的推导公式都是热学上的，自己也不太熟悉，感兴趣的可以去看原论文，引用量超7000吶。我这里只介绍一下最终结论用到的公式。主要迭代方程如下：I就是图像了，因为是个迭代公式，
阅读全文
摘要：自动聚焦需要评判图像的模糊或者清晰程度。当然，真正聚焦还需要配合硬件，我现在已经不搞硬件了，所以也就不去关心那一方面了。主要有三种评判方法：1.灰度方差法评判图像灰度的变化程度，我先使用了全局灰度方差，不过效果不好。又实验了局部灰度方差，不过由于运算量的问题，也没做出什么结果。2.边缘评判法这个很简单了，通过对图像进行边缘检测，比如sobel，prewitt，roberts等等各种边缘检测吧，然后找出边缘最丰富的那一张图片。3.能量法图像细节越丰富，傅里叶分解后的频率肯定越多，所以能量也能作为一种评判标准。我这里的实验图像序列在这里可以下到。实验结果，图中红色是sobel评判，绿色是二阶拉普拉
阅读全文
摘要：关于小波变换我只是有一个很朴素了理解。不过小波变换可以和傅里叶变换结合起来理解。傅里叶变换是用一系列不同频率的正余弦函数去分解原函数，变换后得到是原函数在正余弦不同频率下的系数。小波变换使用一系列的不同尺度的小波去分解原函数，变换后得到的是原函数在不同尺度小波下的系数。不同的小波通过平移与尺度变换分解，平移是为了得到原函数的时间特性，尺度变换是为了得到原函数的频率特性。小波变换步骤：1.把小波w(t)和原函数f(t)的开始部分进行比较，计算系数C。系数C表示该部分函数与小波的相似程度。2.把小波向右移k单位，得到小波w(t-k)，重复1。重复该部知道函数f结束.3.扩展小波w(t)，得到小波w
阅读全文
摘要：计算Haar特征时有用到这个。wiki介绍的非常好，我可耻的转了：积分图的每一点（x,y）的值是原图中对应位置的左上角区域的所有值得和：而且，积分图可以只遍历一次图像即可有效的计算出来，因为积分图每一点的（x,y）值是：一旦积分图计算完毕，对任意矩形区域的和的计算就可以在常数时间内完成。如下图中，阴影矩形区域的值：clear all;close all;clc;img=double(imread('lena.jpg'));[m n]=size(img);%计算积分图I=zeros(m,n);for i=1:m for j=1:n if i==1 && j==1
阅读全文
摘要：很多算法都用到了这个矩阵，比如Hessian affine region detector、SURF，虽然这些算法我还没有完全搞透，不过那都是后话，先把这个矩阵搞出来再说，学习是不断迭代的过程。整个矩阵的形成主要由四部分公式决定：g(x,y)就是高斯函数了，没啥可说的：对高斯函数的x和y分别求二阶偏导：求出的模板对原图进行卷积：卷积后的值构成Hessian矩阵：所以这里的H是相当于原图像每一个像素都能求出一个Hessian矩阵。matlab代码如下：clear all;close all;clc;img=double(imread('lena.jpg'));[m n]=size
阅读全文
摘要：今天在看局部熵方面的内容，看论文中介绍的内容感觉局部熵挺容易了，于是就有了实现的想法，结果效果非常糟糕。 得到的几乎是一张空白的图片，就像下面一样： 究其原因是各种论文上都写了这样一个公式： 这里f(i,j)就是在m*n这个局部的像素，这个没问题，不过这里的p是什么东西，按这里的定义p是当前像素灰度
阅读全文
摘要：算法步骤：1.任意找一顶点加入树中。2.寻找所有与树相邻的元素，并取其边权重最小的并且不在树中的元素入树。3.继续第二步，直到所有元素都入树。效果和Kruskal算法是类似的。matlab代码：main.mclear all;close all;clc;G=[0 4 0 0 0 0 0 8 0; 4 0 8 0 0 0 0 11 0; 0 8 0 7 0 4 0 0 2; 0 0 7 0 9 14 0 0 0; 0 0 0 9 0 10 0 0 0; 0 0 4 14 10 0 2 0 0; 0 0 0 0 0 2 0 1 6; 8 11 0 0 0 0 1 0 7...
阅读全文