2015-04-02 14:34:07 gggg_ggg 阅读数 2862
  • 机器学习入门30天实战

    购买课程后,可扫码进入学习群,获取唐宇迪老师答疑 系列课程包含Python机器学习库,机器学习经典算法原理推导,基于真实数据集案例实战3大模块。从入门开始进行机器学习原理推导,以通俗易懂为基础形象解读晦涩难懂的机器学习算法工作原理,案例实战中使用Python工具库从数据预处理开始一步步完成整个建模工作!具体内容涉及Python必备机器学习库、线性回归算法原理推导、Python实现逻辑回归与梯度下降、案例实战,信用卡欺诈检测、决策树与集成算法、支持向量机原理推导、SVM实例与贝叶斯算法、机器学习常规套路与Xgboost算法、神经网络。

    7865 人正在学习 去看看 唐宇迪


    在数字图像处理过程中,经常会遇到求梯度后,重新构建图像的问题。一般情况下,都是通过解泊松方程(还有其他方式重构图像,具体算法如下图所示,),利用拉普拉斯算子求解;但有一点请注意泊松方程求出的只是近似值,无法求出精确的原始值。


                                                                  常用图像重建算法


2010-09-20 10:37:00 fengbingchun 阅读数 5720
  • 机器学习入门30天实战

    购买课程后,可扫码进入学习群,获取唐宇迪老师答疑 系列课程包含Python机器学习库,机器学习经典算法原理推导,基于真实数据集案例实战3大模块。从入门开始进行机器学习原理推导,以通俗易懂为基础形象解读晦涩难懂的机器学习算法工作原理,案例实战中使用Python工具库从数据预处理开始一步步完成整个建模工作!具体内容涉及Python必备机器学习库、线性回归算法原理推导、Python实现逻辑回归与梯度下降、案例实战,信用卡欺诈检测、决策树与集成算法、支持向量机原理推导、SVM实例与贝叶斯算法、机器学习常规套路与Xgboost算法、神经网络。

    7865 人正在学习 去看看 唐宇迪

转自:http://blog.sina.com.cn/s/blog_4b9b714a0100c9f7.html

 

梯度、边缘和角点

Sobel

使用扩展 Sobel 算子计算一阶、二阶、三阶或混合图像差分

 

void cvSobel( const CvArr* src, CvArr* dst, int xorder, int yorder, int aperture_size=3 );
src
输入图像.
dst
输出图像.
xorder
x 方向上的差分阶数
yorder
y 方向上的差分阶数
aperture_size
扩展 Sobel 核的大小,必须是 1, 3, 5 或 7。 除了尺寸为 1, 其它情况下, aperture_size ×aperture_size 可分离内核将用来计算差分。对 aperture_size=1的情况, 使用 3x1 或 1x3 内核(不进行高斯平滑操作)。这里有一个特殊变量 CV_SCHARR (=-1),对应 3x3 Scharr 滤波器,可以给出比 3x3 Sobel 滤波更精确的结果。Scharr 滤波器系数是:

/begin{bmatrix} -3 & 0 & 3 // -10 & 0 & 10 // -3 & 0 & 3 /end{bmatrix}

对 x-方向 以及转置矩阵对 y-方向。

函数 cvSobel 通过对图像用相应的内核进行卷积操作来计算图像差分:

dst(x,y) = /frac{d^{xorder+yorder}src} {dx^{xorder} dy^{yorder}} |(x,y)

由于Sobel 算子结合了 Gaussian 平滑和微分,所以,其结果或多或少对噪声有一定的鲁棒性。通常情况,函数调用采用如下参数 (xorder=1, yorder=0, aperture_size=3) 或 (xorder=0, yorder=1, aperture_size=3) 来计算一阶 x- 或 y- 方向的图像差分。第一种情况对应:

/begin{bmatrix} -1 & 0 & 1 // -2 & 0 & 2 // -1 & 0 & 1 /end{bmatrix} 核。

第二种对应:

/begin{bmatrix} -1 & -2 & -1 // 0 & 0 & 0 // 1 & 2 & 1 /end{bmatrix}

或者

/begin{bmatrix} 1 & 2 & 1 // 0 & 0 & 0 // -1 & -2 & -1 /end{bmatrix}

核的选则依赖于图像原点的定义 (origin 来自 IplImage 结构的定义)。由于该函数不进行图像尺度变换,所以和输入图像(数组)相比,输出图像(数组)的元素通常具有更大的绝对数值(译者注:即象素的深度)。为防止溢出,当输入图像是 8 位的,要求输出图像是 16 位的。当然可以用函数 cvConvertScale 或 cvConvertScaleAbs 转换为 8 位的。除了 8-比特 图像,函数也接受 32-位浮点数图像。所有输入和输出图像都必须是单通道的,并且具有相同的图像尺寸或者ROI尺寸。

 

Laplace

计算图像的 Laplacian 变换

void cvLaplace( const CvArr* src, CvArr* dst, int aperture_size=3 );
src
输入图像.
dst
输出图像.
aperture_size
核大小 (与 cvSobel 中定义一样).

函数 cvLaplace 计算输入图像的 Laplacian变换,方法是先用 sobel 算子计算二阶 x- 和 y- 差分,再求和:

dst(x,y) = d2src/dx2 + d2src/dy2

对 aperture_size=1 则给出最快计算结果,相当于对图像采用如下内核做卷积:

/begin{bmatrix} 0 & 1 & 0 // 1 & -4 & 1 // 0 & 1 & 0 /end{bmatrix}

类似于 cvSobel 函数,该函数也不作图像的尺度变换,所支持的输入、输出图像类型的组合和cvSobel一致。

 

Canny

采用Canny算法做边缘检测

 

void cvCanny( const CvArr* image, CvArr* edges, double threshold1, double threshold2, int aperture_size=3 );
image
输入图像.
edges
输出的边缘图像
threshold1
第一个阈值
threshold2
第二个阈值
aperture_size
Sobel 算子内核大小 (见 cvSobel).

函数 cvCanny 采用 CANNY 算法发现输入图像的边缘而且在输出图像中标识这些边缘。threshold1和threshold2 当中的

小阈值用来控制边缘连接,大的阈值用来控制强边缘的初始分割。

 

PreCornerDetect

计算用于角点检测的特征图,

void cvPreCornerDetect( const CvArr* image, CvArr* corners, int aperture_size=3 );
image
输入图像.
corners
保存候选角点的特征图
aperture_size
Sobel 算子的核大小(见cvSobel).

函数 cvPreCornerDetect 计算函数D_x^2D_{yy}+D_y^2D_{xx} - 2D_xD_yD_{xy} 其中 D_{/cdot} 表示一阶图像差分,D_{/cdot /cdot} 表示二阶图像差分。 角点被认为是函数的局部最大值:

// 假设图像格式为浮点数 IplImage* corners = cvCloneImage(image); IplImage* dilated_corners = cvCloneImage(image); 
IplImage* corner_mask = cvCreateImage( cvGetSize(image), 8, 1 ); cvPreCornerDetect( image, corners, 3 ); 
cvDilate( corners, dilated_corners, 0, 1 ); cvSubS( corners, dilated_corners, corners ); 
cvCmpS( corners, 0, corner_mask, CV_CMP_GE ); 
cvReleaseImage( &corners ); cvReleaseImage( &dilated_corners );

CornerEigenValsAndVecs

计算图像块的特征值和特征向量,用于角点检测
void cvCornerEigenValsAndVecs( const CvArr* image, CvArr* eigenvv, int block_size, int aperture_size=3 );
image
输入图像.
eigenvv
保存结果的数组。必须比输入图像宽 6 倍。
block_size
邻域大小 (见讨论).
aperture_size
Sobel 算子的核尺寸(见 cvSobel).
对每个象素,函数 cvCornerEigenValsAndVecs 考虑 block_size × block_size 大小的邻域 S(p),
然后在邻域上计算图像差分的相关矩阵:
M=/begin{bmatrix} /sum_{S(p)}(/frac {dI}{dx})^2 & /sum_{S(p)}(/frac{dI}{dx} /cdot /frac{dI}{dy}) // /sum_{S(p)}(/frac{dI}{dx} /cdot /frac{dI}{dy}) & /sum_{S(p)}(/frac {dI}{dy})^2 /end{bmatrix}

然后它计算矩阵的特征值和特征向量,并且按如下方式(λ1, λ2, x1, y1, x2, y2)存储这些值到输出图像中,其中
λ1, λ2 - M 的特征值,没有排序
(x1, y1) - 特征向量,对 λ1
(x2, y2) - 特征向量,对 λ2

CornerMinEigenVal

计算梯度矩阵的最小特征值,用于角点检测
void cvCornerMinEigenVal( const CvArr* image, CvArr* eigenval, int block_size, int aperture_size=3 );
image
输入图像.
eigenval
保存最小特征值的图像. 与输入图像大小一致
block_size
邻域大小 (见讨论 cvCornerEigenValsAndVecs).
aperture_size
Sobel 算子的核尺寸(见 cvSobel). 当输入图像是浮点数格式时,该参数表示用来计算差分固定的浮点滤波器的个数.
函数 cvCornerMinEigenVal 与 cvCornerEigenValsAndVecs 类似,但是它仅仅计算和存储每个象素点差分相关矩阵的最小特征值,
即前一个函数的 min(λ1, λ2)

CornerHarris

哈里斯(Harris)角点检测
void cvCornerHarris( const CvArr* image, CvArr* harris_responce, int block_size, int aperture_size=3, double k=0.04 );
image
输入图像。
harris_responce
存储哈里斯(Harris)检测responces的图像。与输入图像等大。
block_size
邻域大小(见关于cvCornerEigenValsAndVecs的讨论)。
aperture_size
扩展 Sobel 核的大小(见 cvSobel)。格式. 当输入图像是浮点数格式时,该参数表示用来计算差分固定的浮点滤波器的个数。
k
Harris detector free parameter. See the formula below.
harris 检测器的自由参数。请看如下公式。
The function cvCornerHarris runs the Harris edge detector on image. Similarly to cvCornerMinEigenVal and
cvCornerEigenValsAndVecs, for each pixel it calculates 2x2 gradient covariation matrix M
over block_size×block_size neighborhood. Then, it stores
det(M) - k*trace(M)2 to the destination image. Corners in the image can be found as local maxima of the
 destination image.
函数 cvCornerHarris 对输入图像进行 Harris 边界检测。类似于 cvCornerMinEigenVal 和 cvCornerEigenValsAndVecs。
对每个像素,在 block_size*block_size 大小的邻域上,计算其2*2梯度共变矩阵(或相关异变矩阵)M。
然后,将 det(M) - k*trace(M)2 (此公式有待考证,最后的“2”是否应为平方符号?这里2应该是平方)保存到输出图像中。
输入图像中的角点在输出图像中由局部最大值表示。

FindCornerSubPix

精确角点位置

void cvFindCornerSubPix( const CvArr* image, CvPoint2D32f* corners, int count, CvSize win, 
                                         CvSize zero_zone, CvTermCriteria criteria );
image
输入图像.
corners
输入角点的初始坐标,也存储精确的输出坐标
count
角点数目
win
搜索窗口的一半尺寸。如果 win=(5,5) 那么使用 5*2+1 × 5*2+1 = 11 × 11 大小的搜索窗口
zero_zone
死区的一半尺寸,死区为不对搜索区的中央位置做求和运算的区域。它是用来避免自相关矩阵出现的某些可能的奇异性。
当值为 (-1,-1) 表示没有死区。
criteria
求角点的迭代过程的终止条件。即角点位置的确定,要么迭代数大于某个设定值,或者是精确度达到某个设定值。
criteria 可以是最大迭代数目,或者是设定的精确度,也可以是它们的组合。

函数 cvFindCornerSubPix 通过迭代来发现具有子象素精度的角点位置,或如图所示的放射鞍点(radial saddle points)。

Image:Cornersubpix.png

子象素级角点定位的实现是基于对向量正交性的观测而实现的,即从中央点q到其邻域点p 的向量和p点处的

图像梯度正交(服从图像和测量噪声)。考虑以下的表达式:

εi=DIpiT•(q-pi)

其中,DIpi表示在q的一个邻域点pi处的图像梯度,q的值通过最小化εi得到。通过将εi设为0,可以建立系统方程如下:

sumi(DIpi•DIpiT)•q - sumi(DIpi•DIpiT•pi) = 0

其中q的邻域(搜索窗)中的梯度被累加。调用第一个梯度参数G和第二个梯度参数b,得到:

q=G-1•b

该算法将搜索窗的中心设为新的中心q,然后迭代,直到找到低于某个阈值点的中心位置。

 

GoodFeaturesToTrack

确定图像的强角点

void cvGoodFeaturesToTrack( const CvArr* image, CvArr* eig_image, CvArr* temp_image, 
CvPoint2D32f* corners, int* corner_count, double quality_level, double min_distance, const CvArr* mask=NULL );
image
输入图像,8-位或浮点32-比特,单通道
eig_image
临时浮点32-位图像,尺寸与输入图像一致
temp_image
另外一个临时图像,格式与尺寸与 eig_image 一致
corners
输出参数,检测到的角点
corner_count
输出参数,检测到的角点数目
quality_level
最大最小特征值的乘法因子。定义可接受图像角点的最小质量因子。
min_distance
限制因子。得到的角点的最小距离。使用 Euclidian 距离
mask
ROI:感兴趣区域。函数在ROI中计算角点,如果 mask 为 NULL,则选择整个图像。必须为单通道的灰度图,大小与输入图像相同。
mask对应的点不为0,表示计算该点。

函数 cvGoodFeaturesToTrack 在图像中寻找具有大特征值的角点。该函数,首先用cvCornerMinEigenVal

计算输入图像的每一个象素点的最小特征值,

并将结果存储到变量 eig_image 中。然后进行非最大值抑制(仅保留3x3邻域中的局部最大值)。

下一步将最小特征值小于 quality_level•max(eig_image(x,y)) 排除掉。最后,函数确保所有发现的角点之间具有足够的距离,

(最强的角点第一个保留,然后检查新的角点与已有角点之间的距离大于 min_distance )。

2015-02-11 18:57:48 gggg_ggg 阅读数 6192
  • 机器学习入门30天实战

    购买课程后,可扫码进入学习群,获取唐宇迪老师答疑 系列课程包含Python机器学习库,机器学习经典算法原理推导,基于真实数据集案例实战3大模块。从入门开始进行机器学习原理推导,以通俗易懂为基础形象解读晦涩难懂的机器学习算法工作原理,案例实战中使用Python工具库从数据预处理开始一步步完成整个建模工作!具体内容涉及Python必备机器学习库、线性回归算法原理推导、Python实现逻辑回归与梯度下降、案例实战,信用卡欺诈检测、决策树与集成算法、支持向量机原理推导、SVM实例与贝叶斯算法、机器学习常规套路与Xgboost算法、神经网络。

    7865 人正在学习 去看看 唐宇迪


对于有些人,看这些枯燥的公式符号是件痛苦的事情;但痛苦后总会有所欣喜,如果你充分利用它的话,你更能体会到他的美妙;先来几张效果图,激发你学习数学的欲望:

                  注释:图像融合效果,分别应用了不同的算法

在图像图形处理中, 梯度、散度和旋度 有很重要的作用,比如图像修复中的解泊松方程,目标跟踪等等,可以说是他们无处不在。

来句废话:可能有些人,对于数学符号里面倒三角 正三角 符号的意思?与读法感到迷惑,现稍作解释;

△一般指拉普拉斯算子

▽读Nabla,奈不拉,也可以读作“Del” 这是场论中的符号,是矢量微分算符.高等数学中的梯度,散度,旋度都会用到这个算符.其二阶导数中旋度的散度又称Laplace算符


继续核心内容:

梯度、散度和旋度是矢量分析里的重要概念。之所以是“分析”,因为三者是三种偏导数计算形式。这里假设读者已经了解了三者的定义。它们的符号分别记作如下:

梯度、散度和旋度                                 

                                 

                                

从符号中可以获得这样的信息:

求梯度是针对一个标量函数,求梯度的结果是得到一个矢量函数。这里φ称为势函数

求散度则是针对一个矢量函数,得到的结果是一个标量函数,跟求梯度是反一下的;

求旋度是针对一个矢量函数,得到的还是一个矢量函数。

这三种关系可以从定义式很直观地看出,因此可以求“梯度的散度”、“散度的梯度”、“梯度的旋度”、“旋度的散度”和“旋度的旋度”,只有旋度可以连续作用两次,而一维波动方程具有如下的形式

梯度、散度和旋度                               (1)

其中a为一实数,于是可以设想,对于一个矢量函数来说,要求得它的波动方程,只有求它的“旋度的旋度”才能得到。下面先给出梯度、散度和旋度的计算式:

梯度、散度和旋度                         (2)

梯度、散度和旋度                          (3)

梯度、散度和旋度                           (4)

旋度公式略显复杂。这里结合麦克斯韦电磁场理论,来讨论前面几个“X度的X度”。

 

I.梯度的散度:

根据麦克斯韦方程有:

梯度、散度和旋度                                  

梯度、散度和旋度                                  (5)

则电势的梯度的散度为

梯度、散度和旋度                            

这是一个三维空间上的标量函数,常记作

梯度、散度和旋度                                 (6)

称为泊松方程,而算符2称为拉普拉斯算符。事实上因为定义

梯度、散度和旋度                                

所以有

梯度、散度和旋度                               

当然,这只是一种记忆方式。

当空间内无电荷分布时,即ρ=0,则称为拉普拉斯方程

梯度、散度和旋度                                      

当我们仅需要考虑一维情况时,比如电荷均匀分布的无限大平行板电容器之间(不包含极板)的电场,我们知道该电场只有一个指向,场强处处相等,于是该电场满足一维拉普拉斯方程,即

梯度、散度和旋度                                     

这就是说如果那边平行板电容器的负极板接地,则板间一点处的电压与该点距负极板的距离呈线性关系。

 

II.散度的梯度:

散度的梯度,从上面的公式中可以看到结果会比较复杂,但是它的物理意义却是很明确的,因为从麦克斯韦方程可以看出空间某点处电场的散度是该点处的电荷密度,那么再求梯度就是空间中电荷密度的梯度。这就好比说清水中滴入一滴红墨水,起初水面红色浓度最高,杯底浓度最低,这样水面与杯底形成一个浓度梯度,红墨水由水面向杯底扩散,最后均匀。在半导体中,载流子分布的不均匀会导致扩散电流。

散度的梯度这个概念其实不常用,因为计算复杂,但在后面讲用它来推导一个矢量恒等式。

 

III.梯度的旋度:

对于梯度的旋度,直接把(2)式代入(4)式中,有

梯度、散度和旋度

由于势函数在空间一点的领域内往往是有二阶连续混合偏导数的,因此上式的结果为0.所以说梯度的旋度为零,它的物理意义也是很明确的。

比如一个人从海平面爬到一座山上,无论它是从山的陡坡爬上去还是从缓坡爬上去,亦或者坐直升机上去,重力对他所做的功总是相等的,即力场的做工只与位移有关而与路径无关,这样的场称为保守场,而保守场是无旋场。再比如绘有等高线的地图,如果某点只有一个一根等高线穿过,那么该点有一个确定的相对高度。如果该点有两条或以上的等高线穿过,则这个点处在悬崖边上,这个点处是不可微,也就没有求梯度的意义。

 

IV.旋度的散度:

求旋度的散度也是将(4)式代入(3)式即可。若令

梯度、散度和旋度                            (7)

梯度、散度和旋度                                      

                                

                                

从而

梯度、散度和旋度                                   

                             

                             

将上面三式相加结果也为零。所以说旋度的散度为零,这就意味着一个散度场任意叠加上一个有旋场不会改变其散度,也就是说光凭矢量场的散度无法唯一地确定这个矢量场。而光凭矢量场的旋度也无法唯一地确定这个矢量,这是因为有旋场可以叠加上这么一个矢量场而不改变其旋度,而这个矢量场是一个标量函数的梯度。

 

V.旋度的旋度:

旋度的旋度将是本文的重点。若所研究的空间范围内是无源的,即ρ=0J=0,则根据麦克斯韦方程有:

梯度、散度和旋度                                (8)

梯度、散度和旋度                             (9)

梯度、散度和旋度                                   (10)

梯度、散度和旋度                               (11)

(9)式两端取旋度

梯度、散度和旋度                          (12)

再将(8)式代入(12)式有

梯度、散度和旋度                             (13)

看到这里容易让人想到式(1),前面说式(1)的方程为一维波动方程,那么跟(13)式有什么联系呢?棘手的问题是算旋度已经够复杂了,算旋度的旋度岂不是更费周折?幸好有矢量恒等式可以利用来帮助简化计算,这里要用到前面所讲的散度的梯度。即有:

梯度、散度和旋度                          (14)

这里拉普拉斯算子作用于一个矢量函数时,意义变得不明确了,它和前面的几个“X度的X度”都不一样,实际上它有这样的定义:

梯度、散度和旋度                        (15)

为了验证式(14)还是要对计算“旋度的旋度”,但以后可以直接利用该式。还是做(7)式那样的处理,即令

梯度、散度和旋度                                

梯度、散度和旋度                                    

                              

                              

于是

梯度、散度和旋度               (16)

而令

梯度、散度和旋度                                     

梯度、散度和旋度             (17)

两式相减有

梯度、散度和旋度               (18)

类似地有

梯度、散度和旋度                                    

梯度、散度和旋度                                    

由于所关心的空间内是无源的,所以式(13)变成

梯度、散度和旋度                             (19)

这个方程很重要,称为三维波动方程,这也从理论上揭示了电磁波的存在。它的各分量展开后比较复杂,实际上我们无法绘制出一个向四面八方传播的波的振动图像,但好在可以画出一维和二维的波,从而了解波的性质。有些事物我们无法在现实世界中呈现,或绘制出图形,但是数学上却可以计算且有确切的物理意义,比如高于三维的空间,不得不感叹数学的神奇,感叹我们生活的世界的神奇。

 

VI.几个矢量恒等式:

前面已经介绍了一个矢量恒等式,还有其他几个重要的恒等式。由于三种“度”是三种不同微分算法,虽然有些场合可以把▽当做一个普通的矢量来处理,但并不总是正确的,这一点需要引起注意。

梯度、散度和旋度

 

梯度、散度和旋度

 

梯度、散度和旋度

 

这里“×”乘的优先级高于“·”乘对于普通三个不共面的矢量ABC则有A·B×C=C·A×B=B·C×A。得到的结果是令三个矢量共起点,以三个矢量的模为棱构成的六面体的体积或它的负值。但是对于▽算子,则一般

   梯度、散度和旋度                                

但是一般有

梯度、散度和旋度                                

实际上上面的矢量恒等式就是上式的扩展

梯度、散度和旋度

梯度、散度和旋度

上两式相减有

梯度、散度和旋度

记忆上式的方法是记住下标的顺序是xyzyzxzxy

 

梯度、散度和旋度

 

这个等式相对容易证明,但前提是要在直角坐标下。

梯度、散度和旋度

2012-06-14 23:05:48 jia20003 阅读数 47856
  • 机器学习入门30天实战

    购买课程后,可扫码进入学习群,获取唐宇迪老师答疑 系列课程包含Python机器学习库,机器学习经典算法原理推导,基于真实数据集案例实战3大模块。从入门开始进行机器学习原理推导,以通俗易懂为基础形象解读晦涩难懂的机器学习算法工作原理,案例实战中使用Python工具库从数据预处理开始一步步完成整个建模工作!具体内容涉及Python必备机器学习库、线性回归算法原理推导、Python实现逻辑回归与梯度下降、案例实战,信用卡欺诈检测、决策树与集成算法、支持向量机原理推导、SVM实例与贝叶斯算法、机器学习常规套路与Xgboost算法、神经网络。

    7865 人正在学习 去看看 唐宇迪

图像处理之图像梯度效果


基本思想:

利用X方向与Y方向分别实现一阶微分,求取振幅,实现图像梯度效果。关于如何计算图像

一阶微分参见这里:http://blog.csdn.net/jia20003/article/details/7562092

使用的两种微分算子分别为Prewitt与Sobel,其中Soble在X, Y两个方向算子分别为:


Prewitt在X, Y方向上梯度算子分别为:


二:程序思路及实现

梯度滤镜提供了两个参数:

– 方向,用来要决定图像完成X方向梯度计算, Y方向梯度计算,或者是振幅计算

– 算子类型,用来决定是使用sobel算子或者是prewitt算子。

计算振幅的公式可以参见以前《图像处理之一阶微分应用》的文章


 三:运行效果

原图像如下:


基于Prewitt与sobel算子的XY方向振幅效果如下:



该滤镜的源代码如下:

package com.process.blur.study;

import java.awt.image.BufferedImage;
/**
 * 
 * @author gloomy-fish
 * @date 2012-06-11
 * 
 * prewitt operator 
 * X-direction
 * -1, 0, 1
 * -1, 0, 1
 * -1, 0, 1
 * 
 * Y-direction
 * -1, -1, -1
 *  0,  0,  0
 *  1,  1,  1
 *  
 * sobel operator
 * X-direction
 * -1, 0, 1
 * -2, 0, 2
 * -1, 0, 1
 * 
 * Y-direction
 * -1, -2, -1
 *  0,  0,  0
 *  1,  2,  1
 *
 */
public class GradientFilter extends AbstractBufferedImageOp {

	// prewitt operator
	public final static int[][] PREWITT_X = new int[][]{{-1, 0, 1}, {-1, 0, 1}, {-1, 0, 1}};
	public final static int[][] PREWITT_Y = new int[][]{{-1, -1, -1}, {0,  0,  0}, {1,  1,  1}};
	
	// sobel operator
	public final static int[][] SOBEL_X = new int[][]{{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};
	public final static int[][] SOBEL_Y = new int[][]{{-1, -2, -1}, {0,  0,  0}, {1,  2,  1}};
	
	// direction parameter
	public final static int X_DIRECTION = 0;
	public final static int Y_DIRECTION = 2;
	public final static int XY_DIRECTION = 4;
	private int direction;
	private boolean isSobel;
	
	public GradientFilter() {
		direction = XY_DIRECTION;
		isSobel = true;
	}
	
	public void setSoble(boolean sobel) {
		this.isSobel = sobel;
	}

	public int getDirection() {
		return direction;
	}

	public void setDirection(int direction) {
		this.direction = direction;
	}
	
	@Override
	public BufferedImage filter(BufferedImage src, BufferedImage dest) {
		int width = src.getWidth();
        int height = src.getHeight();

        if (dest == null )
        	dest = createCompatibleDestImage( src, null );

        int[] inPixels = new int[width*height];
        int[] outPixels = new int[width*height];
        getRGB( src, 0, 0, width, height, inPixels );
        int index = 0, index2 = 0;
        double xred = 0, xgreen = 0, xblue = 0;
        double yred = 0, ygreen = 0, yblue = 0;
        int newRow, newCol;
        for(int row=0; row<height; row++) {
        	int ta = 255, tr = 0, tg = 0, tb = 0;
        	for(int col=0; col<width; col++) {
        		index = row * width + col;
        		for(int subrow = -1; subrow <= 1; subrow++) {
        			for(int subcol = -1; subcol <= 1; subcol++) {
        				newRow = row + subrow;
        				newCol = col + subcol;
        				if(newRow < 0 || newRow >= height) {
        					newRow = row;
        				}
        				if(newCol < 0 || newCol >= width) {
        					newCol = col;
        				}
        				index2 = newRow * width + newCol;
                        tr = (inPixels[index2] >> 16) & 0xff;
                        tg = (inPixels[index2] >> 8) & 0xff;
                        tb = inPixels[index2] & 0xff;
                        
                        if(isSobel) {
                        	xred += (SOBEL_X[subrow + 1][subcol + 1] * tr);
                        	xgreen +=(SOBEL_X[subrow + 1][subcol + 1] * tg);
                        	xblue +=(SOBEL_X[subrow + 1][subcol + 1] * tb);
                        	
                        	yred += (SOBEL_Y[subrow + 1][subcol + 1] * tr);
                        	ygreen +=(SOBEL_Y[subrow + 1][subcol + 1] * tg);
                        	yblue +=(SOBEL_Y[subrow + 1][subcol + 1] * tb);
                        } else {
                        	xred += (PREWITT_X[subrow + 1][subcol + 1] * tr);
                        	xgreen +=(PREWITT_X[subrow + 1][subcol + 1] * tg);
                        	xblue +=(PREWITT_X[subrow + 1][subcol + 1] * tb);
                        	
                        	yred += (PREWITT_Y[subrow + 1][subcol + 1] * tr);
                        	ygreen +=(PREWITT_Y[subrow + 1][subcol + 1] * tg);
                        	yblue +=(PREWITT_Y[subrow + 1][subcol + 1] * tb);
                        }
        			}
        		}
        		
                double mred = Math.sqrt(xred * xred + yred * yred);
                double mgreen = Math.sqrt(xgreen * xgreen + ygreen * ygreen);
                double mblue = Math.sqrt(xblue * xblue + yblue * yblue);
                if(XY_DIRECTION == direction) 
                {
                	outPixels[index] = (ta << 24) | (clamp((int)mred) << 16) | (clamp((int)mgreen) << 8) | clamp((int)mblue);
                } 
                else if(X_DIRECTION == direction)
                {
                	outPixels[index] = (ta << 24) | (clamp((int)yred) << 16) | (clamp((int)ygreen) << 8) | clamp((int)yblue);
                } 
                else if(Y_DIRECTION == direction) 
                {
                	outPixels[index] = (ta << 24) | (clamp((int)xred) << 16) | (clamp((int)xgreen) << 8) | clamp((int)xblue);
                } 
                else 
                {
                	// as default, always XY gradient
                	outPixels[index] = (ta << 24) | (clamp((int)mred) << 16) | (clamp((int)mgreen) << 8) | clamp((int)mblue);
                }
                
                // cleanup for next loop
                newRow = newCol = 0;
                xred = xgreen = xblue = 0;
                yred = ygreen = yblue = 0;
                
        	}
        }
        setRGB(dest, 0, 0, width, height, outPixels );
        return dest;
	}
	
	public static int clamp(int value) {
		return value < 0 ? 0 : (value > 255 ? 255 : value);
	}

}

转载文章请务必注明出自本博客!

精彩内容查看 - 《数字图像处理-空间域卷积视频教程》

2014-06-07 10:21:47 a8039974 阅读数 1261
  • 机器学习入门30天实战

    购买课程后,可扫码进入学习群,获取唐宇迪老师答疑 系列课程包含Python机器学习库,机器学习经典算法原理推导,基于真实数据集案例实战3大模块。从入门开始进行机器学习原理推导,以通俗易懂为基础形象解读晦涩难懂的机器学习算法工作原理,案例实战中使用Python工具库从数据预处理开始一步步完成整个建模工作!具体内容涉及Python必备机器学习库、线性回归算法原理推导、Python实现逻辑回归与梯度下降、案例实战,信用卡欺诈检测、决策树与集成算法、支持向量机原理推导、SVM实例与贝叶斯算法、机器学习常规套路与Xgboost算法、神经网络。

    7865 人正在学习 去看看 唐宇迪

形状特征

  (一)特点:各种基于形状特征的检索方法都可以比较有效地利用图像中感兴趣的目标来进行检索,但它们也有一些共同的问题,包括:①目前基于形状的检索方法还缺乏比较完善的数学模型;②如果目标有变形时检索结果往往不太可靠;③许多形状特征仅描述了目标局部的性质,要全面描述目标常对计算时间和存储量有较高的要求;④许多形状特征所反映的目标形状信息与人的直观感觉不完全一致,或者说,特征空间的相似性与人视觉系统感受到的相似性有差别。另外,从 2-D 图像中表现的 3-D 物体实际上只是物体在空间某一平面的投影,从 2-D 图像中反映出来的形状常不是 3-D 物体真实的形状,由于视点的变化,可能会产生各种失真。

 

(二)常用的特征提取与匹配方法

 

Ⅰ几种典型的形状特征描述方法

 

通常情况下,形状特征有两类表示方法,一类是轮廓特征,另一类是区域特征。图像的轮廓特征主要针对物体的外边界,而图像的区域特征则关系到整个形状区域。

 

几种典型的形状特征描述方法:

 

(1)边界特征法该方法通过对边界特征的描述来获取图像的形状参数。其中Hough 变换检测平行直线方法和边界方向直方图方法是经典方法。Hough 变换是利用图像全局特性而将边缘像素连接起来组成区域封闭边界的一种方法,其基本思想是点—线的对偶性;边界方向直方图法首先微分图像求得图像边缘,然后,做出关于边缘大小和方向的直方图,通常的方法是构造图像灰度梯度方向矩阵。

 

(2)傅里叶形状描述符法

 

傅里叶形状描述符(Fourier shape descriptors)基本思想是用物体边界的傅里叶变换作为形状描述,利用区域边界的封闭性和周期性,将二维问题转化为一维问题。

 

由边界点导出三种形状表达,分别是曲率函数、质心距离、复坐标函数。

 

(3)几何参数法

 

形状的表达和匹配采用更为简单的区域特征描述方法,例如采用有关形状定量测度(如矩、面积、周长等)的形状参数法(shape factor)。在 QBIC 系统中,便是利用圆度、偏心率、主轴方向和代数不变矩等几何参数,进行基于形状特征的图像检索。

 

需要说明的是,形状参数的提取,必须以图像处理及图像分割为前提,参数的准确性必然受到分割效果的影响,对分割效果很差的图像,形状参数甚至无法提取。

 

(4)形状不变矩法

 

利用目标所占区域的矩作为形状描述参数。

 

(5)其它方法

 

近年来,在形状的表示和匹配方面的工作还包括有限元法(Finite Element Method 或 FEM)、旋转函数(Turning Function)和小波描述符(Wavelet Descriptor)等方法。

 

实际上,只是提取物体的形状,这并不难,最难的是这些特征该怎么用!

特征嘛,自然是讲此物区分彼物的特点。

那么假如,给出了一系列不同形状物体的轮廓该如何识别出他们呢?正方形,圆形,矩形,椭圆,不规则图形,再进一步,这些图形由于受到信号的干扰,有噪声存在时,该如何去识别他们呢?

那就可以使用形状的特征了,我们定义一些参数,来描述这些形状。

1 矩形度:R = A0/A; A0为区域面积,A为区域最小外接矩形面积。

那么R = 1 时,为矩形的概率很大,R = PI/4时为圆的可能性最大

2 体态比 T = a/b;

a b 分别为区域最小外接矩形的长和宽。

T = 1 为正方形或者圆形,

T>1 为细长图形

3 球状性 S = Ri/Rc

Ri Rc分别为内切圆和外接圆半径,圆心都在中心上

4 球状性 C = Ur/Pr

Ur 为区域重心到轮廓点的平均距离

Pr 为区域重心到轮廓点的均方差

5 中心矩

这一特征,使用颇为频繁,OpenCV有专门的函数求解(p,q)次矩

6 长轴 短轴

最小外接矩形的长轴和短轴

7面积

一般会作为阈值使用,判定某个区域的面积在两个阈值之间才判定有效

 

 下面给出各个形状特征的求法:

  1. //图像的形状特征分析  
  2. #include <cv.h>  
  3. #include <cxcore.h>  
  4. #include <highgui.h>  
  5. #include <iostream>  
  6. using namespace std;  
  7.   
  8. int main()  
  9. {  
  10.     IplImage *src = cvLoadImage("E:\\image\\mapleleaf.tif",0);  
  11.     IplImage *image = cvCreateImage(cvGetSize(src),8,3);   
  12.     image = cvCloneImage(src);  
  13.     cvNamedWindow("src",1);  
  14.     cvNamedWindow("dst",1);  
  15.     cvShowImage("src",src);  
  16.   
  17.     CvMemStorage *storage = cvCreateMemStorage(0);  
  18.     CvSeq * seq = cvCreateSeq(CV_SEQ_ELTYPE_POINT, sizeof(CvSeq), sizeof(CvPoint), storage);  
  19.     CvSeq * tempSeq = cvCreateSeq(CV_SEQ_ELTYPE_POINT, sizeof(CvSeq), sizeof(CvPoint), storage);  
  20.     //新图,将轮廓绘制到dst  
  21.     IplImage *dst = cvCreateImage(cvGetSize(src),8,3);   
  22.     cvZero(dst);//赋值为0  
  23.     double length,area;  
  24.   
  25.     //获取轮廓  
  26.     int cnt = cvFindContours(src,storage,&seq);//返回轮廓的数目  
  27.     cout<<"number of contours   "<<cnt<<endl;  
  28.       
  29.     //计算边界序列的参数 长度 面积 矩形 最小矩形   
  30.     //并输出每个边界的参数  
  31.     CvRect rect;  
  32.     CvBox2D box;  
  33.     double axislong,axisShort;//长轴和短轴  
  34.     double temp1= 0.0,temp2 = 0.0;  
  35.     double Rectangle_degree;//矩形度  
  36.     double long2short;//体态比  
  37.     double x0,y0;  
  38.     long sumX  = 0 ,sumY = 0;  
  39.     double sum =0.0;  
  40.     int i,j,m,n;  
  41.     unsigned char* ptr;  
  42.   
  43.     double UR;//区域重心到轮廓的平均距离  
  44.     double PR;//区域重心到轮廓点的均方差  
  45.     CvPoint * contourPoint;   
  46.     int count = 0;  
  47.     double CDegree;//圆形性  
  48.   
  49.     CvPoint *A,*B,*C;  
  50.     double AB,BC,AC;  
  51.     double cosA,sinA;  
  52.     double tempR,inscribedR;  
  53.       
  54.     for (tempSeq = seq;tempSeq != NULL; tempSeq = tempSeq->h_next)  
  55.     {  
  56.         //tempSeq = seq->h_next;  
  57.         length = cvArcLength(tempSeq);  
  58.         area = cvContourArea(tempSeq);  
  59.         cout<<"Length = "<<length<<endl;  
  60.         cout<<"Area = "<<area<<endl;  
  61.         cout<<"num of point "<<tempSeq->total<<endl;  
  62.         //外接矩形  
  63.         rect = cvBoundingRect(tempSeq,1);  
  64.           
  65.   
  66.         //绘制轮廓和外接矩形  
  67.         cvDrawContours(dst,tempSeq,CV_RGB(255,0,0),CV_RGB(255,0,0),0);  
  68.         cvRectangleR(dst,rect,CV_RGB(0,255,0));  
  69.         cvShowImage("dst",dst);  
  70.         //cvWaitKey();  
  71.   
  72.         //绘制轮廓的最小外接圆   
  73.         CvPoint2D32f center;//亚像素精度 因此需要使用浮点数  
  74.         float radius;  
  75.         cvMinEnclosingCircle(tempSeq,¢er,&radius);  
  76.         cvCircle(dst,cvPointFrom32f(center),cvRound(radius),CV_RGB(100,100,100));  
  77.         cvShowImage("dst",dst);  
  78.         //cvWaitKey();  
  79.   
  80.         //寻找近似的拟合椭圆 可以使斜椭圆  
  81.         CvBox2D ellipse = cvFitEllipse2(tempSeq);  
  82.         cvEllipseBox(dst,ellipse,CV_RGB(255,255,0));  
  83.         cvShowImage("dst",dst);  
  84.         //cvWaitKey();  
  85.   
  86.   
  87.         //绘制外接最小矩形  
  88.         CvPoint2D32f pt[4];  
  89.         box = cvMinAreaRect2(tempSeq,0);  
  90.         cvBoxPoints(box,pt);  
  91.         for(int i = 0;i<4;++i){  
  92.             cvLine(dst,cvPointFrom32f(pt[i]),cvPointFrom32f(pt[((i+1)%4)?(i+1):0]),CV_RGB(0,0,255));  
  93.         }  
  94.         cvShowImage("dst",dst);  
  95.         //cvWaitKey();  
  96.   
  97.   
  98.         //下面开始分析图形的形状特征   
  99.         //长轴 短轴  
  100.         temp1 = sqrt(pow(pt[1].x -pt[0].x,2) + pow(pt[1].y -pt[0].y,2));  
  101.         temp2 = sqrt(pow(pt[2].x -pt[1].x,2) + pow(pt[2].y -pt[1].y,2));  
  102.   
  103.         if (temp1 > temp2)  
  104.         {  
  105.             axislong = temp1;  
  106.             axisShort=temp2;  
  107.         }   
  108.         else  
  109.         {  
  110.             axislong = temp2;  
  111.             axisShort=temp1;  
  112.         }  
  113.           
  114.         cout<<"long axis: "<<axislong<<endl;  
  115.         cout<<"short axis: "<<axisShort<<endl;  
  116.         //矩形度 轮廓面积和最小外接矩形面积(可以是斜矩形)之比  
  117.         Rectangle_degree = (double)area/(axisShort*axislong);  
  118.   
  119.         cout<<"Rectangle degree :"<<Rectangle_degree<<endl;  
  120.         //体态比or长宽比 最下外接矩形的长轴和短轴的比值  
  121.         long2short = axislong/axisShort;  
  122.         cout<<"ratio of long axis to short axis: "<<long2short<<endl;  
  123.         //球状性 由于轮廓的内切圆暂时无法求出先搁置   
  124.         //先求内切圆半径 枚举任意轮廓上的三个点,半径最小的就是内切圆的半径  
  125.         //以下的最大内切圆半径求法有误 待改进  
  126.           
  127.         /* 
  128.         for (int i = 0 ; i< tempSeq->total -2;i++) 
  129.         { 
  130.             for (int j= i+1; j<tempSeq->total-1;j++) 
  131.             { 
  132.                 for (int m = j+1; m< tempSeq->total; m++) 
  133.                 { 
  134.                     //已知圆上三点,求半径 
  135.                     A = (CvPoint*)cvGetSeqElem(tempSeq ,i); 
  136.                     B = (CvPoint*)cvGetSeqElem(tempSeq ,j); 
  137.                     C = (CvPoint*)cvGetSeqElem(tempSeq,m); 
  138.                     AB = sqrt(pow((double)A->x - B->x,2)+ pow((double)A->y - B->y,2)); 
  139.                     AC =sqrt(pow((double)A->x - C->x,2) + pow((double)A->y - C->y,2)); 
  140.                     BC = sqrt(pow((double)B->x - C->x,2)+ pow((double)B->y - C->y,2)); 
  141.  
  142.                     cosA = ((B->x - A->x)*(C->x - A->x) + (B->y - A->y)*(C->y - A->y))/(AB*AC); 
  143.                     sinA = sqrt(1 - pow(cosA,2)); 
  144.                     tempR = BC/(2*sinA); 
  145.  
  146.                     if (m == 2) 
  147.                     { 
  148.                         inscribedR = tempR; 
  149.                     }  
  150.                     else 
  151.                     { 
  152.                         if (tempR < inscribedR) 
  153.                         { 
  154.                             inscribedR = tempR; 
  155.                         } 
  156.                     } 
  157.  
  158.  
  159.                 } 
  160.             } 
  161.         } 
  162.  
  163.         //输出最大内切圆半径 
  164.         cout<<"radius of max inscribed circle  "<<inscribedR<<endl; 
  165.         */  
  166.         //圆形性 假设轮廓内是实心的  
  167.         //球区域中心x0 y0  
  168.         sumX = 0;  
  169.         sumY = 0;  
  170.         src = cvCloneImage(image);  
  171.         for (int i = 0 ; i< src->height;i++)  
  172.         {  
  173.             for (int j = 0; j< src->width;j++)  
  174.             {  
  175.                 ptr = (unsigned char *)src->imageData + i*src->widthStep + j;  
  176.                 if ((*ptr) > 128)  
  177.                 {  
  178.                     sumX += (long)j;  
  179.                     sumY += (long)i;  
  180.                 }  
  181.   
  182.             }  
  183.         }  
  184.         x0 = sumX/area;  
  185.         y0 = sumY/area;  
  186.         cout<<"center of gravity "<<x0<<"   "<<y0<<endl;  
  187.      //求区域到重心的平均距离  
  188.         sum = 0;  
  189.         count = 0;  
  190.         for (m = 0 ; m< tempSeq->total;m++)  
  191.         {  
  192.             contourPoint = (CvPoint*)cvGetSeqElem(tempSeq,m);  
  193.             sum += sqrt(pow(contourPoint->x - x0,2)+ pow(contourPoint->y - y0,2));  
  194.             count++;  
  195.         }  
  196.         UR = sum/count;  
  197.         cout<<"mean distance to center of gravity"<<UR<<endl;  
  198.         //求区域重心到轮廓点的均方差  
  199.         sum = 0;  
  200.         for (m = 0 ; m< tempSeq->total;m++)  
  201.         {  
  202.             contourPoint = (CvPoint*)cvGetSeqElem(tempSeq,m);  
  203.             temp1 = sqrt(pow(contourPoint->x - x0,2)+ pow(contourPoint->y - y0,2));  
  204.             sum += pow(temp1 - UR,2);  
  205.         }  
  206.         PR = sum/count;  
  207.         cout<<"mean square error of distance to center of gravity"<<PR<<endl;  
  208.         //圆形性   
  209.         CDegree= UR/PR;  
  210.         cout<<"degree of circle "<<CDegree<<endl;  
  211.         //中心距  
  212.           
  213.         cvWaitKey(0);  
  214.     }  
  215.       
  216.       
  217.     cvReleaseImage(&src);  
  218.     cvReleaseImage(&dst);  
  219.     cvReleaseMemStorage(&storage);    
  220.       
  221.     return 0;  
  222. }  



 

没有更多推荐了,返回首页