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

}