harris corner detector java实现

/**

 * Millie : Multifunctional Library For Image Processing

 * 

 * (c) Copyright 2009 by Humbert Florent

 * 

 *      This program is free software; you can redistribute it and/or modify  

 *      it under the terms of the GNU General Public License as published by  

 *      the Free Software Foundation; only version 2 of the License.          

 *                                                                            

 *      This program is distributed in the hope that it will be useful,       

 *      but WITHOUT ANY WARRANTY; without even the implied warranty of        

 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         

 *      GNU General Public License for more details.                          

 *                                                                            

 *      You should have received a copy of the GNU General Public License     

 *      along with this program; if not, write to the Free Software           

 *      Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA             

 *      02111-1307, USA.                                                      

 */

package millie.operator.detection;

import java.util.ArrayList;

import java.util.Iterator;

import java.util.List;

import millie.image.Image;

import millie.operator.commons.component.ComponentSeparableOp;

/**

 * Harris Corner Detector

 * 

 *  k = det(A) - lambda * trace(A)^2

 * 

 *  Where A is the second-moment matrix 

 * 

 *            | Lx(x+dx,y+dy)    Lx.Ly(x+dx,y+dy) |

 *  A =  Sum  |                                    | * Gaussian(dx,dy)

 *      dx,dy | Lx.Ly(x+dx,y+dy)  Ly(x+dx,y+dy)   |

 * 

 *  and lambda = 0.06  (totaly empirical :-)

 *  

 *  

 * @author Xavier Philippeau

 */

public class HarrisFastDetectionOperator extends ComponentSeparableOp {

// image info

int width,height;

// precomputed values of the derivatives

private double[][] Lx2,Ly2,Lxy;

// private class that represents a corner

private class Corner {

int x,y;

double measure;

public Corner(int x, int y, double measure) {

this.x = x;

this.y = y;

this.measure = measure;

}

}

// list of corners

private List<List<Corner>> allCorners;

// corner filtering

private int radius = 0;

private double gaussiansigma = 0;

private int minDistance = 0;

private int minMeasure = 0;

public HarrisFastDetectionOperator(int radius, double gaussiansigma, int minMeasure, int minDistance) {

this.radius = radius;

this.gaussiansigma = gaussiansigma;

this.minDistance = minDistance;

this.minMeasure = minMeasure;

}

public void setup(Image in, int canal) {

allCorners = new ArrayList<List<Corner>>();

for(int i=0; i<canal; i++)

allCorners.add(new ArrayList<Corner>());

}

@Override

public void computeComponent(Image out, Image in, int canal) {

filter(out, in, canal);

}

private List<Corner> getCorner(int canal) {

return allCorners.get(canal);

}

/**

* Gaussian function

*/

private double gaussian(double x, double y, double sigma2) {

double t = (x*x+y*y)/(2*sigma2);

double u = 1.0/(2*Math.PI*sigma2);

double e = u*Math.exp( -t );

return e;

}

/**

* Sobel gradient 3x3

*/

private double[] sobel(Image image, int canal, int x, int y) {

int v00=0,v01=0,v02=0,v10=0,v12=0,v20=0,v21=0,v22=0;

int x0 = x-1, x1 = x, x2 = x+1;

int y0 = y-1, y1 = y, y2 = y+1;

if (x0<0) x0=0;

if (y0<0) y0=0;

if (x2>=width) x2=width-1;

if (y2>=height) y2=height-1;

v00=(int)getInsidePixel(image, x0, y0, canal);

v10=(int)getInsidePixel(image, x1, y0, canal);

v20=(int)getInsidePixel(image, x2, y0, canal);

v01=(int)getInsidePixel(image, x0, y1, canal);

v21=(int)getInsidePixel(image, x2, y1, canal);

v02=(int)getInsidePixel(image, x0, y2, canal);

v12=(int)getInsidePixel(image, x1, y2, canal);

v22=(int)getInsidePixel(image, x2, y2, canal);

double sx = (v20+2*v21+v22)-(v00+2*v01+v02);

double sy = (v02+2*v12+v22)-(v00+2*v10+v20);

return new double[] {sx/4,sy/4};

}

 

/**

* Compute the 3 arrays Ix, Iy and Ixy

*/

private void computeDerivatives(Image image, int canal, int radius, double sigma){

this.Lx2 = new double[width][height];

this.Ly2 = new double[width][height];

this.Lxy = new double[width][height];

 

// gradient values: Gx,Gy

double[][][] grad = new double[width][height][];

for (int y=0; y<this.height; y++)

for (int x=0; x<this.width; x++)

grad[x][y]= sobel(image, canal, x, y);

// precompute the coefficients of the gaussian filter 

double[][] filter = new double[2*radius+1][2*radius+1];

double filtersum = 0;

for(int j=-radius;j<=radius;j++) {

for(int i=-radius;i<=radius;i++) {

double g = gaussian(i,j,sigma);

filter[i+radius][j+radius]=g;

filtersum+=g;

}

}

// Convolve gradient with gaussian filter:

//

// Ix2 = (F) * (Gx^2)

// Iy2 = (F) * (Gy^2)

// Ixy = (F) * (Gx.Gy)

//

for (int y=0; y<this.height; y++) {

for (int x=0; x<this.width; x++) {

for(int dy=-radius;dy<=radius;dy++) {

for(int dx=-radius;dx<=radius;dx++) {

int xk = x + dx;

int yk = y + dy;

if (xk<0 || xk>=this.width) continue;

if (yk<0 || yk>=this.height) continue;

// filter weight

double f = filter[dx+radius][dy+radius];

// convolution

this.Lx2[x][y]+=f*grad[xk][yk][0]*grad[xk][yk][0];

this.Ly2[x][y]+=f*grad[xk][yk][1]*grad[xk][yk][1];

this.Lxy[x][y]+=f*grad[xk][yk][0]*grad[xk][yk][1];

}

}

this.Lx2[x][y]/=filtersum;

this.Ly2[x][y]/=filtersum;

this.Lxy[x][y]/=filtersum;

}

}

}

/**

* compute harris measure for a pixel

*/

private double harrisMeasure(int x, int y) {

// matrix elements (normalized)

double m00 = this.Lx2[x][y]; 

double m01 = this.Lxy[x][y];

double m10 = this.Lxy[x][y];

double m11 = this.Ly2[x][y];

// Harris corner measure = det(M)-lambda.trace(M)^2

 

return m00*m11 - m01*m10 - 0.06*(m00+m11)*(m00+m11);

}

/**

* return true if the measure at pixel (x,y) is a local spatial Maxima

*/

private boolean isSpatialMaxima(double[][] hmap, int x, int y) {

int n=8;

int[] dx = new int[] {-1,0,1,1,1,0,-1,-1};

int[] dy = new int[] {-1,-1,-1,0,1,1,1,0};

double w =  hmap[x][y];

for(int i=0;i<n;i++) {

double wk = hmap[x+dx[i]][y+dy[i]];

if (wk>=w) return false;

}

return true;

}

/**

* compute the Harris measure for each pixel of the image

*/

private double[][] computeHarrisMap() {

 

// Harris measure map

double[][] harrismap = new double[width][height];

double max=0;

// for each pixel in the image

for (int y=0; y<this.height; y++) {

for (int x=0; x<this.width; x++) {

// compute ans store the harris measure

harrismap[x][y]=harrisMeasure(x,y);

if (harrismap[x][y]>max) max=harrismap[x][y];

}

}

// rescale measures in 0-100

for (int y=0; y<this.height; y++) {

for (int x=0; x<this.width; x++) {

double h=harrismap[x][y];

if (h<0) h=0; else h = 100 * Math.log(1+h) / Math.log(1+max);

harrismap[x][y]=h;

}

}

return harrismap;

}

/**

* Perfom Harris Corner Detection

*/

public void filter(Image output, Image input, int canal) {

this.width = input.getWidth();

this.height = input.getHeight();

 

// precompute derivatives

computeDerivatives(input, canal, this.radius, this.gaussiansigma);

// Harris measure map

double[][] harrismap = computeHarrisMap();

// copy of the original image (a little darker)

output.resize(this.width, this.height); 

for (int y=0; y<this.height; y++)

for (int x=0; x<this.width; x++)

setInsidePixel(output, x,y, canal, (int)(getInsidePixel(input, x,y, canal)*0.80));

// for each pixel in the hmap, keep the local maxima

for (int y=1; y<this.height-1; y++) {

for (int x=1; x<this.width-1; x++) {

double h = harrismap[x][y];

if (h<this.minMeasure) continue;

if (!isSpatialMaxima(harrismap, (int)x, (int)y)) continue;

// add the corner to the list

getCorner(canal).add( new Corner(x,y,h) );

}

}

// remove corners to close to each other (keep the highest measure)

Iterator<Corner> iter = getCorner(canal).iterator();

while(iter.hasNext()) {

Corner p = iter.next();

for(Corner n:getCorner(canal)) {

if (n==p) continue;

int dist = (int)Math.sqrt( (p.x-n.x)*(p.x-n.x)+(p.y-n.y)*(p.y-n.y) );

if( dist>this.minDistance) continue;

if (n.measure<p.measure) continue;

iter.remove();

break;

}

}

// Display corners over the image (cross)

for (Corner p:getCorner(canal)) {

for (int dx=-2; dx<=2; dx++) {

if (p.x+dx<0 || p.x+dx>=width) continue;

setInsidePixel(output, (int)p.x+dx, (int)p.y, canal, 255);

}

for (int dy=-2; dy<=2; dy++) {

if (p.y+dy<0 || p.y+dy>=height) continue;

setInsidePixel(output, (int)p.x,(int)p.y+dy, canal, 255);

}

}

}

float getInsidePixel(Image c, int x, int y, int canal) {

if(x<0||y<0||x>=c.getWidth()||y>=c.getHeight()) return 0;

return c.getPixel(x, y, canal);

}

void setInsidePixel(Image c, int x, int y, int canal, float value) {

if(x<0||y<0||x>=c.getWidth()||y>=c.getHeight()) return;

c.setPixel(x,y, canal, value);

}

}