001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.pub.corners;
011
012import java.util.ArrayList;
013import java.util.Collections;
014import java.util.List;
015
016import ij.process.FloatProcessor;
017import ij.process.ImageProcessor;
018import imagingbook.lib.image.Filter;
019import imagingbook.lib.image.ImageMath;
020
021
022/**
023 * This class implements the Harris corner detector.
024 *
025 * @author W. Burger
026 *      @version 2015/07/05
027 */
028public class HarrisCornerDetector {
029        
030
031        /**
032         * Default parameters; a (usually modified) instance of this class
033         * may be passed to constructor of the main class.
034         */
035        public static class Parameters {
036                /** Sensitivity parameter */
037                public double alpha = 0.05;
038                /** Corner response threshold */
039                public double tH = 20000;
040                /** Min. distance between final corners */
041                public double dmin = 10;
042                /** Width of border region ignored in corner search */
043                public int border = 20;
044                /** Set true to perform final corner cleanup */
045                public boolean doCleanUp = true;        
046        }
047        
048        //filter kernels (one-dim. part of separable 2D filters)
049        private final static float[] hp = {2f/9, 5f/9, 2f/9};
050        private final static float[] hd = {0.5f, 0, -0.5f};
051        private final static float[] hb = {1f/64, 6f/64, 15f/64, 20f/64, 15f/64, 6f/64, 1f/64};
052        
053        private final Parameters params;
054        private final int M, N;
055        
056        private FloatProcessor A;                                                       
057        private FloatProcessor B;
058        private FloatProcessor C;
059
060
061        public HarrisCornerDetector(ImageProcessor ip){
062                this(ip, new Parameters());
063        }
064        
065        public HarrisCornerDetector(ImageProcessor ip, Parameters params) {
066                this.M = ip.getWidth();
067                this.N = ip.getHeight();
068                this.params = params;
069                makeDerivatives(ip);
070        }
071        
072        public List<Corner> findCorners() {
073                FloatProcessor Q = makeCrf((float)params.alpha);                                //corner response function (CRF)
074                List<Corner> corners = collectCorners(Q, (float)params.tH, params.border);
075                if (params.doCleanUp) {
076                        corners = cleanupCorners(corners, params.dmin);
077                }
078                return corners;
079        }
080        
081        private void makeDerivatives(ImageProcessor I) {
082                FloatProcessor Ix = I.convertToFloatProcessor(); 
083                FloatProcessor Iy = I.convertToFloatProcessor();
084                
085                Filter.convolveX(Ix, hp);                               // pre-filter Ix horizontally
086                Filter.convolveX(Ix, hd);                               // get horizontal derivative 
087                
088                Filter.convolveY(Iy, hp);                               // pre-filter Iy vertically
089                Filter.convolveY(Iy, hd);                               // get vertical derivative
090                
091                A = ImageMath.sqr(Ix);                                  // A <- Ix^2
092                B = ImageMath.sqr(Iy);                                  // B <- Iy^2
093                C = ImageMath.mult(Ix, Iy);                             // C <- Ix * Iy
094                
095                Filter.convolveXY(A, hb);                               // blur A in x/y
096                Filter.convolveXY(B, hb);                               // blur B in x/y
097                Filter.convolveXY(C, hb);                               // blur C in x/y
098        }
099        
100        private FloatProcessor makeCrf(float alpha) { //corner response function (CRF)
101                FloatProcessor Q = new FloatProcessor(M, N);
102                final float[] pA = (float[]) A.getPixels();
103                final float[] pB = (float[]) B.getPixels();
104                final float[] pC = (float[]) C.getPixels();
105                final float[] pQ = (float[]) Q.getPixels();
106                for (int i = 0; i < M * N; i++) {
107                        float a = pA[i], b = pB[i], c = pC[i];
108                        float det = a * b - c * c;
109                        float trace = a + b;
110                        pQ[i] = det - alpha * (trace * trace);
111                }
112                return Q;
113        }
114        
115        private List<Corner> collectCorners(FloatProcessor Q, float tH, int border) {
116                List<Corner> C = new ArrayList<Corner>();
117                for (int v = border; v < N - border; v++) {
118                        for (int u = border; u < M - border; u++) {
119                                float q = Q.getf(u, v);
120                                if (q > tH && isLocalMax(Q, u, v)) {
121                                        Corner c = new Corner(u, v, q);
122                                        C.add(c);
123                                }
124                        }
125                }
126                return C;
127        }
128        
129        // uses 8-neighborhood
130        private boolean isLocalMax (FloatProcessor Q, int u, int v) {
131                if (u <= 0 || u >= M - 1 || v <= 0 || v >= N - 1) {
132                        return false;
133                } 
134                else {
135                        final float[] q = (float[]) Q.getPixels();
136                        final int i0 = (v - 1) * M + u;
137                        final int i1 = v * M + u;
138                        final int i2 = (v + 1) * M + u;
139                        final float q0 = q[i1];
140                        return  // check 8 neighbors of q0
141                                q0 >= q[i0 - 1] && q0 >= q[i0] && q0 >= q[i0 + 1] &&
142                                q0 >= q[i1 - 1] &&                q0 >= q[i1 + 1] && 
143                                q0 >= q[i2 - 1] && q0 >= q[i2] && q0 >= q[i2 + 1] ;
144                }
145        }
146        
147        private List<Corner> cleanupCorners(List<Corner> C, double dmin){
148                final double dmin2 = dmin * dmin;
149                // sort corners by descending q-value:
150                Collections.sort(C);
151                // we use an array of corners for efficiency reasons:
152                Corner[] Ca = C.toArray(new Corner[C.size()]);
153                List<Corner> Cclean = new ArrayList<Corner>(C.size());
154                for (int i = 0; i < Ca.length; i++) {
155                        Corner c0 = Ca[i];              // get next strongest corner
156                        if (c0 != null) {
157                                Cclean.add(c0);
158                                // delete all remaining corners cj too close to c0:
159                                for (int j = i + 1; j < Ca.length; j++) {
160                                        Corner cj = Ca[j];
161                                        if (cj != null && c0.dist2(cj) < dmin2)
162                                                Ca[j] = null;   //delete corner cj from C
163                                }
164                        }
165                }
166                return Cclean;
167        }
168
169}