圖像處理之Harris角度檢測算法
Harris角度檢測是通過數(shù)學計算在圖像上發(fā)現(xiàn)角度特征的一種算法,而且其具有旋轉(zhuǎn)不
變性的特質(zhì)。OpenCV中的Shi-Tomasi角度檢測就是基于Harris角度檢測改進算法。
基本原理:
角度是一幅圖像上最明顯與重要的特征,對于一階導數(shù)而言,角度在各個方向的變化是
最大的,而邊緣區(qū)域在只是某一方向有明顯變化。一個直觀的圖示如下:
![](http://image68.360doc.com/DownloadImg/2014/01/0620/38063854_1)
數(shù)學原理:
基本數(shù)學公式如下:
![](http://image68.360doc.com/DownloadImg/2014/01/0620/38063854_2)
其中W(x, y)表示移動窗口,I(x, y)表示像素灰度值強度,范圍為0~255。根據(jù)泰勒級數(shù)
計算一階到N階的偏導數(shù),最終得到一個Harris矩陣公式:
![](http://image68.360doc.com/DownloadImg/2014/01/0620/38063854_3)
根據(jù)Harris的矩陣計算矩陣特征值 ,然后計算Harris角度響應值:
![](http://image68.360doc.com/DownloadImg/2014/01/0620/38063854_5)
其中K為系數(shù)值,通常取值范圍為0.04 ~ 0.06之間。
算法詳細步驟
第一步:計算圖像X方向與Y方向的一階高斯偏導數(shù)Ix與Iy
第二步:根據(jù)第一步結果得到Ix^2 , Iy^2與Ix*Iy值
第三步:高斯模糊第二步三個值得到Sxx, Syy, Sxy
第四部:定義每個像素的Harris矩陣,計算出矩陣的兩個特質(zhì)值
第五步:計算出每個像素的R值
第六步:使用3X3或者5X5的窗口,實現(xiàn)非最大值壓制
第七步:根據(jù)角度檢測結果計算,最提取到的關鍵點以綠色標記,顯示在原圖上。
程序關鍵代碼解讀:
第一步計算一階高斯偏導數(shù)的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);
關于如何計算高斯一階與二階偏導數(shù)請看這里:
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);
- }
- }
- }
-
- }
運行效果:
![](http://image68.360doc.com/DownloadImg/2014/01/0620/38063854_6)
程序完整源代碼:
- 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);
- }
-
- // 基于高斯方法,中心點化窗口計算一階導數(shù)和,關鍵一步 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./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);
- }
- }
- }
-
- }
最后注意:
我是把彩色圖像變?yōu)榛叶葓D像來計算,這個計算量小點
處理容易點,此外很多圖像處理軟件都會用標記來顯示
關鍵點像素,我沒有,只是將關鍵點像素改為綠色。
所以可以從這些方面有很大的提高空間。
|