
图像处理之Harris角度检测算法
Harris角度检测是通过数学计算在图像上发现角度特征的一种算法,而且其具有旋转不
变性的特质。OpenCV中的Shi-Tomasi角度检测就是基于Harris角度检测改进算法。
基本原理:
角度是一幅图像上最明显与重要的特征,对于一阶导数而言,角度在各个方向的变化是
最大的,而边缘区域在只是某一方向有明显变化。一个直观的图示如下:
数学原理:
基本数学公式如下:
其中W(x, y)表示移动窗口,I(x, y)表示像素灰度值强度,范围为0~255。根据泰勒级数
计算一阶到N阶的偏导数,最终得到一个Harris矩阵公式:
根据Harris的矩阵计算矩阵特征值,然后计算Harris角度响应值:
其中K为系数值,通常取值范围为0.04 ~ 0.06之间。
算法详细步骤
第一步:计算图像X方向与Y方向的一阶高斯偏导数Ix与Iy
第二步:根据第一步结果得到Ix^2 , Iy^2与Ix*Iy值
第三步:高斯模糊第二步三个值得到Sxx, Syy, Sxy
第四部:定义每个像素的Harris矩阵,计算出矩阵的两个特质值
第五步:计算出每个像素的R值
第六步:使用3X3或者5X5的窗口,实现非最大值压制
第七步:根据角度检测结果计算,最提取到的关键点以绿色标记,显示在原图上。
程序关键代码解读:
第一步计算一阶高斯偏导数的Ix与Iy值代码如下:
filter.setDirectionType(GaussianDerivativeFilter.X_DIRECTION); BufferedImage xImage = filter.filter(grayImage, null); getRGB( xImage, 0, 0, width, height, inPixels ); extractPixelData(inPixels, GaussianDerivativeFilter.X_DIRECTION, height, width); filter.setDirectionType(GaussianDerivativeFilter.Y_DIRECTION); BufferedImage yImage = filter.filter(grayImage, null); getRGB( yImage, 0, 0, width, height, inPixels ); extractPixelData(inPixels, GaussianDerivativeFilter.Y_DIRECTION, height, width);
关于如何计算高斯一阶与二阶偏导数请看这里:
http://blog.csdn.net/jia20003/article/details/16369143
http://blog.csdn.net/jia20003/article/details/7664777
第三步:分别对第二步计算出来的三个值,单独进行高斯
模糊计算,代码如下:
private void calculateGaussianBlur(int width, int height) { int index = 0; int radius = (int)window_radius; double[][] gw = get2DKernalData(radius, sigma); double sumxx = 0, sumyy = 0, sumxy = 0; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { for(int subrow =-radius; subrow<=radius; subrow++) { for(int subcol=-radius; subcol<=radius; subcol++) { int nrow = row + subrow; int ncol = col + subcol; if(nrow >= height || nrow < 0) { nrow = 0; } if(ncol >= width || ncol < 0) { ncol = 0; } int index2 = nrow * width + ncol; HarrisMatrix whm = harrisMatrixList.get(index2); sumxx += (gw[subrow + radius][subcol + radius] * whm.getXGradient()); sumyy += (gw[subrow + radius][subcol + radius] * whm.getYGradient()); sumxy += (gw[subrow + radius][subcol + radius] * whm.getIxIy()); } } index = row * width + col; HarrisMatrix hm = harrisMatrixList.get(index); hm.setXGradient(sumxx); hm.setYGradient(sumyy); hm.setIxIy(sumxy); // clean up for next loop sumxx = 0; sumyy = 0; sumxy = 0; } } }
第六步:非最大信号压制(non-max value suppression)
这个在边源检测中是为了得到一个像素宽的边缘,在这里则
是为了得到准确的一个角点像素,去掉非角点值。代码如下:
/*** * we still use the 3*3 windows to complete the non-max response value suppression */ private void nonMaxValueSuppression(int width, int height) { int index = 0; int radius = (int)window_radius; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { index = row * width + col; HarrisMatrix hm = harrisMatrixList.get(index); double maxR = hm.getR(); boolean isMaxR = true; for(int subrow =-radius; subrow<=radius; subrow++) { for(int subcol=-radius; subcol<=radius; subcol++) { int nrow = row + subrow; int ncol = col + subcol; if(nrow >= height || nrow < 0) { nrow = 0; } if(ncol >= width || ncol < 0) { ncol = 0; } int index2 = nrow * width + ncol; HarrisMatrix hmr = harrisMatrixList.get(index2); if(hmr.getR() > maxR) { isMaxR = false; } } } if(isMaxR) { hm.setMax(maxR); } } } }
运行效果:
程序完整源代码:
package com.gloomyfish.image.harris.corner; import java.awt.image.BufferedImage; import java.util.ArrayList; import java.util.List; import com.gloomyfish.filter.study.GrayFilter; public class HarrisCornerDetector extends GrayFilter { private GaussianDerivativeFilter filter; private List<HarrisMatrix> harrisMatrixList; private double lambda = 0.04; // scope : 0.04 ~ 0.06 // i hard code the window size just keep it' size is same as // first order derivation Gaussian window size private double sigma = 1; // always private double window_radius = 1; // always public HarrisCornerDetector() { filter = new GaussianDerivativeFilter(); harrisMatrixList = new ArrayList<HarrisMatrix>(); } @Override public BufferedImage filter(BufferedImage src, BufferedImage dest) { int width = src.getWidth(); int height = src.getHeight(); initSettings(height, width); if ( dest == null ) dest = createCompatibleDestImage( src, null ); BufferedImage grayImage = super.filter(src, null); int[] inPixels = new int[width*height]; // first step - Gaussian first-order Derivatives (3 × 3) - X - gradient, (3 × 3) - Y - gradient filter.setDirectionType(GaussianDerivativeFilter.X_DIRECTION); BufferedImage xImage = filter.filter(grayImage, null); getRGB( xImage, 0, 0, width, height, inPixels ); extractPixelData(inPixels, GaussianDerivativeFilter.X_DIRECTION, height, width); filter.setDirectionType(GaussianDerivativeFilter.Y_DIRECTION); BufferedImage yImage = filter.filter(grayImage, null); getRGB( yImage, 0, 0, width, height, inPixels ); extractPixelData(inPixels, GaussianDerivativeFilter.Y_DIRECTION, height, width); // second step - calculate the Ix^2, Iy^2 and Ix^Iy for(HarrisMatrix hm : harrisMatrixList) { double Ix = hm.getXGradient(); double Iy = hm.getYGradient(); hm.setIxIy(Ix * Iy); hm.setXGradient(Ix*Ix); hm.setYGradient(Iy*Iy); } // 基于高斯方法,中心点化窗口计算一阶导数和,关键一步 SumIx2, SumIy2 and SumIxIy, 高斯模糊 calculateGaussianBlur(width, height); // 求取Harris Matrix 特征值 // 计算角度相应值R R= Det(H) - lambda * (Trace(H))^2 harrisResponse(width, height); // based on R, compute non-max suppression nonMaxValueSuppression(width, height); // match result to original image and highlight the key points int[] outPixels = matchToImage(width, height, src); // return result image setRGB( dest, 0, 0, width, height, outPixels ); return dest; } private int[] matchToImage(int width, int height, BufferedImage src) { int[] inPixels = new int[width*height]; int[] outPixels = new int[width*height]; getRGB( src, 0, 0, width, height, inPixels ); int index = 0; for(int row=0; row<height; row++) { int ta = 0, tr = 0, tg = 0, tb = 0; for(int col=0; col<width; col++) { index = row * width + col; ta = (inPixels[index] >> 24) & 0xff; tr = (inPixels[index] >> 16) & 0xff; tg = (inPixels[index] >> 8) & 0xff; tb = inPixels[index] & 0xff; HarrisMatrix hm = harrisMatrixList.get(index); if(hm.getMax() > 0) { tr = 0; tg = 255; // make it as green for corner key pointers tb = 0; outPixels[index] = (ta << 24) | (tr << 16) | (tg << 8) | tb; } else { outPixels[index] = (ta << 24) | (tr << 16) | (tg << 8) | tb; } } } return outPixels; } /*** * we still use the 3*3 windows to complete the non-max response value suppression */ private void nonMaxValueSuppression(int width, int height) { int index = 0; int radius = (int)window_radius; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { index = row * width + col; HarrisMatrix hm = harrisMatrixList.get(index); double maxR = hm.getR(); boolean isMaxR = true; for(int subrow =-radius; subrow<=radius; subrow++) { for(int subcol=-radius; subcol<=radius; subcol++) { int nrow = row + subrow; int ncol = col + subcol; if(nrow >= height || nrow < 0) { nrow = 0; } if(ncol >= width || ncol < 0) { ncol = 0; } int index2 = nrow * width + ncol; HarrisMatrix hmr = harrisMatrixList.get(index2); if(hmr.getR() > maxR) { isMaxR = false; } } } if(isMaxR) { hm.setMax(maxR); } } } } /*** * 计算两个特征值,然后得到R,公式如下,可以自己推导,关于怎么计算矩阵特征值,请看这里: * http://www.sosmath.com/matrix/eigen1/eigen1.html * * A = Sxx; * B = Syy; * C = Sxy*Sxy*4; * lambda = 0.04; * H = (A*B - C) - lambda*(A+B)^2; * * @param width * @param height */ private void harrisResponse(int width, int height) { int index = 0; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { index = row * width + col; HarrisMatrix hm = harrisMatrixList.get(index); double c = hm.getIxIy() * hm.getIxIy(); double ab = hm.getXGradient() * hm.getYGradient(); double aplusb = hm.getXGradient() + hm.getYGradient(); double response = (ab -c) - lambda * Math.pow(aplusb, 2); hm.setR(response); } } } private void calculateGaussianBlur(int width, int height) { int index = 0; int radius = (int)window_radius; double[][] gw = get2DKernalData(radius, sigma); double sumxx = 0, sumyy = 0, sumxy = 0; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { for(int subrow =-radius; subrow<=radius; subrow++) { for(int subcol=-radius; subcol<=radius; subcol++) { int nrow = row + subrow; int ncol = col + subcol; if(nrow >= height || nrow < 0) { nrow = 0; } if(ncol >= width || ncol < 0) { ncol = 0; } int index2 = nrow * width + ncol; HarrisMatrix whm = harrisMatrixList.get(index2); sumxx += (gw[subrow + radius][subcol + radius] * whm.getXGradient()); sumyy += (gw[subrow + radius][subcol + radius] * whm.getYGradient()); sumxy += (gw[subrow + radius][subcol + radius] * whm.getIxIy()); } } index = row * width + col; HarrisMatrix hm = harrisMatrixList.get(index); hm.setXGradient(sumxx); hm.setYGradient(sumyy); hm.setIxIy(sumxy); // clean up for next loop sumxx = 0; sumyy = 0; sumxy = 0; } } } public double[][] get2DKernalData(int n, double sigma) { int size = 2*n +1; double sigma22 = 2*sigma*sigma; double sigma22PI = Math.PI * sigma22; double[][] kernalData = new double[size][size]; int row = 0; for(int i=-n; i<=n; i++) { int column = 0; for(int j=-n; j<=n; j++) { double xDistance = i*i; double yDistance = j*j; kernalData[row][column] = Math.exp(-(xDistance + yDistance)/sigma22)/sigma22PI; column++; } row++; } // for(int i=0; i<size; i++) { // for(int j=0; j<size; j++) { // System.out.print("\t" + kernalData[i][j]); // } // System.out.println(); // System.out.println("\t ---------------------------"); // } return kernalData; } private void extractPixelData(int[] pixels, int type, int height, int width) { int index = 0; for(int row=0; row<height; row++) { int ta = 0, tr = 0, tg = 0, tb = 0; for(int col=0; col<width; col++) { index = row * width + col; ta = (pixels[index] >> 24) & 0xff; tr = (pixels[index] >> 16) & 0xff; tg = (pixels[index] >> 8) & 0xff; tb = pixels[index] & 0xff; HarrisMatrix matrix = harrisMatrixList.get(index); if(type == GaussianDerivativeFilter.X_DIRECTION) { matrix.setXGradient(tr); } if(type == GaussianDerivativeFilter.Y_DIRECTION) { matrix.setYGradient(tr); } } } } private void initSettings(int height, int width) { int index = 0; for(int row=0; row<height; row++) { for(int col=0; col<width; col++) { index = row * width + col; HarrisMatrix matrix = new HarrisMatrix(); harrisMatrixList.add(index, matrix); } } } }
数据分析咨询请扫描二维码
若不方便扫码,搜微信号:CDAshujufenxi
SASEM 决策树:理论与实践应用 在复杂的决策场景中,如何从海量数据中提取有效信息并制定科学决策,是各界关注的焦点。SASEM 决 ...
2025-07-30SPSS 语法使用详解 在当今数据驱动的时代,SPSS( Statistical Package for the Social Sciences)作为一款功能强大的统计分析软 ...
2025-07-30人工智能对CDA数据分析领域的影响 人工智能对 CDA(Certified Data Analyst,注册数据分析师)数据分析领域的影响是全方位、多层 ...
2025-07-30MySQL执行计划中rows的计算逻辑:从原理到实践 MySQL 执行计划中 rows 的计算逻辑:从原理到实践 在 MySQL 数据库的查询优化中 ...
2025-07-29左偏态分布转正态分布:方法、原理与实践 左偏态分布转正态分布:方法、原理与实践 在统计分析、数据建模和科学研究中,正态分 ...
2025-07-29CDA 数据分析师的职业生涯规划:从入门到卓越的成长之路 在数字经济蓬勃发展的当下,数据已成为企业核心竞争力的重要来源,而 CD ...
2025-07-29CDA数据分析师证书考取全攻略 一、了解 CDA 数据分析师认证 CDA 数据分析师认证是一套科学化、专业化、国际化的人才考核标准, ...
2025-07-29解析神经网络中 Softmax 函数的核心作用 在神经网络的发展历程中,激活函数扮演着至关重要的角色,它们为网络赋予了非线性能力, ...
2025-07-29解析 response.text 与 response.content 的核心区别 在网络数据请求与处理的场景中,开发者经常需要从服务器返回的响应中提取数 ...
2025-07-29鸢尾花判别分析:机器学习中的经典实践案例 在机器学习的世界里,有一个经典的数据集如同引路明灯,为无数初学者打开了模式识别 ...
2025-07-29用 Python 开启数据分析之旅:从基础到实践的完整指南 在数据驱动决策的时代,数据分析已成为各行业不可或缺的核心能力。而 Pyt ...
2025-07-29从 CDA LEVEL II 考试题型看 Python 数据分析要点 在数据科学领域蓬勃发展的当下,CDA(Certified Data Analyst)认证成为众多从 ...
2025-07-29CDA 数据分析师的工作范围解析 在数字化时代的浪潮下,数据已成为企业发展的核心资产之一。CDA(Certified Data Analyst)数据分 ...
2025-07-29解析 insert into select 是否会锁表:原理、场景与应对策略 在数据库操作中,insert into select 是一种常用的批量数据插入语句 ...
2025-07-29用 Power BI 制作地图热力图:基于经纬度数据的实践指南 在数据可视化领域,地图热力图凭借直观呈现地理数据分布密度的优势,成 ...
2025-07-29从数据到决策:CDA 数据分析师如何重塑职场竞争力与行业价值 在数字经济席卷全球的今天,数据已从 “辅助工具” 升级为 “核心资 ...
2025-07-292025 年 CDA 数据分析师考纲焕新,引领行业人才新标准 在数字化浪潮奔涌向前的当下,数据已成为驱动各行业发展的核心要素。作为 ...
2025-07-29PyTorch 核心机制:损失函数与反向传播如何驱动模型进化 在深度学习的世界里,模型从 “一无所知” 到 “精准预测” 的蜕变,离 ...
2025-07-29t 检验与 Wilcoxon 检验:数据差异分析的两大核心方法 在数据分析的广阔领域中,判断两组或多组数据之间是否存在显著差异是一项 ...
2025-07-29PowerBI 添加索引列全攻略 在使用 PowerBI 进行数据处理与分析时,添加索引列是一项极为实用的操作技巧。索引列能为数据表中的每 ...
2025-07-29