圖像處理之Harris角度檢測(cè)算法

Harris角度檢測(cè)是通過數(shù)學(xué)計(jì)算在圖像上發(fā)現(xiàn)角度特征的一種算法，而且其具有旋轉(zhuǎn)不

變性的特質(zhì)。OpenCV中的Shi-Tomasi角度檢測(cè)就是基于Harris角度檢測(cè)改進(jìn)算法。

基本原理：

角度是一幅圖像上最明顯與重要的特征，對(duì)于一階導(dǎo)數(shù)而言，角度在各個(gè)方向的變化是

最大的，而邊緣區(qū)域在只是某一方向有明顯變化。一個(gè)直觀的圖示如下：

數(shù)學(xué)原理：

基本數(shù)學(xué)公式如下：

其中W(x, y)表示移動(dòng)窗口，I(x, y)表示像素灰度值強(qiáng)度，范圍為0～255。根據(jù)泰勒級(jí)數(shù)

計(jì)算一階到N階的偏導(dǎo)數(shù)，最終得到一個(gè)Harris矩陣公式：

根據(jù)Harris的矩陣計(jì)算矩陣特征值,然后計(jì)算Harris角度響應(yīng)值：

其中K為系數(shù)值，通常取值范圍為0.04 ~ 0.06之間。

算法詳細(xì)步驟

第一步：計(jì)算圖像X方向與Y方向的一階高斯偏導(dǎo)數(shù)Ix與Iy

第二步：根據(jù)第一步結(jié)果得到Ix^2 , Iy^2與Ix*Iy值

第三步：高斯模糊第二步三個(gè)值得到Sxx, Syy, Sxy

第四部：定義每個(gè)像素的Harris矩陣，計(jì)算出矩陣的兩個(gè)特質(zhì)值

第五步：計(jì)算出每個(gè)像素的R值

第六步：使用3X3或者5X5的窗口，實(shí)現(xiàn)非最大值壓制

第七步：根據(jù)角度檢測(cè)結(jié)果計(jì)算，最提取到的關(guān)鍵點(diǎn)以綠色標(biāo)記，顯示在原圖上。

程序關(guān)鍵代碼解讀：

第一步計(jì)算一階高斯偏導(dǎo)數(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);

關(guān)于如何計(jì)算高斯一階與二階偏導(dǎo)數(shù)請(qǐng)看這里：

http://blog.csdn.net/jia20003/article/details/16369143

http://blog.csdn.net/jia20003/article/details/7664777

第三步：分別對(duì)第二步計(jì)算出來的三個(gè)值，單獨(dú)進(jìn)行高斯

模糊計(jì)算，代碼如下：

	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;         	}         }		 	}

第六步：非最大信號(hào)壓制(non-max value suppression)

這個(gè)在邊源檢測(cè)中是為了得到一個(gè)像素寬的邊緣，在這里則

是為了得到準(zhǔn)確的一個(gè)角點(diǎn)像素，去掉非角點(diǎn)值。代碼如下：

	/*** 	 * 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);         		}         	}         } 		 	}

運(yùn)行效果：

程序完整源代碼：

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); 		} 		 		// 基于高斯方法，中心點(diǎn)化窗口計(jì)算一階導(dǎo)數(shù)和，關(guān)鍵一步 SumIx2, SumIy2 and SumIxIy, 高斯模糊 		calculateGaussianBlur(width, height);  		// 求取Harris Matrix 特征值  		// 計(jì)算角度相應(yīng)值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);         		}         	}         } 		 	} 	 	/*** 	 * 計(jì)算兩個(gè)特征值，然后得到R，公式如下，可以自己推導(dǎo)，關(guān)于怎么計(jì)算矩陣特征值，請(qǐng)看這里： 	 * 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);         	}         } 	}  }